Changeset fddfa64 for src/pitch/pitchyinfft.c
- Timestamp:
- Nov 3, 2009, 4:14:03 PM (14 years ago)
- Branches:
- feature/autosink, feature/cnn, feature/cnn_org, feature/constantq, feature/crepe, feature/crepe_org, feature/pitchshift, feature/pydocstrings, feature/timestretch, fix/ffmpeg5, master, pitchshift, sampler, timestretch, yinfft+
- Children:
- bafe71d
- Parents:
- 63f3c70
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/pitch/pitchyinfft.c
r63f3c70 rfddfa64 27 27 28 28 /** pitch yinfft structure */ 29 struct _aubio_pitchyinfft_t { 30 fvec_t * win; /**< temporal weighting window */ 31 fvec_t * winput; /**< windowed spectrum */ 32 cvec_t * res; /**< complex vector to compute square difference function */ 33 fvec_t * sqrmag; /**< square difference function */ 34 fvec_t * weight; /**< spectral weighting window (psychoacoustic model) */ 35 cvec_t * fftout; /**< Fourier transform output */ 36 aubio_fft_t * fft; /**< fft object to compute square difference function */ 37 fvec_t * yinfft; /**< Yin function */ 29 struct _aubio_pitchyinfft_t 30 { 31 fvec_t *win; /**< temporal weighting window */ 32 fvec_t *winput; /**< windowed spectrum */ 33 cvec_t *res; /**< complex vector to compute square difference function */ 34 fvec_t *sqrmag; /**< square difference function */ 35 fvec_t *weight; /**< spectral weighting window (psychoacoustic model) */ 36 cvec_t *fftout; /**< Fourier transform output */ 37 aubio_fft_t *fft; /**< fft object to compute square difference function */ 38 fvec_t *yinfft; /**< Yin function */ 38 39 smpl_t tol; /**< Yin tolerance */ 39 40 }; 40 41 41 static const smpl_t freqs[] = { 0., 20., 25., 31.5, 40., 50., 63., 80., 100.,42 static const smpl_t freqs[] = { 0., 20., 25., 31.5, 40., 50., 63., 80., 100., 42 43 125., 160., 200., 250., 315., 400., 500., 630., 800., 1000., 1250., 43 44 1600., 2000., 2500., 3150., 4000., 5000., 6300., 8000., 9000., 10000., 44 12500., 15000., 20000., 25100}; 45 12500., 15000., 20000., 25100 46 }; 45 47 46 static const smpl_t weight[] = { -75.8, -70.1, -60.8, -52.1, -44.2, -37.5,48 static const smpl_t weight[] = { -75.8, -70.1, -60.8, -52.1, -44.2, -37.5, 47 49 -31.3, -25.6, -20.9, -16.5, -12.6, -9.6, -7.0, -4.7, -3.0, -1.8, -0.8, 48 50 -0.2, -0.0, 0.5, 1.6, 3.2, 5.4, 7.8, 8.1, 5.3, -2.4, -11.1, -12.8, 49 -12.2, -7.4, -17.8, -17.8, -17.8}; 51 -12.2, -7.4, -17.8, -17.8, -17.8 52 }; 50 53 51 aubio_pitchyinfft_t * new_aubio_pitchyinfft (uint_t bufsize) 54 aubio_pitchyinfft_t * 55 new_aubio_pitchyinfft (uint_t bufsize) 52 56 { 53 aubio_pitchyinfft_t * p = AUBIO_NEW(aubio_pitchyinfft_t);54 p->winput = new_fvec (bufsize,1);55 p->fft = new_aubio_fft(bufsize, 1);56 p->fftout = new_cvec (bufsize,1);57 p->sqrmag = new_fvec (bufsize,1);58 p->res = new_cvec(bufsize,1);59 p->yinfft = new_fvec (bufsize/2+1,1);60 p->tol 61 p->win = new_aubio_window("hanningz", bufsize);62 p->weight = new_fvec(bufsize/2+1,1);57 aubio_pitchyinfft_t *p = AUBIO_NEW (aubio_pitchyinfft_t); 58 p->winput = new_fvec (bufsize, 1); 59 p->fft = new_aubio_fft (bufsize, 1); 60 p->fftout = new_cvec (bufsize, 1); 61 p->sqrmag = new_fvec (bufsize, 1); 62 p->res = new_cvec (bufsize, 1); 63 p->yinfft = new_fvec (bufsize / 2 + 1, 1); 64 p->tol = 0.85; 65 p->win = new_aubio_window ("hanningz", bufsize); 66 p->weight = new_fvec (bufsize / 2 + 1, 1); 63 67 { 64 68 uint_t i = 0, j = 1; 65 69 smpl_t freq = 0, a0 = 0, a1 = 0, f0 = 0, f1 = 0; 66 for (i =0; i<p->weight->length; i++) {67 freq = (smpl_t) i/(smpl_t)bufsize*(smpl_t)44100.;70 for (i = 0; i < p->weight->length; i++) { 71 freq = (smpl_t) i / (smpl_t) bufsize *(smpl_t) 44100.; 68 72 while (freq > freqs[j]) { 69 j += 1;70 71 a0 = weight[j -1];72 f0 = freqs[j -1];73 73 j += 1; 74 } 75 a0 = weight[j - 1]; 76 f0 = freqs[j - 1]; 77 a1 = weight[j]; 74 78 f1 = freqs[j]; 75 if (f0 == f1) { // just in case79 if (f0 == f1) { // just in case 76 80 p->weight->data[0][i] = a0; 77 } else if (f0 == 0) { // y = ax+b78 p->weight->data[0][i] = (a1 -a0)/f1*freq + a0;81 } else if (f0 == 0) { // y = ax+b 82 p->weight->data[0][i] = (a1 - a0) / f1 * freq + a0; 79 83 } else { 80 p->weight->data[0][i] = (a1 -a0)/(f1-f0)*freq +81 (a0 - (a1 - a0)/(f1/f0 - 1.));84 p->weight->data[0][i] = (a1 - a0) / (f1 - f0) * freq + 85 (a0 - (a1 - a0) / (f1 / f0 - 1.)); 82 86 } 83 87 while (freq > freqs[j]) { 84 j += 1;88 j += 1; 85 89 } 86 90 //AUBIO_DBG("%f\n",p->weight->data[0][i]); 87 p->weight->data[0][i] = DB2LIN (p->weight->data[0][i]);91 p->weight->data[0][i] = DB2LIN (p->weight->data[0][i]); 88 92 //p->weight->data[0][i] = SQRT(DB2LIN(p->weight->data[0][i])); 89 93 } … … 92 96 } 93 97 94 void aubio_pitchyinfft_do (aubio_pitchyinfft_t * p, fvec_t * input, fvec_t * output) { 98 void 99 aubio_pitchyinfft_do (aubio_pitchyinfft_t * p, fvec_t * input, fvec_t * output) 100 { 95 101 uint_t i, tau, l; 96 102 uint_t halfperiod; 97 103 smpl_t tmp, sum; 98 cvec_t * res = (cvec_t *)p->res; 99 fvec_t * yin = (fvec_t *)p->yinfft; 100 for (i=0; i < input->channels; i++){ 101 l = 0; tmp = 0.; sum = 0.; 102 for (l=0; l < input->length; l++){ 103 p->winput->data[0][l] = p->win->data[0][l] * input->data[i][l]; 104 } 105 aubio_fft_do(p->fft,p->winput,p->fftout); 106 for (l=0; l < p->fftout->length; l++){ 107 p->sqrmag->data[0][l] = SQR(p->fftout->norm[0][l]); 108 p->sqrmag->data[0][l] *= p->weight->data[0][l]; 109 } 110 for (l=1; l < p->fftout->length; l++){ 111 p->sqrmag->data[0][(p->fftout->length-1)*2-l] = 112 SQR(p->fftout->norm[0][l]); 113 p->sqrmag->data[0][(p->fftout->length-1)*2-l] *= 114 p->weight->data[0][l]; 115 } 116 for (l=0; l < p->sqrmag->length/2+1; l++) { 117 sum += p->sqrmag->data[0][l]; 118 } 119 sum *= 2.; 120 aubio_fft_do(p->fft,p->sqrmag,res); 121 yin->data[0][0] = 1.; 122 for (tau=1; tau < yin->length; tau++) { 123 yin->data[0][tau] = sum - 124 res->norm[0][tau]*COS(res->phas[0][tau]); 125 tmp += yin->data[0][tau]; 126 yin->data[0][tau] *= tau/tmp; 127 } 128 tau = fvec_min_elem(yin); 129 if (yin->data[0][tau] < p->tol) { 130 /* no interpolation */ 131 //return tau; 132 /* 3 point quadratic interpolation */ 133 //return fvec_quadint_min(yin,tau,1); 134 /* additional check for (unlikely) octave doubling in higher frequencies */ 135 if (tau>35) { 136 output->data[i][0] = fvec_quadint(yin,tau,i); 104 cvec_t *res = (cvec_t *) p->res; 105 fvec_t *yin = (fvec_t *) p->yinfft; 106 for (i = 0; i < input->channels; i++) { 107 l = 0; 108 tmp = 0.; 109 sum = 0.; 110 for (l = 0; l < input->length; l++) { 111 p->winput->data[0][l] = p->win->data[0][l] * input->data[i][l]; 112 } 113 aubio_fft_do (p->fft, p->winput, p->fftout); 114 for (l = 0; l < p->fftout->length; l++) { 115 p->sqrmag->data[0][l] = SQR (p->fftout->norm[0][l]); 116 p->sqrmag->data[0][l] *= p->weight->data[0][l]; 117 } 118 for (l = 1; l < p->fftout->length; l++) { 119 p->sqrmag->data[0][(p->fftout->length - 1) * 2 - l] = 120 SQR (p->fftout->norm[0][l]); 121 p->sqrmag->data[0][(p->fftout->length - 1) * 2 - l] *= 122 p->weight->data[0][l]; 123 } 124 for (l = 0; l < p->sqrmag->length / 2 + 1; l++) { 125 sum += p->sqrmag->data[0][l]; 126 } 127 sum *= 2.; 128 aubio_fft_do (p->fft, p->sqrmag, res); 129 yin->data[0][0] = 1.; 130 for (tau = 1; tau < yin->length; tau++) { 131 yin->data[0][tau] = sum - res->norm[0][tau] * COS (res->phas[0][tau]); 132 tmp += yin->data[0][tau]; 133 yin->data[0][tau] *= tau / tmp; 134 } 135 tau = fvec_min_elem (yin); 136 if (yin->data[0][tau] < p->tol) { 137 /* no interpolation */ 138 //return tau; 139 /* 3 point quadratic interpolation */ 140 //return fvec_quadint_min(yin,tau,1); 141 /* additional check for (unlikely) octave doubling in higher frequencies */ 142 if (tau > 35) { 143 output->data[i][0] = fvec_quadint (yin, tau, i); 144 } else { 145 /* should compare the minimum value of each interpolated peaks */ 146 halfperiod = FLOOR (tau / 2 + .5); 147 if (yin->data[0][halfperiod] < p->tol) 148 output->data[i][0] = fvec_quadint (yin, halfperiod, i); 149 else 150 output->data[i][0] = fvec_quadint (yin, tau, i); 151 } 137 152 } else { 138 /* should compare the minimum value of each interpolated peaks */ 139 halfperiod = FLOOR(tau/2+.5); 140 if (yin->data[0][halfperiod] < p->tol) 141 output->data[i][0] = fvec_quadint(yin,halfperiod,i); 142 else 143 output->data[i][0] = fvec_quadint(yin,tau,i); 153 output->data[i][0] = 0.; 144 154 } 145 } else {146 output->data[i][0] = 0.;147 }148 155 } 149 156 } 150 157 151 void del_aubio_pitchyinfft(aubio_pitchyinfft_t *p){ 152 del_fvec(p->win); 153 del_aubio_fft(p->fft); 154 del_fvec(p->yinfft); 155 del_fvec(p->sqrmag); 156 del_cvec(p->res); 157 del_cvec(p->fftout); 158 del_fvec(p->winput); 159 del_fvec(p->weight); 160 AUBIO_FREE(p); 158 void 159 del_aubio_pitchyinfft (aubio_pitchyinfft_t * p) 160 { 161 del_fvec (p->win); 162 del_aubio_fft (p->fft); 163 del_fvec (p->yinfft); 164 del_fvec (p->sqrmag); 165 del_cvec (p->res); 166 del_cvec (p->fftout); 167 del_fvec (p->winput); 168 del_fvec (p->weight); 169 AUBIO_FREE (p); 161 170 } 162 171 163 uint_t aubio_pitchyinfft_set_tolerance (aubio_pitchyinfft_t * p, smpl_t tol) { 172 uint_t 173 aubio_pitchyinfft_set_tolerance (aubio_pitchyinfft_t * p, smpl_t tol) 174 { 164 175 p->tol = tol; 165 176 return 0; 166 177 } 167 178 168 smpl_t aubio_pitchyinfft_get_tolerance (aubio_pitchyinfft_t * p) { 179 smpl_t 180 aubio_pitchyinfft_get_tolerance (aubio_pitchyinfft_t * p) 181 { 169 182 return p->tol; 170 183 }
Note: See TracChangeset
for help on using the changeset viewer.