窗函数的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