Stretch.cpp (10349B)
1 /* 2 Copyright (C) 2006-2011 Nasca Octavian Paul 3 Author: Nasca Octavian Paul 4 5 This program is free software; you can redistribute it and/or modify 6 it under the terms of version 2 of the GNU General Public License 7 as published by the Free Software Foundation. 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 (version 2) for more details. 13 14 You should have received a copy of the GNU General Public License (version 2) 15 along with this program; if not, write to the Free Software Foundation, 16 Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 17 */ 18 19 #include "Stretch.h" 20 #include <stdlib.h> 21 #include <math.h> 22 23 unsigned int FFT::start_rand_seed=1; 24 25 FFT::FFT(int nsamples_){ 26 nsamples=nsamples_; 27 if (nsamples%2!=0) { 28 nsamples+=1; 29 printf("WARNING: Odd sample size on FFT::FFT() (%d)",nsamples); 30 }; 31 smp=new REALTYPE[nsamples];for (int i=0;i<nsamples;i++) smp[i]=0.0; 32 freq=new REALTYPE[nsamples/2+1];for (int i=0;i<nsamples/2+1;i++) freq[i]=0.0; 33 window.data=new REALTYPE[nsamples];for (int i=0;i<nsamples;i++) window.data[i]=0.707; 34 window.type=W_RECTANGULAR; 35 36 #ifdef KISSFFT 37 datar=new kiss_fft_scalar[nsamples+2]; 38 for (int i=0;i<nsamples+2;i++) datar[i]=0.0; 39 datac=new kiss_fft_cpx[nsamples/2+2]; 40 for (int i=0;i<nsamples/2+2;i++) datac[i].r=datac[i].i=0.0; 41 plankfft = kiss_fftr_alloc(nsamples,0,0,0); 42 plankifft = kiss_fftr_alloc(nsamples,1,0,0); 43 #else 44 data=new REALTYPE[nsamples];for (int i=0;i<nsamples;i++) data[i]=0.0; 45 planfftw=fftwf_plan_r2r_1d(nsamples,data,data,FFTW_R2HC,FFTW_ESTIMATE); 46 planifftw=fftwf_plan_r2r_1d(nsamples,data,data,FFTW_HC2R,FFTW_ESTIMATE); 47 #endif 48 rand_seed=start_rand_seed; 49 start_rand_seed+=161103; 50 }; 51 52 FFT::~FFT(){ 53 delete []smp; 54 delete []freq; 55 delete []window.data; 56 #ifdef KISSFFT 57 delete []datar; 58 delete []datac; 59 free(plankfft); 60 free(plankifft); 61 #else 62 delete []data; 63 fftwf_destroy_plan(planfftw); 64 fftwf_destroy_plan(planifftw); 65 #endif 66 }; 67 68 void FFT::smp2freq(){ 69 #ifdef KISSFFT 70 for (int i=0;i<nsamples;i++) datar[i]=smp[i]; 71 kiss_fftr(plankfft,datar,datac); 72 #else 73 for (int i=0;i<nsamples;i++) data[i]=smp[i]; 74 fftwf_execute(planfftw); 75 #endif 76 77 for (int i=1;i<nsamples/2;i++) { 78 #ifdef KISSFFT 79 REALTYPE c=datac[i].r; 80 REALTYPE s=datac[i].i; 81 #else 82 REALTYPE c=data[i]; 83 REALTYPE s=data[nsamples-i]; 84 #endif 85 freq[i]=sqrt(c*c+s*s); 86 }; 87 freq[0]=0.0; 88 }; 89 90 void FFT::freq2smp(){ 91 REALTYPE inv_2p15_2pi=1.0/16384.0*M_PI; 92 for (int i=1;i<nsamples/2;i++) { 93 rand_seed=(rand_seed*1103515245+12345); 94 unsigned int rand=(rand_seed>>16)&0x7fff; 95 REALTYPE phase=rand*inv_2p15_2pi; 96 #ifdef KISSFFT 97 datac[i].r=freq[i]*cos(phase); 98 datac[i].i=freq[i]*sin(phase); 99 #else 100 data[i]=freq[i]*cos(phase); 101 data[nsamples-i]=freq[i]*sin(phase); 102 #endif 103 }; 104 105 #ifdef KISSFFT 106 datac[0].r=datac[0].i=0.0; 107 kiss_fftri(plankifft,datac,datar); 108 for (int i=0;i<nsamples;i++) smp[i]=datar[i]/nsamples; 109 #else 110 data[0]=data[nsamples/2+1]=data[nsamples/2]=0.0; 111 fftwf_execute(planifftw); 112 for (int i=0;i<nsamples;i++) smp[i]=data[i]/nsamples; 113 #endif 114 }; 115 116 void FFT::applywindow(FFTWindow type){ 117 if (window.type!=type){ 118 window.type=type; 119 switch (type){ 120 case W_RECTANGULAR: 121 for (int i=0;i<nsamples;i++) window.data[i]=0.707; 122 break; 123 case W_HAMMING: 124 for (int i=0;i<nsamples;i++) window.data[i]=0.53836-0.46164*cos(2*M_PI*i/(nsamples+1.0)); 125 break; 126 case W_HANN: 127 for (int i=0;i<nsamples;i++) window.data[i]=0.5*(1.0-cos(2*M_PI*i/(nsamples-1.0))); 128 break; 129 case W_BLACKMAN: 130 for (int i=0;i<nsamples;i++) window.data[i]=0.42-0.5*cos(2*M_PI*i/(nsamples-1.0))+0.08*cos(4*M_PI*i/(nsamples-1.0)); 131 break; 132 case W_BLACKMAN_HARRIS: 133 for (int i=0;i<nsamples;i++) window.data[i]=0.35875-0.48829*cos(2*M_PI*i/(nsamples-1.0))+0.14128*cos(4*M_PI*i/(nsamples-1.0))-0.01168*cos(6*M_PI*i/(nsamples-1.0)); 134 break; 135 136 }; 137 }; 138 for (int i=0;i<nsamples;i++) smp[i]*=window.data[i]; 139 }; 140 141 /*******************************************/ 142 143 144 Stretch::Stretch(REALTYPE rap_,int bufsize_,FFTWindow w,bool bypass_,REALTYPE samplerate_,int stereo_mode_){ 145 freezing=false; 146 onset_detection_sensitivity=0.0; 147 148 samplerate=samplerate_; 149 rap=rap_; 150 bufsize=bufsize_; 151 bypass=bypass_; 152 stereo_mode=stereo_mode_; 153 if (bufsize<8) bufsize=8; 154 155 out_buf=new REALTYPE[bufsize]; 156 old_freq=new REALTYPE[bufsize]; 157 158 very_old_smps=new REALTYPE[bufsize]; 159 new_smps=new REALTYPE[bufsize]; 160 old_smps=new REALTYPE[bufsize]; 161 162 old_out_smps=new REALTYPE[bufsize*2]; 163 for (int i=0;i<bufsize*2;i++) { 164 old_out_smps[i]=0.0; 165 }; 166 for (int i=0;i<bufsize;i++) { 167 old_freq[i]=0.0; 168 new_smps[i]=0.0; 169 old_smps[i]=0.0; 170 very_old_smps[i]=0.0; 171 }; 172 173 174 infft=new FFT(bufsize*2); 175 fft=new FFT(bufsize*2); 176 outfft=new FFT(bufsize*2); 177 remained_samples=0.0; 178 window_type=w; 179 require_new_buffer=false; 180 c_pos_percents=0.0; 181 extra_onset_time_credit=0.0; 182 skip_samples=0; 183 }; 184 185 Stretch::~Stretch(){ 186 delete [] old_freq; 187 delete [] out_buf; 188 delete [] new_smps; 189 delete [] old_smps; 190 delete [] very_old_smps; 191 delete [] old_out_smps; 192 delete fft; 193 delete infft; 194 delete outfft; 195 }; 196 197 void Stretch::set_rap(REALTYPE newrap){ 198 //if ((rap>=1.0)&&(newrap>=1.0)) 199 rap=newrap; 200 }; 201 202 void Stretch::do_analyse_inbuf(REALTYPE *smps){ 203 //get the frequencies 204 for (int i=0;i<bufsize;i++) { 205 infft->smp[i]=old_smps[i]; 206 infft->smp[i+bufsize]=smps[i]; 207 208 old_freq[i]=infft->freq[i]; 209 }; 210 infft->applywindow(window_type); 211 infft->smp2freq(); 212 }; 213 214 void Stretch::do_next_inbuf_smps(REALTYPE *smps){ 215 for (int i=0;i<bufsize;i++) { 216 very_old_smps[i]=old_smps[i]; 217 old_smps[i]=new_smps[i]; 218 new_smps[i]=smps[i]; 219 }; 220 }; 221 222 REALTYPE Stretch::do_detect_onset(){ 223 REALTYPE result=0.0; 224 if (onset_detection_sensitivity>1e-3){ 225 REALTYPE os=0.0,osinc=0.0; 226 REALTYPE osincold=1e-5; 227 int maxk=1+(int)(bufsize*500.0/(samplerate*0.5)); 228 int k=0; 229 for (int i=0;i<bufsize;i++) { 230 osinc+=infft->freq[i]-old_freq[i]; 231 osincold+=old_freq[i]; 232 if (k>=maxk) { 233 k=0; 234 os+=osinc/osincold; 235 osinc=0; 236 }; 237 k++; 238 }; 239 os+=osinc; 240 if (os<0.0) os=0.0; 241 //if (os>1.0) os=1.0; 242 243 REALTYPE os_strength=pow(20.0,1.0-onset_detection_sensitivity)-1.0; 244 REALTYPE os_strength_h=os_strength*0.75; 245 if (os>os_strength_h){ 246 result=(os-os_strength_h)/(os_strength-os_strength_h); 247 if (result>1.0) result=1.0; 248 }; 249 250 if (result>1.0) result=1.0; 251 }; 252 return result; 253 }; 254 255 REALTYPE Stretch::process(REALTYPE *smps,int nsmps){ 256 REALTYPE onset=0.0; 257 if (bypass){ 258 for (int i=0;i<bufsize;i++) out_buf[i]=smps[i]; 259 return 0.0; 260 }; 261 262 if (smps!=NULL){ 263 if ((nsmps!=0)&&(nsmps!=bufsize)&&(nsmps!=get_max_bufsize())){ 264 printf("Warning wrong nsmps on Stretch::process() %d,%d\n",nsmps,bufsize); 265 return 0.0; 266 }; 267 if (nsmps!=0){//new data arrived: update the frequency components 268 do_analyse_inbuf(smps); 269 if (nsmps==get_max_bufsize()) { 270 for (int k=bufsize;k<get_max_bufsize();k+=bufsize) do_analyse_inbuf(smps+k); 271 }; 272 if (onset_detection_sensitivity>1e-3) onset=do_detect_onset(); 273 }; 274 275 276 //move the buffers 277 if (nsmps!=0){//new data arrived: update the frequency components 278 do_next_inbuf_smps(smps); 279 if (nsmps==get_max_bufsize()) { 280 for (int k=bufsize;k<get_max_bufsize();k+=bufsize) do_next_inbuf_smps(smps+k); 281 282 }; 283 }; 284 285 //construct the input fft 286 int start_pos=(int)(floor(remained_samples*bufsize)); 287 if (start_pos>=bufsize) start_pos=bufsize-1; 288 for (int i=0;i<bufsize-start_pos;i++) fft->smp[i]=very_old_smps[i+start_pos]; 289 for (int i=0;i<bufsize;i++) fft->smp[i+bufsize-start_pos]=old_smps[i]; 290 for (int i=0;i<start_pos;i++) fft->smp[i+2*bufsize-start_pos]=new_smps[i]; 291 //compute the output spectrum 292 fft->applywindow(window_type); 293 fft->smp2freq(); 294 for (int i=0;i<bufsize;i++) outfft->freq[i]=fft->freq[i]; 295 296 297 298 //for (int i=0;i<bufsize;i++) outfft->freq[i]=infft->freq[i]*remained_samples+old_freq[i]*(1.0-remained_samples); 299 300 301 process_spectrum(outfft->freq); 302 303 outfft->freq2smp(); 304 305 //make the output buffer 306 REALTYPE tmp=1.0/(float) bufsize*M_PI; 307 REALTYPE hinv_sqrt2=0.853553390593;//(1.0+1.0/sqrt(2))*0.5; 308 309 REALTYPE ampfactor=2.0; 310 311 //remove the resulted unwanted amplitude modulation (caused by the interference of N and N+1 windowed buffer and compute the output buffer 312 for (int i=0;i<bufsize;i++) { 313 REALTYPE a=(0.5+0.5*cos(i*tmp)); 314 REALTYPE out=outfft->smp[i+bufsize]*(1.0-a)+old_out_smps[i]*a; 315 out_buf[i]=out*(hinv_sqrt2-(1.0-hinv_sqrt2)*cos(i*2.0*tmp))*ampfactor; 316 }; 317 318 //copy the current output buffer to old buffer 319 for (int i=0;i<bufsize*2;i++) old_out_smps[i]=outfft->smp[i]; 320 321 }; 322 323 if (!freezing){ 324 long double used_rap=rap*get_stretch_multiplier(c_pos_percents); 325 326 long double r=1.0/used_rap; 327 if (extra_onset_time_credit>0){ 328 REALTYPE credit_get=0.5*r;//must be smaller than r 329 extra_onset_time_credit-=credit_get; 330 if (extra_onset_time_credit<0.0) extra_onset_time_credit=0.0; 331 r-=credit_get; 332 }; 333 334 long double old_remained_samples_test=remained_samples; 335 remained_samples+=r; 336 int result=0; 337 if (remained_samples>=1.0){ 338 skip_samples=(int)(floor(remained_samples-1.0)*bufsize); 339 remained_samples=remained_samples-floor(remained_samples); 340 require_new_buffer=true; 341 }else{ 342 require_new_buffer=false; 343 }; 344 }; 345 // long double rf_test=remained_samples-old_remained_samples_test;//this value should be almost like "rf" (for most of the time with the exception of changing the "ri" value) for extremely long stretches (otherwise the shown stretch value is not accurate) 346 //for stretch up to 10^18x "long double" must have at least 64 bits in the fraction part (true for gcc compiler on x86 and macosx) 347 return onset; 348 }; 349 350 void Stretch::here_is_onset(REALTYPE onset){ 351 if (freezing) return; 352 if (onset>0.5){ 353 require_new_buffer=true; 354 extra_onset_time_credit+=1.0-remained_samples; 355 remained_samples=0.0; 356 skip_samples=0; 357 }; 358 }; 359 360 int Stretch::get_nsamples(REALTYPE current_pos_percents){ 361 if (bypass) return bufsize; 362 if (freezing) return 0; 363 c_pos_percents=current_pos_percents; 364 return require_new_buffer?bufsize:0; 365 }; 366 367 int Stretch::get_nsamples_for_fill(){ 368 return get_max_bufsize(); 369 }; 370 371 int Stretch::get_skip_nsamples(){ 372 if (freezing||bypass) return 0; 373 return skip_samples; 374 }; 375 376 REALTYPE Stretch::get_stretch_multiplier(REALTYPE pos_percents){ 377 return 1.0; 378 }; 379