窗函数的C语言实现

一般的讲数字信号处理的书中都会提到窗函数。大多数只会提及其中的几种。这里我把这些窗都用C语言实现了一下，都不复杂，但如果要自己去弄也挺费时间。所有函数都用Matlab验证了。包括以下窗：

 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;

别的不多说了，直接上干货。

 1 /*
 2  *file               WindowFunction.h
 3  *author          Vincent Cui
 4  *e-mail           whcui1987@163.com
 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

WindowFunction.h

  1 /*
  2  *file               WindowFunction.c
  3  *author          Vincent Cui
  4  *e-mail            whcui1987@163.com
  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_16 N, 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         return DSP_ERROR;
 45     sf = retDspDouble;
 46 
 47     /*开出存放系数的空间*/
 48     retDspDouble = (dspDouble *)malloc(sizeof(dspDouble) * N);
 49     if(retDspDouble == NULL)
 50         return DSP_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     return DSP_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_16 n;
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_16 n;
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_16 n;
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_16 n;
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 }

WindowFunction.c

欢迎多交流！

posted @ 2014-11-01 16:51 Vincent.Cui 阅读(9056) 评论(1) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

窗函数的C语言实现

公告