src/pitchyin.c
[aubio.git] / src / pitchyin.c
1 /*
2    Copyright (C) 2003 Paul Brossier
3
4    This program is free software; you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation; either version 2 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program; if not, write to the Free Software
16    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
17 */
18
19 /* This algorithm was developped by A. de Cheveigne and H. Kawahara and
20  * published in:
21  * 
22  * de Cheveign√©, A., Kawahara, H. (2002) "YIN, a fundamental frequency
23  * estimator for speech and music", J. Acoust. Soc. Am. 111, 1917-1930.  
24  *
25  * see http://recherche.ircam.fr/equipes/pcm/pub/people/cheveign.html
26  */
27
28 #include "aubio_priv.h"
29 #include "sample.h"
30 #include "mathutils.h"
31 #include "pitchyin.h"
32
33 /* outputs the difference function */
34 void aubio_pitchyin_diff(fvec_t * input, fvec_t * yin){
35         uint_t c,j,tau;
36         smpl_t tmp;
37         for (c=0;c<input->channels;c++)
38         {
39                 for (tau=0;tau<yin->length;tau++)
40                 {
41                         yin->data[c][tau] = 0.;
42                 }
43                 for (tau=1;tau<yin->length;tau++)
44                 {
45                         for (j=0;j<yin->length;j++)
46                         {
47                                 tmp = input->data[c][j] - input->data[c][j+tau];
48                                 yin->data[c][tau] += SQR(tmp);
49                         }
50                 }
51         }
52 }
53
54 /* cumulative mean normalized difference function */
55 void aubio_pitchyin_getcum(fvec_t * yin) {
56         uint_t c,tau;
57         smpl_t tmp;
58         for (c=0;c<yin->channels;c++)
59         {
60                 tmp = 0.;
61                 yin->data[c][0] = 1.;
62                 //AUBIO_DBG("%f\t",yin->data[c][0]);
63                 for (tau=1;tau<yin->length;tau++)
64                 {
65                         tmp += yin->data[c][tau];
66                         yin->data[c][tau] *= tau/tmp;
67                         //AUBIO_DBG("%f\t",yin->data[c][tau]);
68                 }
69                 //AUBIO_DBG("\n");
70         }
71 }
72
73 uint_t aubio_pitchyin_getpitch(fvec_t * yin) {
74         uint_t c=0,tau=1;
75         do 
76         {
77                 if(yin->data[c][tau] < 0.1) { 
78                         while (yin->data[c][tau+1] < yin->data[c][tau]) {
79                                 tau++;
80                         }
81                         return tau;
82                 }
83                 tau++;
84         } while (tau<yin->length);
85         //AUBIO_DBG("No pitch found");
86         return 0;
87 }
88
89
90 /* all the above in one */
91 void aubio_pitchyin_getpitchfast(fvec_t * input, fvec_t * yin){
92         uint_t c,j,tau;
93         smpl_t tmp, tmp2;
94         for (c=0;c<input->channels;c++)
95         {
96                 for (tau=0;tau<yin->length;tau++)
97                 {
98                         yin->data[c][tau] = 0.;
99                 }
100                 for (tau=1;tau<yin->length;tau++)
101                 {
102                         for (j=0;j<yin->length;j++)
103                         {
104                                 tmp = input->data[c][j] - input->data[c][j+tau];
105                                 yin->data[c][tau] += SQR(tmp);
106                         }
107                 }
108                 tmp2 = 0.;
109                 yin->data[c][0] = 1.;
110                 for (tau=1;tau<yin->length;tau++)
111                 {
112                         tmp += yin->data[c][tau];
113                         yin->data[c][tau] *= tau/tmp;
114                 }
115         }
116         /* should merge the following with above */
117         do 
118         {
119                 if(yin->data[c][tau] < 0.1) { 
120                         while (yin->data[c][tau+1] < yin->data[c][tau]) {
121                                 tau++;
122                         }
123                         return tau;
124                 }
125                 tau++;
126         } while (tau<yin->length);
127         return 0;
128 }
129