README.md: fix appveyor link
[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_specdiff(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
34 void aubio_specdesc_kl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
35 void aubio_specdesc_mkl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
36 void aubio_specdesc_specflux(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
37
38 extern void aubio_specdesc_centroid (aubio_specdesc_t * o, const cvec_t * spec,
39     fvec_t * desc);
40 extern void aubio_specdesc_spread (aubio_specdesc_t * o, const cvec_t * spec,
41     fvec_t * desc);
42 extern void aubio_specdesc_skewness (aubio_specdesc_t * o, const cvec_t * spec,
43     fvec_t * desc);
44 extern void aubio_specdesc_kurtosis (aubio_specdesc_t * o, const cvec_t * spec,
45     fvec_t * desc);
46 extern void aubio_specdesc_slope (aubio_specdesc_t * o, const cvec_t * spec,
47     fvec_t * desc);
48 extern void aubio_specdesc_decrease (aubio_specdesc_t * o, const cvec_t * spec,
49     fvec_t * desc);
50 extern void aubio_specdesc_rolloff (aubio_specdesc_t * o, const cvec_t * spec,
51     fvec_t * desc);
52
53 /** onsetdetection types */
54 typedef enum {
55         aubio_onset_energy,         /**< energy based */          
56         aubio_onset_specdiff,       /**< spectral diff */         
57         aubio_onset_hfc,            /**< high frequency content */
58         aubio_onset_complex,        /**< complex domain */        
59         aubio_onset_phase,          /**< phase fast */            
60         aubio_onset_kl,             /**< Kullback Liebler */
61         aubio_onset_mkl,            /**< modified Kullback Liebler */
62         aubio_onset_specflux,       /**< spectral flux */
63         aubio_specmethod_centroid,  /**< spectral centroid */
64         aubio_specmethod_spread,    /**< spectral spread */
65         aubio_specmethod_skewness,  /**< spectral skewness */
66         aubio_specmethod_kurtosis,  /**< spectral kurtosis */
67         aubio_specmethod_slope,     /**< spectral kurtosis */
68         aubio_specmethod_decrease,  /**< spectral decrease */
69         aubio_specmethod_rolloff,   /**< spectral rolloff */
70         aubio_onset_default = aubio_onset_hfc, /**< default mode, set to hfc */
71 } aubio_specdesc_type;
72
73 /** structure to store object state */
74 struct _aubio_specdesc_t {
75   aubio_specdesc_type onset_type; /**< onset detection type */
76   /** Pointer to aubio_specdesc_<type> function */
77   void (*funcpointer)(aubio_specdesc_t *o,
78       const cvec_t * fftgrain, fvec_t * onset);
79   smpl_t threshold;      /**< minimum norm threshold for phase and specdiff */
80   fvec_t *oldmag;        /**< previous norm vector */
81   fvec_t *dev1 ;         /**< current onset detection measure vector */
82   fvec_t *theta1;        /**< previous phase vector, one frame behind */
83   fvec_t *theta2;        /**< previous phase vector, two frames behind */
84   aubio_hist_t * histog; /**< histogram */
85 };
86
87
88 /* Energy based onset detection function */
89 void aubio_specdesc_energy  (aubio_specdesc_t *o UNUSED,
90     const cvec_t * fftgrain, fvec_t * onset) {
91   uint_t j;
92   onset->data[0] = 0.;
93   for (j=0;j<fftgrain->length;j++) {
94     onset->data[0] += SQR(fftgrain->norm[j]);
95   }
96 }
97
98 /* High Frequency Content onset detection function */
99 void aubio_specdesc_hfc(aubio_specdesc_t *o UNUSED,
100     const cvec_t * fftgrain, fvec_t * onset){
101   uint_t j;
102   onset->data[0] = 0.;
103   for (j=0;j<fftgrain->length;j++) {
104     onset->data[0] += (j+1)*fftgrain->norm[j];
105   }
106 }
107
108
109 /* Complex Domain Method onset detection function */
110 void aubio_specdesc_complex (aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset) {
111   uint_t j;
112   uint_t nbins = fftgrain->length;
113   onset->data[0] = 0.;
114   for (j=0;j<nbins; j++)  {
115     // compute the predicted phase
116     o->dev1->data[j] = 2. * o->theta1->data[j] - o->theta2->data[j];
117     // compute the euclidean distance in the complex domain
118     // sqrt ( r_1^2 + r_2^2 - 2 * r_1 * r_2 * \cos ( \phi_1 - \phi_2 ) )
119     onset->data[0] +=
120       SQRT (ABS (SQR (o->oldmag->data[j]) + SQR (fftgrain->norm[j])
121             - 2 * o->oldmag->data[j] * fftgrain->norm[j]
122             * COS (o->dev1->data[j] - fftgrain->phas[j])));
123     /* swap old phase data (need to remember 2 frames behind)*/
124     o->theta2->data[j] = o->theta1->data[j];
125     o->theta1->data[j] = fftgrain->phas[j];
126     /* swap old magnitude data (1 frame is enough) */
127     o->oldmag->data[j] = fftgrain->norm[j];
128   }
129 }
130
131
132 /* Phase Based Method onset detection function */
133 void aubio_specdesc_phase(aubio_specdesc_t *o, 
134     const cvec_t * fftgrain, fvec_t * onset){
135   uint_t j;
136   uint_t nbins = fftgrain->length;
137   onset->data[0] = 0.0;
138   o->dev1->data[0]=0.;
139   for ( j=0;j<nbins; j++ )  {
140     o->dev1->data[j] = 
141       aubio_unwrap2pi(
142           fftgrain->phas[j]
143           -2.0*o->theta1->data[j]
144           +o->theta2->data[j]);
145     if ( o->threshold < fftgrain->norm[j] )
146       o->dev1->data[j] = ABS(o->dev1->data[j]);
147     else 
148       o->dev1->data[j] = 0.0;
149     /* keep a track of the past frames */
150     o->theta2->data[j] = o->theta1->data[j];
151     o->theta1->data[j] = fftgrain->phas[j];
152   }
153   /* apply o->histogram */
154   aubio_hist_dyn_notnull(o->histog,o->dev1);
155   /* weight it */
156   aubio_hist_weight(o->histog);
157   /* its mean is the result */
158   onset->data[0] = aubio_hist_mean(o->histog);  
159   //onset->data[0] = fvec_mean(o->dev1);
160 }
161
162 /* Spectral difference method onset detection function */
163 void aubio_specdesc_specdiff(aubio_specdesc_t *o,
164     const cvec_t * fftgrain, fvec_t * onset){
165   uint_t j;
166   uint_t nbins = fftgrain->length;
167     onset->data[0] = 0.0;
168     for (j=0;j<nbins; j++)  {
169       o->dev1->data[j] = SQRT(
170           ABS(SQR( fftgrain->norm[j])
171             - SQR(o->oldmag->data[j])));
172       if (o->threshold < fftgrain->norm[j] )
173         o->dev1->data[j] = ABS(o->dev1->data[j]);
174       else 
175         o->dev1->data[j] = 0.0;
176       o->oldmag->data[j] = fftgrain->norm[j];
177     }
178
179     /* apply o->histogram (act somewhat as a low pass on the
180      * overall function)*/
181     aubio_hist_dyn_notnull(o->histog,o->dev1);
182     /* weight it */
183     aubio_hist_weight(o->histog);
184     /* its mean is the result */
185     onset->data[0] = aubio_hist_mean(o->histog);  
186 }
187
188 /* Kullback Liebler onset detection function
189  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid 
190  * negative (1.+) and infinite values (+1.e-10) */
191 void aubio_specdesc_kl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){
192   uint_t j;
193     onset->data[0] = 0.;
194     for (j=0;j<fftgrain->length;j++) {
195       onset->data[0] += fftgrain->norm[j]
196         *LOG(1.+fftgrain->norm[j]/(o->oldmag->data[j]+1.e-1));
197       o->oldmag->data[j] = fftgrain->norm[j];
198     }
199     if (isnan(onset->data[0])) onset->data[0] = 0.;
200 }
201
202 /* Modified Kullback Liebler onset detection function
203  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid 
204  * negative (1.+) and infinite values (+1.e-10) */
205 void aubio_specdesc_mkl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){
206   uint_t j;
207     onset->data[0] = 0.;
208     for (j=0;j<fftgrain->length;j++) {
209       onset->data[0] += LOG(1.+fftgrain->norm[j]/(o->oldmag->data[j]+1.e-1));
210       o->oldmag->data[j] = fftgrain->norm[j];
211     }
212     if (isnan(onset->data[0])) onset->data[0] = 0.;
213 }
214
215 /* Spectral flux */
216 void aubio_specdesc_specflux(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){ 
217   uint_t j;
218   onset->data[0] = 0.;
219   for (j=0;j<fftgrain->length;j++) {
220     if (fftgrain->norm[j] > o->oldmag->data[j])
221       onset->data[0] += fftgrain->norm[j] - o->oldmag->data[j];
222     o->oldmag->data[j] = fftgrain->norm[j];
223   }
224 }
225
226 /* Generic function pointing to the choosen one */
227 void 
228 aubio_specdesc_do (aubio_specdesc_t *o, const cvec_t * fftgrain, 
229     fvec_t * onset) {
230   o->funcpointer(o,fftgrain,onset);
231 }
232
233 /* Allocate memory for an onset detection 
234  * depending on the choosen type, allocate memory as needed
235  */
236 aubio_specdesc_t * 
237 new_aubio_specdesc (const char_t * onset_mode, uint_t size){
238   aubio_specdesc_t * o = AUBIO_NEW(aubio_specdesc_t);
239   uint_t rsize = size/2+1;
240   aubio_specdesc_type onset_type;
241   if (strcmp (onset_mode, "energy") == 0)
242       onset_type = aubio_onset_energy;
243   else if (strcmp (onset_mode, "specdiff") == 0)
244       onset_type = aubio_onset_specdiff;
245   else if (strcmp (onset_mode, "hfc") == 0)
246       onset_type = aubio_onset_hfc;
247   else if (strcmp (onset_mode, "complexdomain") == 0)
248       onset_type = aubio_onset_complex;
249   else if (strcmp (onset_mode, "complex") == 0)
250       onset_type = aubio_onset_complex;
251   else if (strcmp (onset_mode, "phase") == 0)
252       onset_type = aubio_onset_phase;
253   else if (strcmp (onset_mode, "mkl") == 0)
254       onset_type = aubio_onset_mkl;
255   else if (strcmp (onset_mode, "kl") == 0)
256       onset_type = aubio_onset_kl;
257   else if (strcmp (onset_mode, "specflux") == 0)
258       onset_type = aubio_onset_specflux;
259   else if (strcmp (onset_mode, "centroid") == 0)
260       onset_type = aubio_specmethod_centroid;
261   else if (strcmp (onset_mode, "spread") == 0)
262       onset_type = aubio_specmethod_spread;
263   else if (strcmp (onset_mode, "skewness") == 0)
264       onset_type = aubio_specmethod_skewness;
265   else if (strcmp (onset_mode, "kurtosis") == 0)
266       onset_type = aubio_specmethod_kurtosis;
267   else if (strcmp (onset_mode, "slope") == 0)
268       onset_type = aubio_specmethod_slope;
269   else if (strcmp (onset_mode, "decrease") == 0)
270       onset_type = aubio_specmethod_decrease;
271   else if (strcmp (onset_mode, "rolloff") == 0)
272       onset_type = aubio_specmethod_rolloff;
273   else if (strcmp (onset_mode, "default") == 0)
274       onset_type = aubio_onset_default;
275   else {
276       AUBIO_ERR("unknown spectral descriptor type %s\n", onset_mode);
277       AUBIO_FREE(o);
278       return NULL;
279   }
280   switch(onset_type) {
281     /* for both energy and hfc, only fftgrain->norm is required */
282     case aubio_onset_energy:
283       break;
284     case aubio_onset_hfc:
285       break;
286       /* the other approaches will need some more memory spaces */
287     case aubio_onset_complex:
288       o->oldmag = new_fvec(rsize);
289       o->dev1   = new_fvec(rsize);
290       o->theta1 = new_fvec(rsize);
291       o->theta2 = new_fvec(rsize);
292       break;
293     case aubio_onset_phase:
294       o->dev1   = new_fvec(rsize);
295       o->theta1 = new_fvec(rsize);
296       o->theta2 = new_fvec(rsize);
297       o->histog = new_aubio_hist(0.0, PI, 10);
298       o->threshold = 0.1;
299       break;
300     case aubio_onset_specdiff:
301       o->oldmag = new_fvec(rsize);
302       o->dev1   = new_fvec(rsize);
303       o->histog = new_aubio_hist(0.0, PI, 10);
304       o->threshold = 0.1;
305       break;
306     case aubio_onset_kl:
307     case aubio_onset_mkl:
308     case aubio_onset_specflux:
309       o->oldmag = new_fvec(rsize);
310       break;
311     default:
312       break;
313   }
314
315   switch(onset_type) {
316     case aubio_onset_energy:
317       o->funcpointer = aubio_specdesc_energy;
318       break;
319     case aubio_onset_hfc:
320       o->funcpointer = aubio_specdesc_hfc;
321       break;
322     case aubio_onset_complex:
323       o->funcpointer = aubio_specdesc_complex;
324       break;
325     case aubio_onset_phase:
326       o->funcpointer = aubio_specdesc_phase;
327       break;
328     case aubio_onset_specdiff:
329       o->funcpointer = aubio_specdesc_specdiff;
330       break;
331     case aubio_onset_kl:
332       o->funcpointer = aubio_specdesc_kl;
333       break;
334     case aubio_onset_mkl:
335       o->funcpointer = aubio_specdesc_mkl;
336       break;
337     case aubio_onset_specflux:
338       o->funcpointer = aubio_specdesc_specflux;
339       break;
340     case aubio_specmethod_centroid:
341       o->funcpointer = aubio_specdesc_centroid;
342       break;
343     case aubio_specmethod_spread:
344       o->funcpointer = aubio_specdesc_spread;
345       break;
346     case aubio_specmethod_skewness:
347       o->funcpointer = aubio_specdesc_skewness;
348       break;
349     case aubio_specmethod_kurtosis:
350       o->funcpointer = aubio_specdesc_kurtosis;
351       break;
352     case aubio_specmethod_slope:
353       o->funcpointer = aubio_specdesc_slope;
354       break;
355     case aubio_specmethod_decrease:
356       o->funcpointer = aubio_specdesc_decrease;
357       break;
358     case aubio_specmethod_rolloff:
359       o->funcpointer = aubio_specdesc_rolloff;
360       break;
361     default:
362       break;
363   }
364   o->onset_type = onset_type;
365   return o;
366 }
367
368 void del_aubio_specdesc (aubio_specdesc_t *o){
369   switch(o->onset_type) {
370     case aubio_onset_energy:
371       break;
372     case aubio_onset_hfc:
373       break;
374     case aubio_onset_complex:
375       del_fvec(o->oldmag);
376       del_fvec(o->dev1);
377       del_fvec(o->theta1);
378       del_fvec(o->theta2);
379       break;
380     case aubio_onset_phase:
381       del_fvec(o->dev1);
382       del_fvec(o->theta1);
383       del_fvec(o->theta2);
384       del_aubio_hist(o->histog);
385       break;
386     case aubio_onset_specdiff:
387       del_fvec(o->oldmag);
388       del_fvec(o->dev1);
389       del_aubio_hist(o->histog);
390       break;
391     case aubio_onset_kl:
392     case aubio_onset_mkl:
393     case aubio_onset_specflux:
394       del_fvec(o->oldmag);
395       break;
396     default:
397       break;
398   }
399   AUBIO_FREE(o);
400 }