CLAM-Development
1.1
|
00001 /* 00002 * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG) 00003 * UNIVERSITAT POMPEU FABRA 00004 * 00005 * 00006 * This program is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program; if not, write to the Free Software 00018 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00019 * 00020 */ 00021 00022 #include "ProcessingData.hxx" 00023 #include "DataTypes.hxx" 00024 #include "Spectrum.hxx" 00025 #include "Audio.hxx" 00026 #include "WindowGenerator.hxx" 00027 #include "ProcessingFactory.hxx" 00028 00029 namespace CLAM 00030 { 00031 00032 namespace Hidden 00033 { 00034 static const char * metadata[] = { 00035 "key", "WindowGenerator", 00036 "category", "Generators", // must change to Analysis? 00037 "description", "WindowGenerator", 00038 0 00039 }; 00040 static FactoryRegistrator<ProcessingFactory, WindowGenerator> reg = metadata; 00041 } 00042 00043 /* Processing object Method implementations */ 00044 00045 WindowGenerator::WindowGenerator(const WindowGeneratorConfig &c) 00046 : mOutput( "Generated window function samples", this ) 00047 , mSize("Size",this) 00048 { 00049 Configure(c); 00050 } 00051 00052 WindowGenerator::~WindowGenerator() 00053 {} 00054 00055 00056 /* Configure the Processing Object according to the Config object */ 00057 00058 bool WindowGenerator::ConcreteConfigure(const ProcessingConfig& c) 00059 { 00060 CopyAsConcreteConfig(mConfig, c); 00061 mSize.DoControl(TControlData(mConfig.GetSize())); 00062 00063 if (!mConfig.HasUseTable()) return true; 00064 if (!mConfig.GetUseTable()) return true; 00065 00066 00067 /* Fill the table */ 00068 00069 mTable.Resize(mConfig.GetSize()); 00070 mTable.SetSize(mConfig.GetSize()); 00071 mSize.DoControl(TControlData(mConfig.GetSize())); 00072 00073 EWindowType type = mConfig.HasType()? 00074 mConfig.GetType() : 00075 EWindowType::eHamming; 00076 00077 CreateTable(mTable,type,mConfig.GetSize()); 00078 00079 return true; 00080 00081 } 00082 00083 /* The supervised Do() function */ 00084 00085 bool WindowGenerator::Do() 00086 { 00087 CLAM_ASSERT( AbleToExecute(), "This processing is not ready to do anything" ); 00088 00089 bool result = Do( mOutput.GetAudio() ); 00090 mOutput.Produce(); 00091 return result; 00092 } 00093 00094 /* The unsupervised Do() function */ 00095 00096 bool WindowGenerator::Do(DataArray& out) 00097 { 00098 const int winsize = (int) mSize.GetLastValue(); 00099 const int audiosize = out.Size(); 00100 00101 bool useTable = mConfig.HasUseTable() && mConfig.GetUseTable(); 00102 00103 if (useTable) 00104 CreateWindowFromTable(out); 00105 else { 00106 EWindowType type = mConfig.HasType()? 00107 mConfig.GetType() : 00108 EWindowType::eHamming; 00109 CreateTable(out,type,winsize); 00110 } 00111 00112 //zero padding is applied if audiosize is greater than window size 00113 if (winsize < audiosize) { 00114 TData* audio = out.GetPtr(); 00115 memset(audio+winsize,0,(audiosize-winsize)*sizeof(TData)); 00116 } 00117 00118 NormalizeWindow(out); 00119 if (mConfig.GetInvert()) 00120 InvertWindow(out); 00121 00122 return true; 00123 } 00124 00125 bool WindowGenerator::Do(Audio& out) 00126 { 00127 Do(out.GetBuffer()); 00128 00129 return true; 00130 } 00131 00132 00133 bool WindowGenerator::Do(Spectrum& out) 00134 { 00135 00136 CLAM_ASSERT(out.HasMagBuffer(), 00137 "WindowGenerator::Do(): Spectral Window exists only for type MagPhase"); 00138 00139 Do(out.GetMagBuffer()); 00140 return true; 00141 } 00142 00143 00144 00145 /*Create Table or window 'on the fly'*/ 00146 void WindowGenerator::CreateTable(DataArray& table,EWindowType windowType, 00147 long windowsize) const 00148 { 00149 switch(windowType)//use mathematical function according to type 00150 { 00151 case EWindowType::eKaiserBessel17: 00152 { 00153 KaiserBessel(windowsize,table,1.7); 00154 break; 00155 } 00156 case EWindowType::eKaiserBessel18: 00157 { 00158 KaiserBessel(windowsize,table,1.8); 00159 break; 00160 } 00161 case EWindowType::eKaiserBessel19: 00162 { 00163 KaiserBessel(windowsize,table,1.9); 00164 break; 00165 } 00166 case EWindowType::eKaiserBessel20: 00167 { 00168 KaiserBessel(windowsize,table,2.0); 00169 break; 00170 } 00171 case EWindowType::eKaiserBessel25: 00172 { 00173 KaiserBessel(windowsize,table,2.5); 00174 break; 00175 } 00176 case EWindowType::eKaiserBessel30: 00177 { 00178 KaiserBessel(windowsize,table,3.0); 00179 break; 00180 } 00181 case EWindowType::eKaiserBessel35: 00182 { 00183 KaiserBessel(windowsize,table,3.5); 00184 break; 00185 } 00186 case EWindowType::eBlackmanHarris62: 00187 { 00188 BlackmanHarris62(windowsize,table); 00189 break; 00190 } 00191 case EWindowType::eBlackmanHarris70: 00192 { 00193 BlackmanHarris70(windowsize,table); 00194 break; 00195 } 00196 case EWindowType::eBlackmanHarris74: 00197 { 00198 BlackmanHarris74(windowsize,table); 00199 break; 00200 } 00201 case EWindowType::eBlackmanHarris92: 00202 { 00203 BlackmanHarris92(windowsize,table); 00204 break; 00205 } 00206 case EWindowType::eHamming: 00207 { 00208 Hamming(windowsize,table); 00209 break; 00210 } 00211 case EWindowType::eTriangular: 00212 { 00213 Triangular(windowsize,table); 00214 break; 00215 } 00216 case EWindowType::eBlackmanHarris92TransMainLobe: 00217 { 00218 BlackmanHarris92TransMainLobe(windowsize,table); 00219 break; 00220 } 00221 case EWindowType::eGaussian: 00222 { 00223 Gaussian(windowsize,table); 00224 break; 00225 } 00226 case EWindowType::eBlackmanHarrisLike: 00227 { 00228 BlackmanHarrisLike(windowsize,table); 00229 break; 00230 } 00231 case EWindowType::eSine: 00232 { 00233 Sine(windowsize, table); 00234 break; 00235 } 00236 case EWindowType::eSquare: 00237 case EWindowType::eNone: 00238 { 00239 Square(windowsize, table); 00240 break; 00241 } 00242 } 00243 } 00244 00245 /*Create window from table*/ 00246 void WindowGenerator::CreateWindowFromTable(DataArray &array) const 00247 { 00248 unsigned int index = 0x00000000; 00249 unsigned int increment = (unsigned int)((double) (0x00010000) * mConfig.GetSize()/ 00250 (mSize.GetLastValue())); 00251 00252 //fix point increment [ 16bit | 16bit ] --> 1 = [ 0x0001 | 0x0000 ] 00253 00254 int size = int(mSize.GetLastValue()); 00255 CLAM_ASSERT(size<=array.Size(),"WindowGenerator::CreateWindowFromTable:output array does not have a valid size"); 00256 for (int i=0;i<size;i++) 00257 { 00258 const TData & val = mTable[index>>16]; 00259 array[i] = val; 00260 index += increment; 00261 } 00262 } 00263 00264 00265 /* function that returns the zero-order modified Bessel function of the first 00266 kind of X*/ 00267 00268 double WindowGenerator::BesselFunction(double x) const 00269 { 00270 double Sum = 0; 00271 double Factorial = 1.0; 00272 double HalfX = x/2.0; 00273 Sum += 1.0; 00274 Sum += HalfX * HalfX; 00275 for(int i=2; i<50; i++) 00276 { 00277 Factorial *= i; 00278 Sum += CLAM_pow( CLAM_pow(HalfX, (double)i) / Factorial, 2.0); 00279 } 00280 return Sum; 00281 } 00282 00283 /* function to create a Kaiser-Bessel window; window size (must be odd)*/ 00284 00285 void WindowGenerator::KaiserBessel(long size,DataArray& window, 00286 double alpha) const 00287 { 00288 TData PiAlpha = TData(PI) * TData(alpha); 00289 long windowsize = size; 00290 00291 TData dHalfsize = windowsize/2.0f; 00292 int iHalfsize = (int)windowsize/2; 00293 00294 // compute window 00295 if (windowsize % 2 != 0) 00296 window[iHalfsize]= TData(BesselFunction(PiAlpha) / BesselFunction(PiAlpha)); 00297 for(int i=0; i<iHalfsize; i++) 00298 { 00299 window[i] = window[windowsize-i-1] =TData( 00300 BesselFunction(PiAlpha * CLAM_sqrt(1.0 - CLAM_pow((double)(i-iHalfsize) / 00301 dHalfsize, 2.0))) / BesselFunction(PiAlpha) ); 00302 } 00303 00304 } 00305 00306 /* function to create a backmanHarris window*/ 00307 void WindowGenerator::BlackmanHarrisX(long size,DataArray& window, 00308 double a0,double a1,double a2,double a3) const 00309 { 00310 double fConst = TWO_PI / (size-1); 00311 /* compute window */ 00312 00313 if(size%2 !=0) 00314 { 00315 window[(int)(size/2)] = a0 - a1 * CLAM_cos(fConst * ((int)(size/2))) + a2 * 00316 CLAM_cos(fConst * 2 * ((int)(size/2))) - a3 * CLAM_cos(fConst * 3 * ((int)(size/2))); 00317 } 00318 for(int i = 0; i < (int)(size/2); i++) 00319 { 00320 window[i] = window[size-i-1] = a0 - a1 * CLAM_cos(fConst * i) + 00321 a2 * CLAM_cos(fConst * 2 * i) - a3 * CLAM_cos(fConst * 3 * i); 00322 } 00323 00324 00325 00326 } 00327 00328 00329 /* function to create a BlackmanHarris window*/ 00330 void WindowGenerator::BlackmanHarris62(long size,DataArray& window) const 00331 { 00332 /* for 3 term -62.05 */ 00333 double a0 = .44959, a1 = .49364, a2 = .05677; 00334 BlackmanHarrisX(size,window,a0,a1,a2); 00335 } 00336 00337 00338 /* function to create a backmanHarris window*/ 00339 00340 void WindowGenerator::BlackmanHarris70(long size,DataArray& window) const 00341 { 00342 /* for 3 term -70.83 */ 00343 double a0 = .42323, a1 = .49755, a2 = .07922; 00344 BlackmanHarrisX(size,window,a0,a1,a2); 00345 } 00346 00347 /* function to create a backmanHarris window*/ 00348 void WindowGenerator::BlackmanHarris74(long size,DataArray& window) const 00349 { 00350 00351 /* for -74dB from the Nuttall paper */ 00352 double a0 = .40217, a1 = .49703, a2 = .09892, a3 = .00188; 00353 00354 BlackmanHarrisX(size,window,a0,a1,a2,a3); 00355 } 00356 00357 /* function to create a backmanHarris window*/ 00358 void WindowGenerator::BlackmanHarris92(long size,DataArray& window) const 00359 { 00360 00361 /* for -92dB */ 00362 double a0 = .35875, a1 = .48829, a2 = .14128, a3 = .01168; 00363 00364 BlackmanHarrisX(size,window,a0,a1,a2,a3); 00365 } 00366 00367 void WindowGenerator::BlackmanHarrisLike(long size, DataArray& window) const 00368 { 00369 TData fSum=0; 00370 // float a0 = .51f, a1 = .42f, a2 = -0.04f, a3 = .03f, a4=0.03f, a5=0.05f; 00371 for(int i=0; i<size; i++) 00372 fSum += window[i] = 00373 + 0.47 00374 - 0.45 * CLAM_cos(TData(TWO_PI/(size-1.0)*i)) 00375 - 0.01 * CLAM_cos(TData(TWO_PI/(size-1.0)*i*2.0)) 00376 - 0.01 * CLAM_cos(TData(TWO_PI/(size-1.0)*i*3.0)); 00377 fSum = fSum/2; 00378 for (int i = 0; i < size; i++) 00379 window[i] = window[i] / fSum; 00380 return; 00381 } 00382 00383 00384 /* function to design a Hamming window*/ 00385 void WindowGenerator::Hamming(long size,DataArray& window) const 00386 { 00387 if(size%2 !=0) 00388 window[(int)(size/2)]= 00389 + TData(0.54) 00390 - TData(0.46)*CLAM_cos(TData(2.0)*TData(PI)*((int)(size/2))/(size-1)); 00391 for(int i = 0; i < (int)(size/2); i++) 00392 window[i] = window[size-i-1] = 00393 + TData(0.54) 00394 - TData(0.46) * CLAM_cos(TData(2.0)*TData(PI)*i/(size-1)); 00395 } 00396 00397 /* function to design a Triangular window*/ 00398 void WindowGenerator::Triangular(long size,DataArray& window) const 00399 { 00400 if(size%2 !=0) 00401 window[(int)(size/2)] = (TData)2*((int)(size/2))/(size-1); 00402 for(int i = 0; i < (int)(size/2); i++) 00403 { 00404 window[i] = window[size-i-1]= (TData)2*i/(size-1); 00405 } 00406 } 00407 00408 /* function to create the (fft-)transformed Mainlobe of a BlackmanHarris 92 dB window*/ 00409 void WindowGenerator::BlackmanHarris92TransMainLobe(long size,DataArray& window) const 00410 { 00411 const short N = 512; 00412 TData fA[4] = {TData(.35875), TData(.48829), TData(.14128), TData(.01168)}, 00413 fMax = 0; 00414 TData fTheta = -TData(4.0) * TData(TWO_PI) / N, 00415 fThetaIncr = (TData(8.0) * TData(TWO_PI) / N) / (size); 00416 00417 00418 for(int i = 0; i < size; i++) 00419 { 00420 window[i] = 0; // init value 00421 for (int m = 0; m < 4; m++) 00422 window[i] += -1 * (fA[m]/2) * 00423 (Sinc (fTheta - m * TData(TWO_PI)/N, N) + 00424 Sinc (fTheta + m * TData(TWO_PI)/N, N)); 00425 fTheta += fThetaIncr; 00426 } 00427 00428 /* normalize window */ 00429 fMax = window[(int) size / 2]; 00430 for (int i = 0; i < size; i++) 00431 window[i] = window[i] / fMax; 00432 00433 } 00434 00435 void WindowGenerator::Gaussian(long size,DataArray& window) const 00436 { 00437 double s = 0.15; 00438 double scale = 1.0 / ( 2.0 * M_PI * 0.15 * 0.15 ); 00439 00440 if(size%2 !=0) 00441 window[size/2] = scale; 00442 for(int i = 0; i < size/2; i++) 00443 { 00444 TData x = (TData)(i-(TData)size/2.)/(TData)(size-1); 00445 window[i] = window[size-i-1]= scale * exp(-(x*x)/(2*s*s)); 00446 } 00447 } 00448 00449 /* function to design a Sine window*/ 00450 void WindowGenerator::Sine(long size,DataArray& window) const 00451 { 00452 double tmp1 = PI/(2.0*float(size)); 00453 double tmp2 = 0.5*(2.0*float(size)); 00454 00455 for (int i=0;i<size;i++) 00456 window[i] = (float)(1+tmp2*CLAM_sin(tmp1*(i+1))); 00457 } 00458 void WindowGenerator::Square(long size,DataArray& window) const 00459 { 00460 for (int i=0;i<size;i++) 00461 window[i] = 1; 00462 } 00463 00464 00465 void WindowGenerator::InvertWindow(const DataArray& originalWindow, 00466 DataArray& invertedWindow) const 00467 { 00468 if(invertedWindow.AllocatedSize()!=originalWindow.AllocatedSize()) 00469 invertedWindow.Resize(originalWindow.AllocatedSize()); 00470 if(invertedWindow.Size()!=originalWindow.Size()) 00471 invertedWindow.SetSize(originalWindow.Size()); 00472 00473 if (originalWindow.Size()%2!=0) 00474 if(originalWindow[(int)(originalWindow.Size()/2)]!=0) 00475 invertedWindow[(int)(originalWindow.Size()/2)]= 00476 1/originalWindow[(int)(originalWindow.Size()/2)]; 00477 for(int i=0;i<(int)(originalWindow.Size()/2);i++) 00478 { 00479 if(originalWindow[i]!=0) 00480 invertedWindow[i]=invertedWindow[originalWindow.Size()-1-i]=1.0/originalWindow[i]; 00481 } 00482 invertedWindow[originalWindow.Size()-1]=0; 00483 } 00484 00485 void WindowGenerator::InvertWindow(DataArray& window) const 00486 { 00487 InvertWindow(window,window); 00488 } 00489 00490 void WindowGenerator::NormalizeWindow(DataArray& window) const 00491 { 00492 if (mConfig.GetNormalize() == EWindowNormalize::eNone) return; 00493 double normalizationFactor=1.0; 00494 00495 const int size = window.Size(); 00496 //Note: We multiply by two because we add the energy of the negative spectrum 00497 switch (mConfig.GetNormalize()) { 00498 case EWindowNormalize::eAnalysis: 00499 { 00500 double sum=0.0; 00501 for(int i=0;i<size;i++) sum+=window[i]; 00502 normalizationFactor=1/(sum/2); 00503 break; 00504 } 00505 case EWindowNormalize::eEnergy: 00506 { 00507 double sum=0.0; 00508 for(int i=0;i<size;i++) sum+=window[i]; 00509 normalizationFactor=size/(sum/2); 00510 break; 00511 } 00512 case EWindowNormalize::eMax: 00513 { 00514 double max=0.0; 00515 for(int i=0;i<size;i++) 00516 if (max<window[i]) max=window[i]; 00517 normalizationFactor = 1.0/max; // be careful with even windows!! 00518 break; 00519 } 00520 default: 00521 CLAM_ASSERT(false, "Unexpected normalization type"); 00522 } 00523 00524 for(int i=0;i<size;i++) 00525 { 00526 window[i]*=normalizationFactor; 00527 } 00528 } 00529 00530 /* internal math functions */ 00531 double WindowGenerator::Sinc(double x, short N) const 00532 { 00533 return (sin ((N/2) * x) / sin (x/2)); 00534 } 00535 00536 } 00537