Merge branch 'master' into awhitening
[aubio.git] / src / spectral / specdesc.c
1 /*
2   Copyright (C) 2003-2009 Paul Brossier <piem@aubio.org>
3
4   This file is part of aubio.
5
6   aubio is free software: you can redistribute it and/or modify
7   it under the terms of the GNU General Public License as published by
8   the Free Software Foundation, either version 3 of the License, or
9   (at your option) any later version.
10
11   aubio is distributed in the hope that it will be useful,
12   but WITHOUT ANY WARRANTY; without even the implied warranty of
13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14   GNU General Public License for more details.
15
16   You should have received a copy of the GNU General Public License
17   along with aubio.  If not, see <http://www.gnu.org/licenses/>.
18
19 */
20
21 #include "aubio_priv.h"
22 #include "fvec.h"
23 #include "cvec.h"
24 #include "spectral/fft.h"
25 #include "spectral/specdesc.h"
26 #include "mathutils.h"
27 #include "utils/hist.h"
28
29 void aubio_specdesc_energy(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
30 void aubio_specdesc_hfc(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
31 void aubio_specdesc_complex(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
32 void aubio_specdesc_phase(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
33 void aubio_specdesc_wphase(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
34 void aubio_specdesc_specdiff(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
35 void aubio_specdesc_kl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
36 void aubio_specdesc_mkl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
37 void aubio_specdesc_specflux(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
38
39 extern void aubio_specdesc_centroid (aubio_specdesc_t * o, const cvec_t * spec,
40     fvec_t * desc);
41 extern void aubio_specdesc_spread (aubio_specdesc_t * o, const cvec_t * spec,
42     fvec_t * desc);
43 extern void aubio_specdesc_skewness (aubio_specdesc_t * o, const cvec_t * spec,
44     fvec_t * desc);
45 extern void aubio_specdesc_kurtosis (aubio_specdesc_t * o, const cvec_t * spec,
46     fvec_t * desc);
47 extern void aubio_specdesc_slope (aubio_specdesc_t * o, const cvec_t * spec,
48     fvec_t * desc);
49 extern void aubio_specdesc_decrease (aubio_specdesc_t * o, const cvec_t * spec,
50     fvec_t * desc);
51 extern void aubio_specdesc_rolloff (aubio_specdesc_t * o, const cvec_t * spec,
52     fvec_t * desc);
53
54 /** onsetdetection types */
55 typedef enum {
56         aubio_onset_energy,         /**< energy based */          
57         aubio_onset_specdiff,       /**< spectral diff */         
58         aubio_onset_hfc,            /**< high frequency content */
59         aubio_onset_complex,        /**< complex domain */        
60         aubio_onset_phase,          /**< phase fast */            
61         aubio_onset_wphase,         /**< weighted phase */
62         aubio_onset_kl,             /**< Kullback Liebler */
63         aubio_onset_mkl,            /**< modified Kullback Liebler */
64         aubio_onset_specflux,       /**< spectral flux */
65         aubio_specmethod_centroid,  /**< spectral centroid */
66         aubio_specmethod_spread,    /**< spectral spread */
67         aubio_specmethod_skewness,  /**< spectral skewness */
68         aubio_specmethod_kurtosis,  /**< spectral kurtosis */
69         aubio_specmethod_slope,     /**< spectral kurtosis */
70         aubio_specmethod_decrease,  /**< spectral decrease */
71         aubio_specmethod_rolloff,   /**< spectral rolloff */
72         aubio_onset_default = aubio_onset_hfc, /**< default mode, set to hfc */
73 } aubio_specdesc_type;
74
75 /** structure to store object state */
76 struct _aubio_specdesc_t {
77   aubio_specdesc_type onset_type; /**< onset detection type */
78   /** Pointer to aubio_specdesc_<type> function */
79   void (*funcpointer)(aubio_specdesc_t *o,
80       const cvec_t * fftgrain, fvec_t * onset);
81   smpl_t threshold;      /**< minimum norm threshold for phase and specdiff */
82   fvec_t *oldmag;        /**< previous norm vector */
83   fvec_t *dev1 ;         /**< current onset detection measure vector */
84   fvec_t *theta1;        /**< previous phase vector, one frame behind */
85   fvec_t *theta2;        /**< previous phase vector, two frames behind */
86   aubio_hist_t * histog; /**< histogram */
87 };
88
89
90 /* Energy based onset detection function */
91 void aubio_specdesc_energy  (aubio_specdesc_t *o UNUSED,
92     const cvec_t * fftgrain, fvec_t * onset) {
93   uint_t j;
94   onset->data[0] = 0.;
95   for (j=0;j<fftgrain->length;j++) {
96     onset->data[0] += SQR(fftgrain->norm[j]);
97   }
98 }
99
100 /* High Frequency Content onset detection function */
101 void aubio_specdesc_hfc(aubio_specdesc_t *o UNUSED,
102     const cvec_t * fftgrain, fvec_t * onset){
103   uint_t j;
104   onset->data[0] = 0.;
105   for (j=0;j<fftgrain->length;j++) {
106     onset->data[0] += (j+1)*fftgrain->norm[j];
107   }
108 }
109
110
111 /* Complex Domain Method onset detection function */
112 void aubio_specdesc_complex (aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset) {
113   uint_t j;
114   uint_t nbins = fftgrain->length;
115   onset->data[0] = 0.;
116   for (j=0;j<nbins; j++)  {
117     // compute the predicted phase
118     o->dev1->data[j] = 2. * o->theta1->data[j] - o->theta2->data[j];
119     // compute the euclidean distance in the complex domain
120     // sqrt ( r_1^2 + r_2^2 - 2 * r_1 * r_2 * \cos ( \phi_1 - \phi_2 ) )
121     onset->data[0] +=
122       SQRT (ABS (SQR (o->oldmag->data[j]) + SQR (fftgrain->norm[j])
123             - 2 * o->oldmag->data[j] * fftgrain->norm[j]
124             * COS (o->dev1->data[j] - fftgrain->phas[j])));
125     /* swap old phase data (need to remember 2 frames behind)*/
126     o->theta2->data[j] = o->theta1->data[j];
127     o->theta1->data[j] = fftgrain->phas[j];
128     /* swap old magnitude data (1 frame is enough) */
129     o->oldmag->data[j] = fftgrain->norm[j];
130   }
131 }
132
133
134 /* Phase Based Method onset detection function */
135 void aubio_specdesc_phase(aubio_specdesc_t *o, 
136     const cvec_t * fftgrain, fvec_t * onset){
137   uint_t j;
138   uint_t nbins = fftgrain->length;
139   onset->data[0] = 0.0;
140   o->dev1->data[0]=0.;
141   for ( j=0;j<nbins; j++ )  {
142     o->dev1->data[j] = 
143       aubio_unwrap2pi(
144           fftgrain->phas[j]
145           -2.0*o->theta1->data[j]
146           +o->theta2->data[j]);
147     if ( o->threshold < fftgrain->norm[j] )
148       o->dev1->data[j] = ABS(o->dev1->data[j]);
149     else 
150       o->dev1->data[j] = 0.0;
151     /* keep a track of the past frames */
152     o->theta2->data[j] = o->theta1->data[j];
153     o->theta1->data[j] = fftgrain->phas[j];
154   }
155   /* apply o->histogram */
156   aubio_hist_dyn_notnull(o->histog,o->dev1);
157   /* weight it */
158   aubio_hist_weight(o->histog);
159   /* its mean is the result */
160   onset->data[0] = aubio_hist_mean(o->histog);  
161   //onset->data[0] = fvec_mean(o->dev1);
162 }
163
164 /* weighted phase */
165 void
166 aubio_specdesc_wphase(aubio_specdesc_t *o,
167     const cvec_t *fftgrain, fvec_t *onset) {
168   uint_t i;
169   aubio_specdesc_phase(o, fftgrain, onset);
170   for (i = 0; i < fftgrain->length; i++) {
171     o->dev1->data[i] *= fftgrain->norm[i];
172   }
173   /* apply o->histogram */
174   aubio_hist_dyn_notnull(o->histog,o->dev1);
175   /* weight it */
176   aubio_hist_weight(o->histog);
177   /* its mean is the result */
178   onset->data[0] = aubio_hist_mean(o->histog);
179 }
180
181 /* Spectral difference method onset detection function */
182 void aubio_specdesc_specdiff(aubio_specdesc_t *o,
183     const cvec_t * fftgrain, fvec_t * onset){
184   uint_t j;
185   uint_t nbins = fftgrain->length;
186     onset->data[0] = 0.0;
187     for (j=0;j<nbins; j++)  {
188       o->dev1->data[j] = SQRT(
189           ABS(SQR( fftgrain->norm[j])
190             - SQR(o->oldmag->data[j])));
191       if (o->threshold < fftgrain->norm[j] )
192         o->dev1->data[j] = ABS(o->dev1->data[j]);
193       else 
194         o->dev1->data[j] = 0.0;
195       o->oldmag->data[j] = fftgrain->norm[j];
196     }
197
198     /* apply o->histogram (act somewhat as a low pass on the
199      * overall function)*/
200     aubio_hist_dyn_notnull(o->histog,o->dev1);
201     /* weight it */
202     aubio_hist_weight(o->histog);
203     /* its mean is the result */
204     onset->data[0] = aubio_hist_mean(o->histog);  
205 }
206
207 /* Kullback Liebler onset detection function
208  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid 
209  * negative (1.+) and infinite values (+1.e-10) */
210 void aubio_specdesc_kl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){
211   uint_t j;
212     onset->data[0] = 0.;
213     for (j=0;j<fftgrain->length;j++) {
214       onset->data[0] += fftgrain->norm[j]
215         *LOG(1.+fftgrain->norm[j]/(o->oldmag->data[j]+1.e-1));
216       o->oldmag->data[j] = fftgrain->norm[j];
217     }
218     if (isnan(onset->data[0])) onset->data[0] = 0.;
219 }
220
221 /* Modified Kullback Liebler onset detection function
222  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid 
223  * negative (1.+) and infinite values (+1.e-10) */
224 void aubio_specdesc_mkl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){
225   uint_t j;
226     onset->data[0] = 0.;
227     for (j=0;j<fftgrain->length;j++) {
228       onset->data[0] += LOG(1.+fftgrain->norm[j]/(o->oldmag->data[j]+1.e-1));
229       o->oldmag->data[j] = fftgrain->norm[j];
230     }
231     if (isnan(onset->data[0])) onset->data[0] = 0.;
232 }
233
234 /* Spectral flux */
235 void aubio_specdesc_specflux(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){ 
236   uint_t j;
237   onset->data[0] = 0.;
238   for (j=0;j<fftgrain->length;j++) {
239     if (fftgrain->norm[j] > o->oldmag->data[j])
240       onset->data[0] += fftgrain->norm[j] - o->oldmag->data[j];
241     o->oldmag->data[j] = fftgrain->norm[j];
242   }
243 }
244
245 /* Generic function pointing to the choosen one */
246 void 
247 aubio_specdesc_do (aubio_specdesc_t *o, const cvec_t * fftgrain, 
248     fvec_t * onset) {
249   o->funcpointer(o,fftgrain,onset);
250 }
251
252 /* Allocate memory for an onset detection 
253  * depending on the choosen type, allocate memory as needed
254  */
255 aubio_specdesc_t * 
256 new_aubio_specdesc (const char_t * onset_mode, uint_t size){
257   aubio_specdesc_t * o = AUBIO_NEW(aubio_specdesc_t);
258   uint_t rsize = size/2+1;
259   aubio_specdesc_type onset_type;
260   if (strcmp (onset_mode, "energy") == 0)
261       onset_type = aubio_onset_energy;
262   else if (strcmp (onset_mode, "specdiff") == 0)
263       onset_type = aubio_onset_specdiff;
264   else if (strcmp (onset_mode, "hfc") == 0)
265       onset_type = aubio_onset_hfc;
266   else if (strcmp (onset_mode, "complexdomain") == 0)
267       onset_type = aubio_onset_complex;
268   else if (strcmp (onset_mode, "complex") == 0)
269       onset_type = aubio_onset_complex;
270   else if (strcmp (onset_mode, "phase") == 0)
271       onset_type = aubio_onset_phase;
272   else if (strcmp (onset_mode, "wphase") == 0)
273       onset_type = aubio_onset_wphase;
274   else if (strcmp (onset_mode, "mkl") == 0)
275       onset_type = aubio_onset_mkl;
276   else if (strcmp (onset_mode, "kl") == 0)
277       onset_type = aubio_onset_kl;
278   else if (strcmp (onset_mode, "specflux") == 0)
279       onset_type = aubio_onset_specflux;
280   else if (strcmp (onset_mode, "centroid") == 0)
281       onset_type = aubio_specmethod_centroid;
282   else if (strcmp (onset_mode, "spread") == 0)
283       onset_type = aubio_specmethod_spread;
284   else if (strcmp (onset_mode, "skewness") == 0)
285       onset_type = aubio_specmethod_skewness;
286   else if (strcmp (onset_mode, "kurtosis") == 0)
287       onset_type = aubio_specmethod_kurtosis;
288   else if (strcmp (onset_mode, "slope") == 0)
289       onset_type = aubio_specmethod_slope;
290   else if (strcmp (onset_mode, "decrease") == 0)
291       onset_type = aubio_specmethod_decrease;
292   else if (strcmp (onset_mode, "rolloff") == 0)
293       onset_type = aubio_specmethod_rolloff;
294   else if (strcmp (onset_mode, "default") == 0)
295       onset_type = aubio_onset_default;
296   else {
297       AUBIO_ERR("unknown spectral descriptor type %s\n", onset_mode);
298       AUBIO_FREE(o);
299       return NULL;
300   }
301   switch(onset_type) {
302     /* for both energy and hfc, only fftgrain->norm is required */
303     case aubio_onset_energy:
304       break;
305     case aubio_onset_hfc:
306       break;
307       /* the other approaches will need some more memory spaces */
308     case aubio_onset_complex:
309       o->oldmag = new_fvec(rsize);
310       o->dev1   = new_fvec(rsize);
311       o->theta1 = new_fvec(rsize);
312       o->theta2 = new_fvec(rsize);
313       break;
314     case aubio_onset_phase:
315     case aubio_onset_wphase:
316       o->dev1   = new_fvec(rsize);
317       o->theta1 = new_fvec(rsize);
318       o->theta2 = new_fvec(rsize);
319       o->histog = new_aubio_hist(0.0, PI, 10);
320       o->threshold = 0.1;
321       break;
322     case aubio_onset_specdiff:
323       o->oldmag = new_fvec(rsize);
324       o->dev1   = new_fvec(rsize);
325       o->histog = new_aubio_hist(0.0, PI, 10);
326       o->threshold = 0.1;
327       break;
328     case aubio_onset_kl:
329     case aubio_onset_mkl:
330     case aubio_onset_specflux:
331       o->oldmag = new_fvec(rsize);
332       break;
333     default:
334       break;
335   }
336
337   switch(onset_type) {
338     case aubio_onset_energy:
339       o->funcpointer = aubio_specdesc_energy;
340       break;
341     case aubio_onset_hfc:
342       o->funcpointer = aubio_specdesc_hfc;
343       break;
344     case aubio_onset_complex:
345       o->funcpointer = aubio_specdesc_complex;
346       break;
347     case aubio_onset_phase:
348       o->funcpointer = aubio_specdesc_phase;
349       break;
350     case aubio_onset_wphase:
351       o->funcpointer = aubio_specdesc_wphase;
352       break;
353     case aubio_onset_specdiff:
354       o->funcpointer = aubio_specdesc_specdiff;
355       break;
356     case aubio_onset_kl:
357       o->funcpointer = aubio_specdesc_kl;
358       break;
359     case aubio_onset_mkl:
360       o->funcpointer = aubio_specdesc_mkl;
361       break;
362     case aubio_onset_specflux:
363       o->funcpointer = aubio_specdesc_specflux;
364       break;
365     case aubio_specmethod_centroid:
366       o->funcpointer = aubio_specdesc_centroid;
367       break;
368     case aubio_specmethod_spread:
369       o->funcpointer = aubio_specdesc_spread;
370       break;
371     case aubio_specmethod_skewness:
372       o->funcpointer = aubio_specdesc_skewness;
373       break;
374     case aubio_specmethod_kurtosis:
375       o->funcpointer = aubio_specdesc_kurtosis;
376       break;
377     case aubio_specmethod_slope:
378       o->funcpointer = aubio_specdesc_slope;
379       break;
380     case aubio_specmethod_decrease:
381       o->funcpointer = aubio_specdesc_decrease;
382       break;
383     case aubio_specmethod_rolloff:
384       o->funcpointer = aubio_specdesc_rolloff;
385       break;
386     default:
387       break;
388   }
389   o->onset_type = onset_type;
390   return o;
391 }
392
393 void del_aubio_specdesc (aubio_specdesc_t *o){
394   switch(o->onset_type) {
395     case aubio_onset_energy:
396       break;
397     case aubio_onset_hfc:
398       break;
399     case aubio_onset_complex:
400       del_fvec(o->oldmag);
401       del_fvec(o->dev1);
402       del_fvec(o->theta1);
403       del_fvec(o->theta2);
404       break;
405     case aubio_onset_phase:
406     case aubio_onset_wphase:
407       del_fvec(o->dev1);
408       del_fvec(o->theta1);
409       del_fvec(o->theta2);
410       del_aubio_hist(o->histog);
411       break;
412     case aubio_onset_specdiff:
413       del_fvec(o->oldmag);
414       del_fvec(o->dev1);
415       del_aubio_hist(o->histog);
416       break;
417     case aubio_onset_kl:
418     case aubio_onset_mkl:
419     case aubio_onset_specflux:
420       del_fvec(o->oldmag);
421       break;
422     default:
423       break;
424   }
425   AUBIO_FREE(o);
426 }