[filterbank] check samplerate in _slaney, use temp variables
[aubio.git] / src / spectral / filterbank_mel.c
1 /*
2   Copyright (C) 2007-2009 Paul Brossier <piem@aubio.org>
3                       and Amaury Hazan <ahazan@iua.upf.edu>
4
5   This file is part of aubio.
6
7   aubio is free software: you can redistribute it and/or modify
8   it under the terms of the GNU General Public License as published by
9   the Free Software Foundation, either version 3 of the License, or
10   (at your option) any later version.
11
12   aubio is distributed in the hope that it will be useful,
13   but WITHOUT ANY WARRANTY; without even the implied warranty of
14   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15   GNU General Public License for more details.
16
17   You should have received a copy of the GNU General Public License
18   along with aubio.  If not, see <http://www.gnu.org/licenses/>.
19
20 */
21
22 #include "aubio_priv.h"
23 #include "fmat.h"
24 #include "fvec.h"
25 #include "cvec.h"
26 #include "spectral/filterbank.h"
27 #include "spectral/filterbank_mel.h"
28 #include "mathutils.h"
29
30 uint_t
31 aubio_filterbank_set_triangle_bands (aubio_filterbank_t * fb,
32     const fvec_t * freqs, smpl_t samplerate)
33 {
34
35   fmat_t *filters = aubio_filterbank_get_coeffs (fb);
36   uint_t n_filters = filters->height, win_s = filters->length;
37   fvec_t *lower_freqs, *upper_freqs, *center_freqs;
38   fvec_t *triangle_heights, *fft_freqs;
39
40   uint_t fn;                    /* filter counter */
41   uint_t bin;                   /* bin counter */
42
43   smpl_t riseInc, downInc;
44
45   /* freqs define the bands of triangular overlapping windows.
46      throw a warning if filterbank object fb is too short. */
47   if (freqs->length - 2 > n_filters) {
48     AUBIO_WRN ("not enough filters, %d allocated but %d requested\n",
49         n_filters, freqs->length - 2);
50   }
51
52   if (freqs->length - 2 < n_filters) {
53     AUBIO_WRN ("too many filters, %d allocated but %d requested\n",
54         n_filters, freqs->length - 2);
55   }
56
57   for (fn = 0; fn < freqs->length; fn++) {
58     if (freqs->data[fn] < 0) {
59       AUBIO_ERR("filterbank_mel: freqs must contain only positive values.\n");
60       return AUBIO_FAIL;
61     } else if (freqs->data[fn] > samplerate / 2) {
62       AUBIO_WRN("filterbank_mel: freqs should contain only "
63           "values < samplerate / 2.\n");
64     } else if (fn > 0 && freqs->data[fn] < freqs->data[fn-1]) {
65       AUBIO_ERR("filterbank_mel: freqs should be a list of frequencies "
66           "sorted from low to high, but freq[%d] < freq[%d-1]\n", fn, fn);
67       return AUBIO_FAIL;
68     } else if (fn > 0 && freqs->data[fn] == freqs->data[fn-1]) {
69       AUBIO_WRN("filterbank_mel: set_triangle_bands received a list "
70           "with twice the frequency %f\n", freqs->data[fn]);
71     }
72   }
73
74   /* convenience reference to lower/center/upper frequency for each triangle */
75   lower_freqs = new_fvec (n_filters);
76   upper_freqs = new_fvec (n_filters);
77   center_freqs = new_fvec (n_filters);
78
79   /* height of each triangle */
80   triangle_heights = new_fvec (n_filters);
81
82   /* lookup table of each bin frequency in hz */
83   fft_freqs = new_fvec (win_s);
84
85   /* fill up the lower/center/upper */
86   for (fn = 0; fn < n_filters; fn++) {
87     lower_freqs->data[fn] = freqs->data[fn];
88     center_freqs->data[fn] = freqs->data[fn + 1];
89     upper_freqs->data[fn] = freqs->data[fn + 2];
90   }
91
92   /* compute triangle heights so that each triangle has unit area */
93   if (aubio_filterbank_get_norm(fb)) {
94     for (fn = 0; fn < n_filters; fn++) {
95       triangle_heights->data[fn] =
96           2. / (upper_freqs->data[fn] - lower_freqs->data[fn]);
97     }
98   } else {
99     fvec_ones (triangle_heights);
100   }
101
102   /* fill fft_freqs lookup table, which assigns the frequency in hz to each bin */
103   for (bin = 0; bin < win_s; bin++) {
104     fft_freqs->data[bin] =
105         aubio_bintofreq (bin, samplerate, (win_s - 1) * 2);
106   }
107
108   /* zeroing of all filters */
109   fmat_zeros (filters);
110
111   /* building each filter table */
112   for (fn = 0; fn < n_filters; fn++) {
113
114     /* skip first elements */
115     for (bin = 0; bin < win_s - 1; bin++) {
116       if (fft_freqs->data[bin] <= lower_freqs->data[fn] &&
117           fft_freqs->data[bin + 1] > lower_freqs->data[fn]) {
118         bin++;
119         break;
120       }
121     }
122
123     /* compute positive slope step size */
124     riseInc = triangle_heights->data[fn]
125       / (center_freqs->data[fn] - lower_freqs->data[fn]);
126
127     /* compute coefficients in positive slope */
128     for (; bin < win_s - 1; bin++) {
129       filters->data[fn][bin] =
130           (fft_freqs->data[bin] - lower_freqs->data[fn]) * riseInc;
131
132       if (fft_freqs->data[bin + 1] >= center_freqs->data[fn]) {
133         bin++;
134         break;
135       }
136     }
137
138     /* compute negative slope step size */
139     downInc = triangle_heights->data[fn]
140       / (upper_freqs->data[fn] - center_freqs->data[fn]);
141
142     /* compute coefficents in negative slope */
143     for (; bin < win_s - 1; bin++) {
144       filters->data[fn][bin] +=
145           (upper_freqs->data[fn] - fft_freqs->data[bin]) * downInc;
146
147       if (filters->data[fn][bin] < 0.) {
148         filters->data[fn][bin] = 0.;
149       }
150
151       if (fft_freqs->data[bin + 1] >= upper_freqs->data[fn])
152         break;
153     }
154     /* nothing else to do */
155
156   }
157
158   /* destroy temporarly allocated vectors */
159   del_fvec (lower_freqs);
160   del_fvec (upper_freqs);
161   del_fvec (center_freqs);
162
163   del_fvec (triangle_heights);
164   del_fvec (fft_freqs);
165
166   return AUBIO_OK;
167 }
168
169 uint_t
170 aubio_filterbank_set_mel_coeffs_slaney (aubio_filterbank_t * fb,
171     smpl_t samplerate)
172 {
173   uint_t retval;
174
175   if (samplerate <= 0) {
176     AUBIO_ERR("filterbank: set_mel_coeffs_slaney samplerate should be > 0\n");
177     return AUBIO_FAIL;
178   }
179
180   /* Malcolm Slaney parameters */
181   smpl_t lowestFrequency = 133.3333;
182   smpl_t linearSpacing = 66.66666666;
183   smpl_t logSpacing = 1.0711703;
184
185   uint_t linearFilters = 13;
186   uint_t logFilters = 27;
187   uint_t n_filters = linearFilters + logFilters;
188
189   uint_t fn;                    /* filter counter */
190
191   smpl_t lastlinearCF;
192
193   /* buffers to compute filter frequencies */
194   fvec_t *freqs = new_fvec (n_filters + 2);
195
196   /* first step: fill all the linear filter frequencies */
197   for (fn = 0; fn < linearFilters; fn++) {
198     freqs->data[fn] = lowestFrequency + fn * linearSpacing;
199   }
200   lastlinearCF = freqs->data[fn - 1];
201
202   /* second step: fill all the log filter frequencies */
203   for (fn = 0; fn < logFilters + 2; fn++) {
204     freqs->data[fn + linearFilters] =
205         lastlinearCF * (POW (logSpacing, fn + 1));
206   }
207
208   /* now compute the actual coefficients */
209   retval = aubio_filterbank_set_triangle_bands (fb, freqs, samplerate);
210
211   /* destroy vector used to store frequency limits */
212   del_fvec (freqs);
213
214   return retval;
215 }
216
217 static uint_t aubio_filterbank_check_freqs (aubio_filterbank_t *fb UNUSED,
218     smpl_t samplerate, smpl_t *freq_min, smpl_t *freq_max)
219 {
220   if (samplerate <= 0) {
221     AUBIO_ERR("filterbank: set_mel_coeffs samplerate should be > 0\n");
222     return AUBIO_FAIL;
223   }
224   if (*freq_max < 0) {
225     AUBIO_ERR("filterbank: set_mel_coeffs freq_max should be > 0\n");
226     return AUBIO_FAIL;
227   } else if (*freq_max == 0) {
228     *freq_max = samplerate / 2.;
229   }
230   if (*freq_min < 0) {
231     AUBIO_ERR("filterbank: set_mel_coeffs freq_min should be > 0\n");
232     return AUBIO_FAIL;
233   }
234   return AUBIO_OK;
235 }
236
237 uint_t
238 aubio_filterbank_set_mel_coeffs (aubio_filterbank_t * fb, smpl_t samplerate,
239     smpl_t freq_min, smpl_t freq_max)
240 {
241   uint_t m, retval;
242   smpl_t start = freq_min, end = freq_max, step;
243   fvec_t *freqs;
244   fmat_t *coeffs = aubio_filterbank_get_coeffs(fb);
245   uint_t n_bands = coeffs->height;
246
247   if (aubio_filterbank_check_freqs(fb, samplerate, &start, &end)) {
248     return AUBIO_FAIL;
249   }
250
251   start = aubio_hztomel(start);
252   end = aubio_hztomel(end);
253
254   freqs = new_fvec(n_bands + 2);
255   step = (end - start) / (n_bands + 1);
256
257   for (m = 0; m < n_bands + 2; m++)
258   {
259     freqs->data[m] = MIN(aubio_meltohz(start + step * m), samplerate/2.);
260   }
261
262   retval = aubio_filterbank_set_triangle_bands (fb, freqs, samplerate);
263
264   /* destroy vector used to store frequency limits */
265   del_fvec (freqs);
266   return retval;
267 }
268
269 uint_t
270 aubio_filterbank_set_mel_coeffs_htk (aubio_filterbank_t * fb, smpl_t samplerate,
271     smpl_t freq_min, smpl_t freq_max)
272 {
273   uint_t m, retval;
274   smpl_t start = freq_min, end = freq_max, step;
275   fvec_t *freqs;
276   fmat_t *coeffs = aubio_filterbank_get_coeffs(fb);
277   uint_t n_bands = coeffs->height;
278
279   if (aubio_filterbank_check_freqs(fb, samplerate, &start, &end)) {
280     return AUBIO_FAIL;
281   }
282
283   start = aubio_hztomel_htk(start);
284   end = aubio_hztomel_htk(end);
285
286   freqs = new_fvec (n_bands + 2);
287   step = (end - start) / (n_bands + 1);
288
289   for (m = 0; m < n_bands + 2; m++)
290   {
291     freqs->data[m] = MIN(aubio_meltohz_htk(start + step * m), samplerate/2.);
292   }
293
294   retval = aubio_filterbank_set_triangle_bands (fb, freqs, samplerate);
295
296   /* destroy vector used to store frequency limits */
297   del_fvec (freqs);
298   return retval;
299 }