拙网论坛

 找回密码
 立即注册
搜索
热搜: 活动 交友 discuz
查看: 232|回复: 0

窗函数的C语言实现

[复制链接]

949

主题

1001

帖子

3736

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3736
发表于 2018-8-2 15:43:11 | 显示全部楼层 |阅读模式
https://www.cnblogs.com/mildsim/p/4067374.html

一般的讲数字信号处理的书中都会提到窗函数。大多数只会提及其中的几种。这里我把这些窗都用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
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
[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
  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_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]


回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Archiver|手机版|小黑屋|抱朴守拙BBS

GMT+8, 2025-6-13 02:41 , Processed in 0.237731 second(s), 18 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表