Mercurial > hg > x

comparison WindowFunctions.cpp @ 1:6422640a802f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first upload
author Wen X <xue.wen@elec.qmul.ac.uk>
date Tue, 05 Oct 2010 10:45:57 +0100
parents
children
comparison
equal deleted inserted replaced
0:9b9f21935f24 1:6422640a802f
1 //---------------------------------------------------------------------------
2 #include <math.h>
3 #include <mem.h>
4 #include "WindowFunctions.h"
5 #include "align8.h"
6 //---------------------------------------------------------------------------
7
8 /*
9 function I0: Bessel function of order zero
10
11 Returns the Bessel function of x.
12 */
13 double I0(double x)
14 {
15 const double eps=1e-9;
16
17 int n=1;
18 double S=1, D=1, T;
19
20 while (D>eps*S)
21 {
22 T=x/(2*n++);
23 D*=T*T;
24 S+=D;
25 }
26 return S;
27 }//I0
28
29 /*
30 function FillWindow: fills a buffer $Result with a window function.
31
32 In: wt:window type
33 Count: window size
34 ips & ups: extra window specificatin parameters
35 Out: newwindow[Count+1]: the window function.
36
37 No return value. This function is designed assuming Count being even, and Count/2 is the symmetry
38 centre of newwindow. newwindow has physical size Count+1 for compatibility purpose. For all vanishing
39 windows (e.g. Hann) newwindow[0]=newwindow[Count]=0.
40 */
41 void FillWindow(double* Result, WindowType wt, int Count, int* ips, double* dps)
42 {
43 switch(wt)
44 {
45 case wtRectangle:
46 for (int i=0; i<=Count; i++)
47 Result[i]=1;
48 break;
49 case wtTriangular:
50 for (int i=0; i<=Count; i++)
51 Result[i]=1-2*fabs(i-Count/2.0)/Count;
52 break;
53 case wtHamming:
54 for (int i=0; i<=Count; i++)
55 Result[i]=(0.54-0.46*cos(2*M_PI*i/Count));
56 break;
57 case wtBlackman:
58 for (int i=0; i<=Count; i++)
59 Result[i]=0.42-0.5*cos(2*M_PI*i/Count)+0.08*cos(4*M_PI*i/Count);
60 break;
61 case wtGaussian:
62 {
63 double num=*dps*M_PI*M_PI/2/Count/Count;
64 for (int i=0; i<=Count; i++)
65 Result[i]=exp(-(i-Count/2.0)*(i-Count/2.0)*num);
66 break;
67 }
68 case wtKaiser:
69 {
70 double ldps=*dps*M_PI*M_PI/2/Count;
71 double den=I0(0.5*Count*ldps);
72 for (int i=0; i<=Count; i++)
73 Result[i]=I0(ldps*sqrt(0.25*Count*Count-(i-Count*0.5)*(i-Count*0.5)))/den;
74 break;
75 }
76 case wtHalfCosine:
77 {
78 for (int i=0; i<=Count; i++)
79 Result[i]=sin(M_PI*i/Count);
80 break;
81 }
82 case wtHann:
83 {
84 for (int i=0; i<=Count; i++)
85 Result[i]=(0.5-cos(M_PI*2*i/Count)/2);
86 break;
87 }
88 case wtHannSqr:
89 {
90 for (int i=0; i<=Count; i++)
91 Result[i]=(3-4*cos(M_PI*2*i/Count)+cos(M_PI*4*i/Count))/8;
92 break;
93 }
94 case wtHann3sqrt:
95 {
96 for (int i=0; i<=Count; i++)
97 {
98 double tmp=sin(M_PI*i/Count);
99 Result[i]=tmp*tmp*tmp;
100 }
101 break;
102 }
103 }
104 }//FillWindow
105
106 /*
107 function NewWindow: creates a window function.
108
109 In: wt: window type
110 Count: window size
111 ips & ups: extra window specificatin parameters
112 Out: newwindow[Count+1]: the window function.
113
114 Returns pointer to newwindow. newwindow is created anew if newwindow=0 is specified on start.
115 */
116 double* NewWindow(WindowType wt, int Count, int* ips, double* dps, double* newwindow)
117 {
118 double* Result=newwindow; if (!Result) Result=new double[Count+1];
119 FillWindow(Result, wt, Count, ips, dps);
120 return Result;
121 }//NewWindow
122
123 /*
124 function NewWindow8: 8-byte aligned version of NewWindow.
125
126 In: wt: window type
127 Count: window size
128 ips & ups: extra window specificatin parameters
129 Out: newwindow[Count+1]: the window function.
130
131 Returns pointer to newwindow. newwindow is created anew and 8-byte aligned if newwindow=0 is
132 specified on start. However, if newwindow is not 8-byte aligned on start, the unaligned buffer is
133 returned.
134 */
135 double* NewWindow8(WindowType wt, int Count, int* ips, double* dps, double* newwindow)
136 {
137 double* Result=newwindow; if (!Result) Result=(double*)malloc8(sizeof(double)*(Count+1));
138 FillWindow(Result, wt, Count, ips, dps);
139 return Result;
140 }//NewWindow8
141
142 /*
143 function NewdWindow: computes the derivative of a window function.
144
145 In: wt: window type
146 Count: window size
147 ips & ups: extra window specificatin parameters
148 Out: newdwindow[Count+1]: the derivative window function.
149
150 Returns pointer to newdwindow. newdwindow is created anew if newdwindow=0 is specified on start.
151 */
152 double* NewdWindow(WindowType wt, int Count, int* ips, double* dps, double* newdwindow)
153 {
154 double* Result=newdwindow; if (!Result) Result=new double[Count+1];
155 double piiCount=M_PI/Count;
156 switch(wt)
157 {
158 case wtRectangle:
159 memset(Result, 0, sizeof(double)*(Count+1));
160 break;
161 case wtTriangular:
162 throw("Trying to differentiate triangular window.");
163 break;
164 case wtHamming:
165 for (int i=0; i<=Count; i++)
166 Result[i]=0.92*piiCount*sin(2*piiCount*i);
167 break;
168 case wtBlackman:
169 for (int i=0; i<=Count; i++)
170 Result[i]=piiCount*sin(2*piiCount*i)-0.32*piiCount*sin(4*piiCount*i);
171 break;
172 case wtGaussian:
173 throw("Trying to differentiate Gaussian window.");
174 break;
175 case wtKaiser:
176 throw("Trying to differentiate Kaiser window.");
177 break;
178 case wtHalfCosine:
179 {
180 for (int i=0; i<=Count; i++)
181 Result[i]=piiCount*cos(piiCount*i);
182 break;
183 }
184 case wtHann:
185 {
186 for (int i=0; i<=Count; i++)
187 Result[i]=piiCount*sin(2*piiCount*i);
188 break;
189 }
190 case wtHannSqr:
191 {
192 for (int i=0; i<=Count; i++)
193 Result[i]=piiCount*sin(2*piiCount*i)-0.5*piiCount*sin(4*piiCount*i);
194 break;
195 }
196 case wtHann3sqrt:
197 {
198 for (int i=0; i<=Count; i++)
199 {
200 double s=sin(M_PI*i/Count), c=cos(M_PI*i/Count);
201 Result[i]=3*piiCount*s*s*c;
202 }
203 }
204 }
205 return Result;
206 }//NewdWindow
207
208 /*
209 function NewddWindow: computes the 2nd-order derivative of a window function.
210
211 In: wt: window type
212 Count: window size
213 ips & ups: extra window specificatin parameters
214 Out: newddwindow[Count+1]: the 2nd-order derivative window function.
215
216 Returns pointer to newddwindow. newddwindow is created anew if newddwindow=0 is specified on start.
217 */
218 double* NewddWindow(WindowType wt, int Count, int* ips, double* dps, double* newddwindow)
219 {
220 double* Result=newddwindow; if (!Result) Result=new double[Count+1];
221 double piiCount=M_PI/Count;
222 double piC2=piiCount*piiCount;
223 switch(wt)
224 {
225 case wtRectangle:
226 memset(Result, 0, sizeof(double)*(Count+1));
227 break;
228 case wtTriangular:
229 throw("Trying to double-differentiate triangular window.");
230 break;
231 case wtHamming:
232 for (int i=0; i<=Count; i++)
233 Result[i]=1.84*piC2*cos(2*piiCount*i);
234 break;
235 case wtBlackman:
236 for (int i=0; i<=Count; i++)
237 Result[i]=2*piC2*cos(2*piiCount*i)-1.28*piC2*cos(4*piiCount*i);
238 break;
239 case wtGaussian:
240 throw("Trying to double-differentiate Gaussian window.");
241 break;
242 case wtKaiser:
243 throw("Trying to double-differentiate Kaiser window.");
244 break;
245 case wtHalfCosine:
246 {
247 for (int i=0; i<=Count; i++)
248 Result[i]=-piC2*sin(piiCount*i);
249 break;
250 }
251 case wtHann:
252 {
253 for (int i=0; i<=Count; i++)
254 Result[i]=2*piC2*cos(2*piiCount*i);
255 break;
256 }
257 case wtHannSqr:
258 {
259 for (int i=0; i<=Count; i++)
260 Result[i]=2*piC2*cos(2*piiCount*i)-2*piC2*cos(4*piiCount*i);
261 break;
262 }
263 case wtHann3sqrt:
264 {
265 for (int i=0; i<=Count; i++)
266 {
267 double s=sin(M_PI*i/Count), c=cos(M_PI*i/Count);
268 Result[i]=3*piC2*(2*c*c-s*s)*s;
269 }
270 break;
271 }
272 }
273 return Result;
274 }//NewddWindow
275
276 /*
277 function NewdddWindow: computes the 3rd-order derivative of a window function.
278 In: wt: window type
279 Count: window size
280 ips & ups: extra window specificatin parameters
281 Out: newdddwindow[Count+1]: the 3rd-order derivative window function.
282
283 Returns pointer to newdddwindow. newdddwindow is created anew if newdddwindow=0 is specified on start.
284 */
285 double* NewdddWindow(WindowType wt, int Count, int* ips, double* dps, double* newdddwindow)
286 {
287 double* Result=newdddwindow; if (!Result) Result=new double[Count+1];
288 double piiCount=M_PI/Count;
289 double piC2=piiCount*piiCount;
290 double piC3=piiCount*piC2;
291 switch(wt)
292 {
293 case wtRectangle:
294 memset(Result, 0, sizeof(double)*(Count+1));
295 break;
296 case wtTriangular:
297 throw("Trying to triple-differentiate triangular window.");
298 break;
299 case wtHamming:
300 for (int i=0; i<=Count; i++)
301 Result[i]=-3.68*piC3*sin(2*piiCount*i);
302 break;
303 case wtBlackman:
304 for (int i=0; i<=Count; i++)
305 Result[i]=-4*piC3*sin(2*piiCount*i)+5.12*piC3*sin(4*piiCount*i);
306 break;
307 case wtGaussian:
308 throw("Trying to triple-differentiate Gaussian window.");
309 break;
310 case wtKaiser:
311 throw("Trying to triple-differentiate Kaiser window.");
312 break;
313 case wtHalfCosine:
314 {
315 for (int i=0; i<=Count; i++)
316 Result[i]=-piC3*cos(piiCount*i);
317 break;
318 }
319 case wtHann:
320 {
321 for (int i=0; i<=Count; i++)
322 Result[i]=-4*piC3*sin(2*piiCount*i);
323 break;
324 }
325 case wtHannSqr:
326 {
327 for (int i=0; i<=Count; i++)
328 Result[i]=-4*piC3*sin(2*piiCount*i)+8*piC3*sin(4*piiCount*i);
329 break;
330 }
331 case wtHann3sqrt:
332 throw("Trying to triple-differentiate Hann^1.5 window.");
333 break;
334 }
335 return Result;
336 }//NewdddWindow
337
338 //---------------------------------------------------------------------------
339 /*
340 function windowspec: computes a few descriptors of a cosine family window. A window function in the
341 cosine window family is the linear combination of a few cosine functions, therefore has a cosine
342 decomposition.
343
344 In: wt: window type
345 N: window size
346 Out: c[M+1], coefficients of cosine decomposition
347 iH2: reciprocal of square of L2 norm,
348 d[2M+1]: self-convolution of c[], optional
349
350 No return value.
351 */
352 void windowspec(WindowType wt, int N, int* M, double* c, double* iH2, double* d)
353 {
354 switch(wt)
355 {
356 case wtRectangle:
357 M[0]=0, c[0]=1, iH2[0]=1.0/N/N;
358 if (d) d[0]=1;
359 break;
360 case wtHamming:
361 M[0]=1, c[0]=0.54, c[1]=0.23, iH2[0]=1.0/(0.3974*N*N);
362 if (d) d[0]=0.3974, d[1]=0.2484, d[2]=0.0529;
363 break;
364 case wtBlackman:
365 M[0]=2, c[0]=0.42, c[1]=0.25, c[2]=0.04, iH2[0]=1.0/(0.3046*N*N);
366 if (d) d[0]=0.3046, d[1]=0.23, d[2]=0.0961, d[4]=0.02, d[5]=0.0016;
367 break;
368 case wtHann:
369 M[0]=1, c[0]=0.5, c[1]=0.25, iH2[0]=1.0/(0.375*N*N);
370 if (d) d[0]=0.375, d[1]=0.25, d[2]=0.0625;
371 break;
372 default:
373 M[0]=1, c[0]=0.5, c[1]=0.25, iH2[0]=1.0/(0.375*N*N);
374 if (d) d[0]=0.375, d[1]=0.25, d[2]=0.0625;
375 break;
376 }
377 }//windowspec