CLAM-Development  1.4.0
WindowGenerator.cxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program 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 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program 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 this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 #include "ProcessingData.hxx"
23 #include "DataTypes.hxx"
24 #include "Spectrum.hxx"
25 #include "Audio.hxx"
26 #include "WindowGenerator.hxx"
27 #include "ProcessingFactory.hxx"
28 
29 namespace CLAM
30 {
31 
32 namespace Hidden
33 {
34  static const char * metadata[] = {
35  "key", "WindowGenerator",
36  "category", "Generators", // must change to Analysis?
37  "description", "WindowGenerator",
38  0
39  };
40  static FactoryRegistrator<ProcessingFactory, WindowGenerator> reg = metadata;
41 }
42 
43 /* Processing object Method implementations */
44 
46  : mOutput( "Generated window function samples", this )
47  , mSize("Size",this)
48 {
49  Configure(c);
50 }
51 
53 {}
54 
55 
56 /* Configure the Processing Object according to the Config object */
57 
58 bool WindowGenerator::ConcreteConfigure(const ProcessingConfig& c)
59 {
60  CopyAsConcreteConfig(mConfig, c);
61  mSize.DoControl(TControlData(mConfig.GetSize()));
62 
63  if (!mConfig.HasUseTable()) return true;
64  if (!mConfig.GetUseTable()) return true;
65 
66 
67  /* Fill the table */
68 
69  mTable.Resize(mConfig.GetSize());
70  mTable.SetSize(mConfig.GetSize());
71  mSize.DoControl(TControlData(mConfig.GetSize()));
72 
73  EWindowType type = mConfig.HasType()?
74  mConfig.GetType() :
76 
77  CreateTable(mTable,type,mConfig.GetSize());
78 
79  return true;
80 
81 }
82 
83 /* The supervised Do() function */
84 
86 {
87  CLAM_ASSERT( AbleToExecute(), "This processing is not ready to do anything" );
88 
89  bool result = Do( mOutput.GetAudio() );
90  mOutput.Produce();
91  return result;
92 }
93 
94 /* The unsupervised Do() function */
95 
97 {
98  const int winsize = (int) mSize.GetLastValue();
99  const int audiosize = out.Size();
100 
101  bool useTable = mConfig.HasUseTable() && mConfig.GetUseTable();
102 
103  if (useTable)
104  CreateWindowFromTable(out);
105  else {
106  EWindowType type = mConfig.HasType()?
107  mConfig.GetType() :
109  CreateTable(out,type,winsize);
110  }
111 
112  //zero padding is applied if audiosize is greater than window size
113  if (winsize < audiosize) {
114  TData* audio = out.GetPtr();
115  memset(audio+winsize,0,(audiosize-winsize)*sizeof(TData));
116  }
117 
118  NormalizeWindow(out);
119  if (mConfig.GetInvert())
120  InvertWindow(out);
121 
122  return true;
123 }
124 
126 {
127  Do(out.GetBuffer());
128 
129  return true;
130 }
131 
132 
134 {
135 
136  CLAM_ASSERT(out.HasMagBuffer(),
137  "WindowGenerator::Do(): Spectral Window exists only for type MagPhase");
138 
139  Do(out.GetMagBuffer());
140  return true;
141 }
142 
143 
144 
145 /*Create Table or window 'on the fly'*/
146 void WindowGenerator::CreateTable(DataArray& table,EWindowType windowType,
147  long windowsize) const
148 {
149  switch(windowType)//use mathematical function according to type
150  {
152  {
153  KaiserBessel(windowsize,table,1.7);
154  break;
155  }
157  {
158  KaiserBessel(windowsize,table,1.8);
159  break;
160  }
162  {
163  KaiserBessel(windowsize,table,1.9);
164  break;
165  }
167  {
168  KaiserBessel(windowsize,table,2.0);
169  break;
170  }
172  {
173  KaiserBessel(windowsize,table,2.5);
174  break;
175  }
177  {
178  KaiserBessel(windowsize,table,3.0);
179  break;
180  }
182  {
183  KaiserBessel(windowsize,table,3.5);
184  break;
185  }
187  {
188  BlackmanHarris62(windowsize,table);
189  break;
190  }
192  {
193  BlackmanHarris70(windowsize,table);
194  break;
195  }
197  {
198  BlackmanHarris74(windowsize,table);
199  break;
200  }
202  {
203  BlackmanHarris92(windowsize,table);
204  break;
205  }
207  {
208  Hamming(windowsize,table);
209  break;
210  }
212  {
213  Triangular(windowsize,table);
214  break;
215  }
217  {
218  BlackmanHarris92TransMainLobe(windowsize,table);
219  break;
220  }
222  {
223  Gaussian(windowsize,table);
224  break;
225  }
227  {
228  BlackmanHarrisLike(windowsize,table);
229  break;
230  }
231  case EWindowType::eSine:
232  {
233  Sine(windowsize, table);
234  break;
235  }
237  case EWindowType::eNone:
238  {
239  Square(windowsize, table);
240  break;
241  }
242  }
243 }
244 
245 /*Create window from table*/
246 void WindowGenerator::CreateWindowFromTable(DataArray &array) const
247 {
248  unsigned int index = 0x00000000;
249  unsigned int increment = (unsigned int)((double) (0x00010000) * mConfig.GetSize()/
250  (mSize.GetLastValue()));
251 
252  //fix point increment [ 16bit | 16bit ] --> 1 = [ 0x0001 | 0x0000 ]
253 
254  int size = int(mSize.GetLastValue());
255  CLAM_ASSERT(size<=array.Size(),"WindowGenerator::CreateWindowFromTable:output array does not have a valid size");
256  for (int i=0;i<size;i++)
257  {
258  const TData & val = mTable[index>>16];
259  array[i] = val;
260  index += increment;
261  }
262 }
263 
264 
265 /* function that returns the zero-order modified Bessel function of the first
266 kind of X*/
267 
268 double WindowGenerator::BesselFunction(double x) const
269 {
270  double Sum = 0;
271  double Factorial = 1.0;
272  double HalfX = x/2.0;
273  Sum += 1.0;
274  Sum += HalfX * HalfX;
275  for(int i=2; i<50; i++)
276  {
277  Factorial *= i;
278  Sum += CLAM_pow( CLAM_pow(HalfX, (double)i) / Factorial, 2.0);
279  }
280  return Sum;
281 }
282 
283 /* function to create a Kaiser-Bessel window; window size (must be odd)*/
284 
285 void WindowGenerator::KaiserBessel(long size,DataArray& window,
286  double alpha) const
287 {
288  TData PiAlpha = TData(PI) * TData(alpha);
289  long windowsize = size;
290 
291  TData dHalfsize = windowsize/2.0f;
292  int iHalfsize = (int)windowsize/2;
293 
294  // compute window
295  if (windowsize % 2 != 0)
296  window[iHalfsize]= TData(BesselFunction(PiAlpha) / BesselFunction(PiAlpha));
297  for(int i=0; i<iHalfsize; i++)
298  {
299  window[i] = window[windowsize-i-1] =TData(
300  BesselFunction(PiAlpha * CLAM_sqrt(1.0 - CLAM_pow((double)(i-iHalfsize) /
301  dHalfsize, 2.0))) / BesselFunction(PiAlpha) );
302  }
303 
304 }
305 
306 /* function to create a backmanHarris window*/
307 void WindowGenerator::BlackmanHarrisX(long size,DataArray& window,
308  double a0,double a1,double a2,double a3) const
309 {
310  double fConst = TWO_PI / (size-1);
311  /* compute window */
312 
313  if(size%2 !=0)
314  {
315  window[(int)(size/2)] = a0 - a1 * CLAM_cos(fConst * ((int)(size/2))) + a2 *
316  CLAM_cos(fConst * 2 * ((int)(size/2))) - a3 * CLAM_cos(fConst * 3 * ((int)(size/2)));
317  }
318  for(int i = 0; i < (int)(size/2); i++)
319  {
320  window[i] = window[size-i-1] = a0 - a1 * CLAM_cos(fConst * i) +
321  a2 * CLAM_cos(fConst * 2 * i) - a3 * CLAM_cos(fConst * 3 * i);
322  }
323 
324 
325 
326 }
327 
328 
329 /* function to create a BlackmanHarris window*/
330 void WindowGenerator::BlackmanHarris62(long size,DataArray& window) const
331 {
332  /* for 3 term -62.05 */
333  double a0 = .44959, a1 = .49364, a2 = .05677;
334  BlackmanHarrisX(size,window,a0,a1,a2);
335 }
336 
337 
338 /* function to create a backmanHarris window*/
339 
340 void WindowGenerator::BlackmanHarris70(long size,DataArray& window) const
341 {
342  /* for 3 term -70.83 */
343  double a0 = .42323, a1 = .49755, a2 = .07922;
344  BlackmanHarrisX(size,window,a0,a1,a2);
345 }
346 
347 /* function to create a backmanHarris window*/
348 void WindowGenerator::BlackmanHarris74(long size,DataArray& window) const
349 {
350 
351  /* for -74dB from the Nuttall paper */
352  double a0 = .40217, a1 = .49703, a2 = .09892, a3 = .00188;
353 
354  BlackmanHarrisX(size,window,a0,a1,a2,a3);
355 }
356 
357 /* function to create a backmanHarris window*/
358 void WindowGenerator::BlackmanHarris92(long size,DataArray& window) const
359 {
360 
361  /* for -92dB */
362  double a0 = .35875, a1 = .48829, a2 = .14128, a3 = .01168;
363 
364  BlackmanHarrisX(size,window,a0,a1,a2,a3);
365 }
366 
367 void WindowGenerator::BlackmanHarrisLike(long size, DataArray& window) const
368 {
369  TData fSum=0;
370 // float a0 = .51f, a1 = .42f, a2 = -0.04f, a3 = .03f, a4=0.03f, a5=0.05f;
371  for(int i=0; i<size; i++)
372  fSum += window[i] =
373  + 0.47
374  - 0.45 * CLAM_cos(TData(TWO_PI/(size-1.0)*i))
375  - 0.01 * CLAM_cos(TData(TWO_PI/(size-1.0)*i*2.0))
376  - 0.01 * CLAM_cos(TData(TWO_PI/(size-1.0)*i*3.0));
377  fSum = fSum/2;
378  for (int i = 0; i < size; i++)
379  window[i] = window[i] / fSum;
380  return;
381 }
382 
383 
384 /* function to design a Hamming window*/
385 void WindowGenerator::Hamming(long size,DataArray& window) const
386 {
387  if(size%2 !=0)
388  window[(int)(size/2)]=
389  + TData(0.54)
390  - TData(0.46)*CLAM_cos(TData(2.0)*TData(PI)*((int)(size/2))/(size-1));
391  for(int i = 0; i < (int)(size/2); i++)
392  window[i] = window[size-i-1] =
393  + TData(0.54)
394  - TData(0.46) * CLAM_cos(TData(2.0)*TData(PI)*i/(size-1));
395 }
396 
397 /* function to design a Triangular window*/
398 void WindowGenerator::Triangular(long size,DataArray& window) const
399 {
400  if(size%2 !=0)
401  window[(int)(size/2)] = (TData)2*((int)(size/2))/(size-1);
402  for(int i = 0; i < (int)(size/2); i++)
403  {
404  window[i] = window[size-i-1]= (TData)2*i/(size-1);
405  }
406 }
407 
408 /* function to create the (fft-)transformed Mainlobe of a BlackmanHarris 92 dB window*/
409 void WindowGenerator::BlackmanHarris92TransMainLobe(long size,DataArray& window) const
410 {
411  const short N = 512;
412  TData fA[4] = {TData(.35875), TData(.48829), TData(.14128), TData(.01168)},
413  fMax = 0;
414  TData fTheta = -TData(4.0) * TData(TWO_PI) / N,
415  fThetaIncr = (TData(8.0) * TData(TWO_PI) / N) / (size);
416 
417 
418  for(int i = 0; i < size; i++)
419  {
420  window[i] = 0; // init value
421  for (int m = 0; m < 4; m++)
422  window[i] += -1 * (fA[m]/2) *
423  (Sinc (fTheta - m * TData(TWO_PI)/N, N) +
424  Sinc (fTheta + m * TData(TWO_PI)/N, N));
425  fTheta += fThetaIncr;
426  }
427 
428  /* normalize window */
429  fMax = window[(int) size / 2];
430  for (int i = 0; i < size; i++)
431  window[i] = window[i] / fMax;
432 
433 }
434 
435 void WindowGenerator::Gaussian(long size,DataArray& window) const
436 {
437  double s = 0.15;
438  double scale = 1.0 / ( 2.0 * M_PI * 0.15 * 0.15 );
439 
440  if(size%2 !=0)
441  window[size/2] = scale;
442  for(int i = 0; i < size/2; i++)
443  {
444  TData x = (TData)(i-(TData)size/2.)/(TData)(size-1);
445  window[i] = window[size-i-1]= scale * exp(-(x*x)/(2*s*s));
446  }
447 }
448 
449 /* function to design a Sine window*/
450 void WindowGenerator::Sine(long size,DataArray& window) const
451 {
452  double tmp1 = PI/(2.0*float(size));
453  double tmp2 = 0.5*(2.0*float(size));
454 
455  for (int i=0;i<size;i++)
456  window[i] = (float)(1+tmp2*CLAM_sin(tmp1*(i+1)));
457 }
458 void WindowGenerator::Square(long size,DataArray& window) const
459 {
460  for (int i=0;i<size;i++)
461  window[i] = 1;
462 }
463 
464 
465 void WindowGenerator::InvertWindow(const DataArray& originalWindow,
466  DataArray& invertedWindow) const
467 {
468  if(invertedWindow.AllocatedSize()!=originalWindow.AllocatedSize())
469  invertedWindow.Resize(originalWindow.AllocatedSize());
470  if(invertedWindow.Size()!=originalWindow.Size())
471  invertedWindow.SetSize(originalWindow.Size());
472 
473  if (originalWindow.Size()%2!=0)
474  if(originalWindow[(int)(originalWindow.Size()/2)]!=0)
475  invertedWindow[(int)(originalWindow.Size()/2)]=
476  1/originalWindow[(int)(originalWindow.Size()/2)];
477  for(int i=0;i<(int)(originalWindow.Size()/2);i++)
478  {
479  if(originalWindow[i]!=0)
480  invertedWindow[i]=invertedWindow[originalWindow.Size()-1-i]=1.0/originalWindow[i];
481  }
482  invertedWindow[originalWindow.Size()-1]=0;
483 }
484 
485 void WindowGenerator::InvertWindow(DataArray& window) const
486 {
487  InvertWindow(window,window);
488 }
489 
490 void WindowGenerator::NormalizeWindow(DataArray& window) const
491 {
492  if (mConfig.GetNormalize() == EWindowNormalize::eNone) return;
493  double normalizationFactor=1.0;
494 
495  const int size = window.Size();
496  //Note: We multiply by two because we add the energy of the negative spectrum
497  switch (mConfig.GetNormalize()) {
499  {
500  double sum=0.0;
501  for(int i=0;i<size;i++) sum+=window[i];
502  normalizationFactor=1/(sum/2);
503  break;
504  }
506  {
507  double sum=0.0;
508  for(int i=0;i<size;i++) sum+=window[i];
509  normalizationFactor=size/(sum/2);
510  break;
511  }
513  {
514  double max=0.0;
515  for(int i=0;i<size;i++)
516  if (max<window[i]) max=window[i];
517  normalizationFactor = 1.0/max; // be careful with even windows!!
518  break;
519  }
520  default:
521  CLAM_ASSERT(false, "Unexpected normalization type");
522  }
523 
524  for(int i=0;i<size;i++)
525  {
526  window[i]*=normalizationFactor;
527  }
528 }
529 
530 /* internal math functions */
531  double WindowGenerator::Sinc(double x, short N) const
532 {
533  return (sin ((N/2) * x) / sin (x/2));
534 }
535 
536 }
537