一般的讲数字信号处理的书中都会提到窗函数。大多数只会提及其中的几种。这里我把这些窗都用C语言实现了一下,都不复杂,但如果要自己去弄也挺费时间。所有函数都用Matlab验证了。包括以下窗:
[url=]file:///C:/Users/zhang/AppData/Local/Temp/msohtmlclip1/01/clip_image001.gif[/url]
1 /*窗类型*/
2 typedef enum
3 {
4 Bartlett = 0,
5 BartLettHann,
6 BlackMan,
7 BlackManHarris,
8 Bohman,
9 Chebyshev,
10 FlatTop,
11 Gaussian,
12 Hamming,
13 Hann,
14 Kaiser,
15 Nuttal,
16 Parzen,
17 Rectangular,
18 Taylor,
19 Triangular,
20 Tukey
21 }winType;
[url=]file:///C:/Users/zhang/AppData/Local/Temp/msohtmlclip1/01/clip_image001.gif[/url]
别的不多说了,直接上干货。
file:///C:/Users/zhang/AppData/Local/Temp/msohtmlclip1/01/clip_image002.gif
[url=]file:///C:/Users/zhang/AppData/Local/Temp/msohtmlclip1/01/clip_image001.gif[/url]
1 /*
2 *file WindowFunction.h
3 *author Vincent Cui
5 *version 0.3
6 *data 31-Oct-2014
7 *brief 各种窗函数的C语言实现
8 */
9
10
11 #ifndef _WINDOWFUNCTION_H_
12 #define _WINDOWFUNCTION_H_
13
14 #include "GeneralConfig.h"
15
16 #define BESSELI_K_LENGTH 10
17
18 #define FLATTOPWIN_A0 0.215578995
19 #define FLATTOPWIN_A1 0.41663158
20 #define FLATTOPWIN_A2 0.277263158
21 #define FLATTOPWIN_A3 0.083578947
22 #define FLATTOPWIN_A4 0.006947368
23
24 #define NUTTALL_A0 0.3635819
25 #define NUTTALL_A1 0.4891775
26 #define NUTTALL_A2 0.1365995
27 #define NUTTALL_A3 0.0106411
28
29 #define BLACKMANHARRIS_A0 0.35875
30 #define BLACKMANHARRIS_A1 0.48829
31 #define BLACKMANHARRIS_A2 0.14128
32 #define BLACKMANHARRIS_A3 0.01168
33
34
35
36 dspErrorStatus taylorWin(dspUint_16 N, dspUint_16 nbar,dspDouble sll, dspDouble **w);
37 dspErrorStatus triangularWin(dspUint_16 N, dspDouble **w);
38 dspErrorStatus tukeyWin(dspUint_16 N, dspDouble r,dspDouble **w);
39 dspErrorStatus bartlettWin(dspUint_16 N, dspDouble **w);
40 dspErrorStatus bartLettHannWin(dspUint_16 N, dspDouble**w);
41 dspErrorStatus blackManWin(dspUint_16 N, dspDouble **w);
42 dspErrorStatus blackManHarrisWin(dspUint_16 N, dspDouble**w);
43 dspErrorStatus bohmanWin(dspUint_16 N, dspDouble **w);
44 dspErrorStatus chebyshevWin(dspUint_16 N, dspDouble r,dspDouble **w);
45 dspErrorStatus flatTopWin(dspUint_16 N, dspDouble **w);
46 dspErrorStatus gaussianWin(dspUint_16 N, dspDouble alpha,dspDouble **w);
47 dspErrorStatus hammingWin(dspUint_16 N, dspDouble **w);
48 dspErrorStatus hannWin(dspUint_16 N, dspDouble **w);
49 dspErrorStatus kaiserWin(dspUint_16 N, dspDouble beta,dspDouble **w);
50 dspErrorStatus nuttalWin(dspUint_16 N, dspDouble **w);
51 dspErrorStatus parzenWin(dspUint_16 N, dspDouble **w);
52 dspErrorStatus rectangularWin(dspUint_16 N, dspDouble**w);
53
54
55
56 #endif
[url=]file:///C:/Users/zhang/AppData/Local/Temp/msohtmlclip1/01/clip_image001.gif[/url]
file:///C:/Users/zhang/AppData/Local/Temp/msohtmlclip1/01/clip_image002.gif
[url=]file:///C:/Users/zhang/AppData/Local/Temp/msohtmlclip1/01/clip_image001.gif[/url]
1 /*
2 *file WindowFunction.c
3 *author Vincent Cui
5 *version 0.3
6 *data 31-Oct-2014
7 *brief 各种窗函数的C语言实现
8 */
9
10
11 #include "WindowFunction.h"
12 #include "GeneralConfig.h"
13 #include "MathReplenish.h"
14 #include "math.h"
15 #include <stdlib.h>
16 #include <stdio.h>
17 #include <string.h>
18
19 /*函数名:taylorWin
20 *说明:计算泰勒窗。泰勒加权函数
21 *输入:
22 *输出:
23 *返回:
24 *调用:prod()连乘函数
25 *其它:用过以后,需要手动释放掉*w的内存空间
26 * 调用示例:ret = taylorWin(99,4, 40, &w); 注意此处的40是正数 表示-40dB
27 */
28 dspErrorStatus taylorWin(dspUint_16N, dspUint_16 nbar, dspDouble sll, dspDouble **w)
29 {
30 dspDouble A;
31 dspDouble *retDspDouble;
32 dspDouble *sf;
33 dspDouble *result;
34 dspDouble alpha,beta,theta;
35 dspUint_16 i,j;
36
37 /*A = R cosh(PI,A) = R*/
38 A =(dspDouble)acosh(pow((dspDouble)10.0,(dspDouble)sll/20.0)) /PI;
39 A = A * A;
40
41 /*开出存放系数的空间*/
42 retDspDouble = (dspDouble*)malloc(sizeof(dspDouble) * (nbar - 1));
43 if(retDspDouble== NULL)
44 returnDSP_ERROR;
45 sf = retDspDouble;
46
47 /*开出存放系数的空间*/
48 retDspDouble = (dspDouble*)malloc(sizeof(dspDouble) * N);
49 if(retDspDouble== NULL)
50 returnDSP_ERROR;
51 result = retDspDouble;
52
53 alpha = prod(1, 1, (nbar - 1));
54 alpha *= alpha;
55 beta = (dspDouble)nbar / sqrt(A + pow((nbar - 0.5), 2) );
56 for(i = 1; i <= (nbar - 1); i++)
57 {
58 *(sf + i - 1) = prod(1,1,(nbar -1 + i)) * prod(1,1,(nbar -1 - i));
59 theta = 1;
60 for(j = 1; j <= (nbar - 1); j++)
61 {
62 theta *= 1 - (dspDouble)(i * i) / ( beta * beta * ( A + (j - 0.5) * (j - 0.5)) );
63 }
64 *(sf + i - 1) = alpha * (dspDouble)theta / (*(sf + i - 1));
65 }
66
67 /*奇数阶*/
68 if((N % 2) == 1)
69 {
70 for(i = 0; i < N; i++)
71 {
72 alpha = 0;
73 for(j = 1; j <= (nbar - 1);j++)
74 {
75 alpha += (*(sf + j- 1)) * cos( 2 * PI * j * (dspDouble)(i - ((N-1)/2))/N );
76 }
77 *(result + i) = 1 + 2 * alpha;
78 }
79 }
80 /*偶数阶*/
81 else
82 {
83 for(i = 0; i < N; i++)
84 {
85 alpha = 0;
86 for(j = 1; j <= (nbar - 1);j++)
87 {
88 alpha += (*(sf + j- 1)) * cos( PI * j * (dspDouble)(2 * (i- (N/2)) + 1) / N );
89 }
90 *(result + i) = 1 + 2 * alpha;
91
92 }
93 }
94 *w = result;
95 free(sf);
96
97 returnDSP_SUCESS;
98
99 }
100
101
102 /*
103 *函数名:triangularWin
104 *说明:计算三角窗函数
105 *输入:
106 *输出:
107 *返回:
108 *调用:
109 *其它:用过以后,需要手动释放掉*w的内存空间
110 * 调用示例:ret = triangularWin(99, &w);
111 */
112 dspErrorStatus triangularWin(dspUint_16 N, dspDouble **w)
113 {
114 dspDouble *ptr;
115 dspUint_16 i;
116
117 ptr = (dspDouble *)malloc(N * sizeof(dspDouble));
118 if(ptr == NULL)
119 return DSP_ERROR;
120
121
122 /*阶数为奇*/
123 if((N % 2) == 1)
124 {
125 for(i = 0; i< ((N - 1)/2); i++)
126 {
127 *(ptr + i) = 2 * (dspDouble)(i + 1) / (N+ 1);
128 }
129 for(i = ((N - 1)/2); i < N; i++)
130 {
131 *(ptr + i) = 2 * (dspDouble)(N - i) / (N + 1);
132 }
133 }
134 /*阶数为偶*/
135 else
136 {
137 for(i = 0; i< (N/2); i++)
138 {
139 *(ptr + i) = (i + i + 1) * (dspDouble)1 / N;
140 }
141 for(i = (N/2); i< N; i++)
142 {
143 *(ptr + i) = *(ptr + N - 1 - i);
144 }
145 }
146 *w= ptr;
147
148 return DSP_SUCESS;
149 }
150
151 /*
152 *函数名:tukeyWin
153 *说明:计算tukey窗函数
154 *输入:
155 *输出:
156 *返回:linSpace()
157 *调用:
158 *其它:用过以后,需要手动释放掉*w的内存空间
159 * 调用示例:ret = tukeyWin(99, 0.5, &w);
160 */
161 dspErrorStatus tukeyWin(dspUint_16 N, dspDouble r,dspDouble **w)
162 {
163 dspErrorStatus retErrorStatus;
164 dspUint_16 index;
165 dspDouble *x,*result,*retPtr;
166 dspDouble alpha;
167
168 retErrorStatus = linSpace(0, 1, N, &x);
169 if(retErrorStatus == DSP_ERROR)
170 return DSP_ERROR;
171
172 result = (dspDouble *)malloc(N * sizeof(dspDouble));
173 if(result == NULL)
174 return DSP_ERROR;
175
176 /*r <= 0 就是矩形窗*/
177 if(r <= 0)
178 {
179 retErrorStatus = rectangularWin(N, &retPtr);
180 if(retErrorStatus == DSP_ERROR)
181 return DSP_ERROR;
182 /*将数据拷出来以后,释放调用的窗函数的空间*/
183 memcpy(result, retPtr, ( N * sizeof(dspDouble)));
184 free(retPtr);
185 }
186 /*r >= 1 就是汉宁窗*/
187 else if(r >= 1)
188 {
189 retErrorStatus = hannWin(N, &retPtr);
190 if(retErrorStatus == DSP_ERROR)
191 return DSP_ERROR;
192 /*将数据拷出来以后,释放调用的窗函数的空间*/
193 memcpy(result, retPtr, ( N * sizeof(dspDouble)));
194 free(retPtr);
195 }
196 else
197 {
198 for(index = 0;index < N; index++)
199 {
200 alpha = *(x + index);
201 if(alpha < (r/2))
202 {
203 *(result + index) =(dspDouble)(1 + cos( 2 * PI* (dspDouble)(alpha - (dspDouble)r/2)/r))/2;
204 }
205 else if((alpha >= (r/2))&& (alpha <(1 - r/2)))
206 {
207 *(result + index) = 1;
208 }
209 else
210 {
211 *(result + index) =(dspDouble)(1 + cos( 2 * PI* (dspDouble)(alpha - 1 + (dspDouble)r/2)/r))/2;
212 }
213
214 }
215 }
216
217 free(x);
218
219 *w= result;
220
221 return DSP_SUCESS;
222
223 }
224
225 /*
226 *函数名:bartlettWin
227 *说明:计算bartlettWin窗函数
228 *输入:
229 *输出:
230 *返回:
231 *调用:
232 *其它:用过以后,需要手动释放掉*w的内存空间
233 * 调用示例:ret = bartlettWin(99, &w);
234 */
235 dspErrorStatus bartlettWin(dspUint_16 N, dspDouble **w)
236 {
237 dspDouble*ret;
238 dspUint_16 n;
239
240 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
241 if(ret == NULL)
242 return DSP_ERROR;
243
244 for(n = 0; n < ( N - 1 ) / 2; n++)
245 {
246 *(ret + n) = 2 * (dspDouble)n / (N - 1);
247 }
248
249 for(n = ( N - 1 ) / 2; n< N; n++)
250 {
251 *(ret + n) = 2 - 2 * (dspDouble)n / ((N - 1 ));
252 }
253
254 *w= ret;
255
256 return DSP_SUCESS;
257 }
258
259 /*
260 *函数名:bartLettHannWin
261 *说明:计算bartLettHannWin窗函数
262 *输入:
263 *输出:
264 *返回:
265 *调用:
266 *其它:用过以后,需要手动释放掉*w的内存空间
267 * 调用示例:ret = bartLettHannWin(99, &w);
268 */
269 dspErrorStatus bartLettHannWin(dspUint_16 N, dspDouble**w)
270 {
271 dspUint_16n;
272 dspDouble *ret;
273
274 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
275 if(ret == NULL)
276 return DSP_ERROR;
277 /*奇*/
278 if(( N % 2 ) == 1)
279 {
280 for(n = 0; n< N; n++)
281 {
282 *(ret + n) = 0.62 - 0.48 * myAbs( ((dspDouble)n / ( N - 1 ) ) - 0.5 ) + 0.38 * cos( 2 * PI * ( ((dspDouble)n / ( N - 1 ) ) - 0.5 ) );
283 }
284 for(n = 0; n< (N-1)/2; n++)
285 {
286 *(ret + n) = *(ret + N - 1 - n);
287 }
288 }
289 /*偶*/
290 else
291 {
292 for(n = 0; n< N; n++)
293 {
294 *(ret + n) = 0.62 - 0.48 * myAbs( ((dspDouble)n / ( N - 1 ) ) - 0.5 ) + 0.38 * cos( 2 * PI * ( ((dspDouble)n / ( N - 1 ) ) - 0.5 ) );
295 }
296 for(n = 0; n< N/2; n++)
297 {
298 *(ret + n) = *(ret + N -1 - n);
299 }
300 }
301
302 *w= ret;
303
304 return DSP_SUCESS;
305
306 }
307
308 /*
309 *函数名:blackManWin
310 *说明:计算blackManWin窗函数
311 *输入:
312 *输出:
313 *返回:
314 *调用:
315 *其它:用过以后,需要手动释放掉*w的内存空间
316 * 调用示例:ret = blackManWin(99, &w);
317 */
318 dspErrorStatus blackManWin(dspUint_16 N, dspDouble **w)
319 {
320 dspUint_16n;
321 dspDouble *ret;
322 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
323 if(ret == NULL)
324 return DSP_ERROR;
325
326 for(n = 0; n < N; n++)
327 {
328 *(ret + n) = 0.42 - 0.5 * cos(2 * PI * (dspDouble)n / ( N - 1 )) + 0.08 * cos( 4 * PI * ( dspDouble )n / ( N - 1 ) );
329 }
330
331 *w= ret;
332
333 return DSP_SUCESS;
334 }
335
336 /*
337 *函数名:blackManHarrisWin
338 *说明:计算blackManHarrisWin窗函数
339 *输入:
340 *输出:
341 *返回:
342 *调用:
343 *其它:用过以后,需要手动释放掉*w的内存空间
344 * 调用示例:ret = blackManHarrisWin(99, &w);
345 * minimum 4-term Blackman-harris window-- From Matlab
346 */
347 dspErrorStatus blackManHarrisWin(dspUint_16 N, dspDouble**w)
348 {
349 dspUint_16n;
350 dspDouble *ret;
351
352 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
353 if(ret == NULL)
354 return DSP_ERROR;
355
356 for(n = 0; n < N; n++)
357 {
358 *(ret + n) = BLACKMANHARRIS_A0 - BLACKMANHARRIS_A1 * cos( 2 * PI * (dspDouble)n/ (N) ) + \
359 BLACKMANHARRIS_A2 * cos( 4 * PI * (dspDouble)n/ (N) ) - \
360 BLACKMANHARRIS_A3 * cos( 6 * PI * (dspDouble)n/ (N) );
361 }
362
363 *w= ret;
364
365 return DSP_SUCESS;
366 }
367
368 /*
369 *函数名:bohmanWin
370 *说明:计算bohmanWin窗函数
371 *输入:
372 *输出:
373 *返回:
374 *调用:
375 *其它:用过以后,需要手动释放掉*w的内存空间
376 * 调用示例:ret = bohmanWin(99, &w);
377 */
378 dspErrorStatus bohmanWin(dspUint_16 N, dspDouble **w)
379 {
380 dspUint_16 n;
381 dspDouble *ret;
382 dspDouble x;
383 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
384 if(ret == NULL)
385 return DSP_ERROR;
386
387 for(n = 0; n < N; n++)
388 {
389 x = -1+ n * (dspDouble)2 / ( N - 1 ) ;
390 /*取绝对值*/
391 x = x >= 0 ? x : ( x * ( -1 ) );
392 *(ret + n) = ( 1 - x ) * cos( PI * x) + (dspDouble)(1 / PI)* sin( PI * x);
393 }
394
395 *w= ret;
396
397 return DSP_SUCESS;
398 }
399
400 /*
401 *函数名:chebyshevWin
402 *说明:计算chebyshevWin窗函数
403 *输入:
404 *输出:
405 *返回:
406 *调用:
407 *其它:用过以后,需要手动释放掉*w的内存空间
408 * 调用示例:ret = chebyshevWin(99,100, &w);
409 */
410 dspErrorStatus chebyshevWin(dspUint_16 N, dspDouble r,dspDouble **w)
411 {
412 dspUint_16 n,index;
413 dspDouble *ret;
414 dspDouble x, alpha, beta, theta, gama;
415
416 ret= (dspDouble *)malloc(N * sizeof(dspDouble));
417 if(ret == NULL)
418 return DSP_ERROR;
419
420
421 /*10^(r/20)*/
422 theta = pow((dspDouble)10, (dspDouble)(myAbs(r)/20));
423 beta = pow(cosh(acosh(theta)/(N - 1)),2);
424 alpha= 1 - (dspDouble)1 / beta;
425
426 if((N % 2) == 1)
427 {
428 /*计算一半的区间*/
429 for( n = 1; n< ( N + 1 ) / 2; n++)
430 {
431 gama = 1;
432 for(index = 1;index < n; index++)
433 {
434 x = index * (dspDouble)( N - 1 - 2 * n + index) /(( n - index ) * (n + 1 -index));
435 gama = gama * alpha * x + 1;
436 }
437 *(ret + n) = (N - 1) * alpha * gama;
438 }
439
440 theta = *( ret + (N - 1)/2 );
441 *ret = 1;
442
443 for(n = 0; n< ( N + 1 ) / 2; n++)
444 {
445 *(ret + n) = (dspDouble)(*(ret + n)) / theta;
446 }
447
448 /*填充另一半*/
449 for(; n < N; n++)
450 {
451 *(ret + n) = ret[N - n - 1];
452 }
453 }
454 else
455 {
456 /*计算一半的区间*/
457 for( n = 1; n< ( N + 1 ) / 2; n++)
458 {
459 gama = 1;
460 for(index = 1;index < n; index++)
461 {
462 x = index * (dspDouble)( N - 1 - 2 * n + index) /(( n - index ) * (n + 1 -index));
463 gama = gama * alpha * x + 1;
464 }
465 *(ret + n) = (N - 1) * alpha * gama;
466 }
467
468 theta = *( ret + (N/2) - 1);
469 *ret = 1;
470
471 for(n = 0; n< ( N + 1 ) / 2; n++)
472 {
473 *(ret + n) = (dspDouble)(*(ret + n)) / theta;
474 }
475
476 /*填充另一半*/
477 for(; n < N; n++)
478 {
479 *(ret + n) = ret[N - n - 1];
480 }
481 }
482
483
484 *w= ret;
485
486 return DSP_SUCESS;
487 }
488
489 /*
490 *函数名:flatTopWin
491 *说明:计算flatTopWin窗函数
492 *输入:
493 *输出:
494 *返回:
495 *调用:
496 *其它:用过以后,需要手动释放掉*w的内存空间
497 * 调用示例:ret = flatTopWin(99, &w);
498 */
499 dspErrorStatus flatTopWin(dspUint_16 N, dspDouble **w)
500 {
501 dspUint_16 n;
502 dspDouble *ret;
503 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
504 if(ret == NULL)
505 return DSP_ERROR;
506
507 for(n = 0; n < N; n++)
508 {
509 *(ret + n) = FLATTOPWIN_A0 - FLATTOPWIN_A1 * cos(2 * PI * (dspDouble)n / (N - 1)) +\
510 FLATTOPWIN_A2 * cos(4 * PI * (dspDouble)n / (N - 1)) -\
511 FLATTOPWIN_A3 * cos(6 * PI * (dspDouble)n / (N - 1)) +\
512 FLATTOPWIN_A4 * cos(8 * PI * (dspDouble)n / (N - 1));
513 }
514
515 *w= ret;
516
517 return DSP_SUCESS;
518 }
519
520
521 /*
522 *函数名:gaussianWin
523 *说明:计算gaussianWin窗函数
524 *输入:
525 *输出:
526 *返回:
527 *调用:
528 *其它:用过以后,需要手动释放掉*w的内存空间
529 * 调用示例:ret = gaussianWin(99,2.5, &w);
530 */
531 dspErrorStatus gaussianWin(dspUint_16 N, dspDouble alpha,dspDouble **w)
532 {
533 dspUint_16 n;
534 dspDouble k, beta, theta;
535 dspDouble *ret;
536
537 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
538 if(ret == NULL)
539 return DSP_ERROR;
540
541 for(n =0; n < N; n++)
542 {
543 if((N % 2) == 1)
544 {
545 k = n - (N - 1)/2;
546 beta = 2 * alpha * (dspDouble)k / (N - 1);
547 }
548 else
549 {
550 k = n - (N)/2;
551 beta = 2 * alpha * (dspDouble)k / (N - 1);
552 }
553
554 theta = pow(beta, 2);
555 *(ret + n) = exp((-1) * (dspDouble)theta / 2);
556 }
557
558 *w= ret;
559
560 return DSP_SUCESS;
561 }
562
563 /*
564 *函数名:hammingWin
565 *说明:计算hammingWin窗函数
566 *输入:
567 *输出:
568 *返回:
569 *调用:
570 *其它:用过以后,需要手动释放掉*w的内存空间
571 * 调用示例:ret = hammingWin(99, &w);
572 */
573 dspErrorStatus hammingWin(dspUint_16 N, dspDouble **w)
574 {
575 dspUint_16 n;
576 dspDouble *ret;
577 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
578 if(ret == NULL)
579 return DSP_ERROR;
580
581 for(n = 0; n < N; n++)
582 {
583 *(ret + n) = 0.54 - 0.46 * cos (2 * PI * ( dspDouble )n / ( N - 1 ) );
584 }
585
586 *w= ret;
587
588 return DSP_SUCESS;
589 }
590
591 /*
592 *函数名:hannWin
593 *说明:计算hannWin窗函数
594 *输入:
595 *输出:
596 *返回:
597 *调用:
598 *其它:用过以后,需要手动释放掉*w的内存空间
599 * 调用示例:ret = hannWin(99, &w);
600 */
601 dspErrorStatus hannWin(dspUint_16 N, dspDouble **w)
602 {
603 dspUint_16 n;
604 dspDouble *ret;
605 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
606 if(ret == NULL)
607 return DSP_ERROR;
608
609 for(n = 0; n < N; n++)
610 {
611 *(ret + n) = 0.5 * ( 1 -cos( 2 * PI * (dspDouble)n / (N - 1)));
612 }
613
614 *w= ret;
615
616 return DSP_SUCESS;
617 }
618
619 /*
620 *函数名:kaiserWin
621 *说明:计算kaiserWin窗函数
622 *输入:
623 *输出:
624 *返回:
625 *调用:besseli()第一类修正贝塞尔函数
626 *其它:用过以后,需要手动释放掉*w的内存空间
627 * 调用示例:ret = kaiserWin(99, 5, &w);
628 */
629 dspErrorStatus kaiserWin(dspUint_16 N, dspDouble beta,dspDouble **w)
630 {
631 dspUint_16n;
632 dspDouble *ret;
633 dspDouble theta;
634
635 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
636 if(ret == NULL)
637 return DSP_ERROR;
638
639 for(n = 0; n < N; n++)
640 {
641 theta = beta * sqrt( 1 - pow( ( (2 *(dspDouble)n/(N -1)) - 1),2 ) );
642 *(ret + n) = (dspDouble)besseli(0,theta, BESSELI_K_LENGTH) / besseli(0, beta,BESSELI_K_LENGTH);
643 }
644
645 *w= ret;
646
647 return DSP_SUCESS;
648 }
649
650 /*
651 *函数名:nuttalWin
652 *说明:计算nuttalWin窗函数
653 *输入:
654 *输出:
655 *返回:
656 *调用:
657 *其它:用过以后,需要手动释放掉*w的内存空间
658 * 调用示例:ret = nuttalWin(99, &w);
659 */
660 dspErrorStatus nuttalWin(dspUint_16 N, dspDouble **w)
661 {
662 dspUint_16 n;
663 dspDouble *ret;
664
665 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
666 if(ret == NULL)
667 return DSP_ERROR;
668
669 for(n = 0; n < N; n++)
670 {
671 *(ret + n) =NUTTALL_A0 - NUTTALL_A1 * cos(2 * PI* (dspDouble)n / (N - 1)) +\
672 NUTTALL_A2 *cos(4 * PI * (dspDouble)n / (N - 1)) -\
673 NUTTALL_A3 *cos(6 * PI * (dspDouble)n / (N - 1));
674
675 }
676
677 *w = ret;
678
679 return DSP_SUCESS;
680 }
681
682 /*
683 *函数名:parzenWin
684 *说明:计算parzenWin窗函数
685 *输入:
686 *输出:
687 *返回:
688 *调用:
689 *其它:用过以后,需要手动释放掉*w的内存空间
690 * 调用示例:ret = parzenWin(99, &w);
691 */
692 dspErrorStatus parzenWin(dspUint_16 N, dspDouble **w)
693 {
694 dspUint_16 n;
695 dspDouble *ret;
696 dspDouble alpha,k;
697
698 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
699 if(ret == NULL)
700 return DSP_ERROR;
701
702 if(( N % 2) == 1)
703 {
704 for(n = 0; n< N; n++)
705 {
706 k = n - (N - 1) / 2;
707 alpha = 2 * (dspDouble)myAbs(k) / N;
708 if(myAbs(k) <= (N - 1) / 4)
709 {
710 *(ret + n) = 1 - 6 * pow(alpha,2) + 6 * pow(alpha, 3);
711 }
712 else
713 {
714 *(ret + n) = 2 * pow( (1 - alpha), 3 );
715 }
716
717 }
718 }
719 else
720 {
721 for(n = 0; n< N; n++)
722 {
723 k = n - (N - 1) / 2;
724 alpha = 2 * (dspDouble)myAbs(k) / N;
725 if(myAbs(k) <= (dspDouble)(N -1) / 4)
726 {
727 *(ret + n) = 1 - 6 * pow(alpha,2) + 6 * pow(alpha, 3);
728 }
729 else
730 {
731 *(ret + n) = 2 * pow( (1 - alpha), 3 );
732 }
733
734 }
735 }
736
737
738
739 *w= ret;
740
741 return DSP_SUCESS;
742 }
743
744 /*
745 *函数名:rectangularWin
746 *说明:计算rectangularWin窗函数
747 *输入:
748 *输出:
749 *返回:
750 *调用:
751 *其它:用过以后,需要手动释放掉*w的内存空间
752 * 调用示例:ret = rectangularWin(99, &w);
753 */
754 dspErrorStatus rectangularWin(dspUint_16 N, dspDouble **w)
755 {
756 dspUint_16 n;
757 dspDouble *ret;
758
759 ret = (dspDouble *)malloc(N * sizeof(dspDouble));
760 if(ret == NULL)
761 return DSP_ERROR;
762
763 for(n = 0; n < N; n++)
764 {
765 *(ret + n) = 1;
766 }
767
768 *w= ret;
769
770 return DSP_SUCESS;
771 }
[url=]file:///C:/Users/zhang/AppData/Local/Temp/msohtmlclip1/01/clip_image001.gif[/url]