窗函数的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           [email protected]
 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            [email protected]
  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

欢迎多交流!

时间: 2024-10-09 08:23:51

窗函数的C语言实现的相关文章

使用R语言计算均值,方差等

R语言对于数值计算很方便,最近用到了计算方差,标准差的功能,特记录. 数据准备 height <- c(6.00, 5.92, 5.58, 5.92) 1 计算均值 mean(height) [1] 5.855 2 计算中位数 median(height) [1] 5.92 3 计算标准差 sd(height) [1] 0.1871719 4 计算方差 var(height) [1] 0.03503333 5 计算两个变量之间的相关系数 cor(height,log(height)) [1] 0

GCC在C语言中内嵌汇编 asm __volatile__ 【转】

转自:http://blog.csdn.net/pbymw8iwm/article/details/8227839 在内嵌汇编中,可以将C语言表达式指定为汇编指令的操作数,而且不用去管如何将C语言表达式的值读入哪个寄存器,以及如何将计算结果写回C 变量,你只要告诉程序中C语言表达式与汇编指令操作数之间的对应关系即可, GCC会自动插入代码完成必要的操作. 1.简单的内嵌汇编 例: __asm__ __volatile__("hlt"); "__asm__"表示后面的

C语言轻松高效学习方法之:多种方法实现

多种方法实现同一个功能,可以调动你学的所有知识去做,有助于你学的融会贯通. 下面举例来看: 实现功能:求一个整数的位数: 实现语言:C语言: 开发环境:Visual Studio 2017 如:3215是4位数 实现原理: 3215/10 = 321 ----1位数 321/10 = 32 ----又是1位数 32/10 = 3 ----又是1位数 3/10 = 0 ----又是1位数 共4位数,且终止计算条件是/10结果为0的时候: 根据这个原理,先写一个最笨的原始方法: 效果: 这种实现方案

轻松学习C语言编程的秘诀:总结+灵感

目前在准备一套C语言的学习教程,所以我这里就以C语言编程的学习来讲.注意,讲的是"轻松学习",那种不注重方法,拼命玩命的方式也有其效果,但不是我提倡的.我讲究的是在方式方法对头.适合你.减轻你学习负担和心里压力的前提下,才适当的抓紧时间. 因此,探索一种很好的学习方法就是我所研究的主要内容. 众所周知,学习C语言并非易事,要学好它更是难上加难.这和你期末考试背会几个题目的答案考上满分没多大关系,也就是说你考试满分也说明不了你学好.学精通了C语言.那么怎么才算学精通C语言?闭着眼睛对自己

详解go语言的array和slice 【二】

上一篇  详解go语言的array和slice [一]已经讲解过,array和slice的一些基本用法,使用array和slice时需要注意的地方,特别是slice需要注意的地方比较多.上一篇的最后讲解到创建新的slice时使用第三个索引来限制slice的容量,在操作新slice时,如果新slice的容量大于长度时,添加新元素依然后使源的相应元素改变.这一篇里我会讲解到如何避免这些问题,以及迭代.和做为方法参数方面的知识点. slice的长度和容量设置为同一个值 如果在创建新的slice时我们把

自动生成小学四则运算题目(C语言)

这个简易四则运算是我在百度上找的博主叫53余雅诗的一篇c语言代码,网址为http://www.cnblogs.com/ys1101/p/4368103.html,功能是可以选择加减乘除进行简易的四则运算,判断对错.我在VS2017上编译没有bug,因为功能只有整数运算,所以我在此基础上加了真分数的四则运算以及统计得分等,最后成功运行程序.我把我的源代码放在github上,地址为https://github.com/xiaofancheng/helloworld.

PAT 1009 说反话 C语言

给定一句英语,要求你编写程序,将句中所有单词的顺序颠倒输出. 输入格式:测试输入包含一个测试用例,在一行内给出总长度不超过80的字符串.字符串由若干单词和若干空格组成,其中单词是由英文字母(大小写有区分)组成的字符串,单词之间用1个空格分开,输入保证句子末尾没有多余的空格. 输出格式:每个测试用例的输出占一行,输出倒序后的句子. 输入样例: Hello World Here I Come 输出样例: Come I Here World Hello 1 #include<stdio.h> 2 #

PAT 1006 换个格式输出 C语言

让我们用字母B来表示"百".字母S表示"十",用"12...n"来表示个位数字n(<10),换个格式来输出任一个不超过3位的正整数.例如234应该被输出为BBSSS1234,因为它有2个"百".3个"十".以及个位的4. 输入格式:每个测试输入包含1个测试用例,给出正整数n(<1000). 输出格式:每个测试用例的输出占一行,用规定的格式输出n. 输入样例1: 234 输出样例1: BBSSS1

Go语言 IDE之Gogland配置使用

Gogland 是 JetBrains 公司推出的 Go 语言集成开发环境.Gogland 同样基于 IntelliJ 平台开发,支持 JetBrains 的插件体系.目前正式版尚未发布.官方:https://www.jetbrains.com/go/.关于使用,即将开始咯! 一.安装Golang 1) 首先到https://golang.org/dl/选择适合你系统的安装包,(墙内:http://golangtc.com/download). 2)下载完成安装到指定目录即可.我这里是(D:\G