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

    欢迎多交流!

  • 相关阅读:
    浅谈ruby中的block及yield
    Node.js使用TCP通讯
    JavaScript垃圾回收机制
    node.js的Promise库-bluebird示例
    使用pkg打包Node.js应用的方法步骤
    windows server 2012 安装 VC14(VC2015) 安装失败解决方案
    pm2命令,端口查询,mongodb服务启动,nginx服务启动,n模块的使用,搭建nodejs服务器环境,搭建oracledb服务器环境 linux的环境搭建
    linux+node.js+redis+mongodb+nginx环境的搭建
    nginx.exe启动错误:bind() to 0.0.0.0:80 failed (10013: An attempt was made to access a socket in a way forbidden by its access permissions)
    windows下nginx的安装及使用
  • 原文地址:https://www.cnblogs.com/mildsim/p/4067374.html
Copyright © 2020-2023  润新知