C语言完成窗口算法

1、C 语言完成窗口算法/* 窗类型 */typedef enumBartlett = 0,BartLettHann,BlackMan,BlackManHarris,Bohman,Chebyshev,FlatTop,Gaussian,Hamming,Hann,Kaiser,Nuttal,Parzen,Rectangular,Taylor,Triangular,TukeywinType;/* *fileWindowFunction.h *authorVincent Cui *e-mailwhcui1987*version0.3 *data31-Oct-2014 *brief各种窗函数的 C 语言实现

2、 */#ifndefWINDOWFUNCTION_H_#defineWINDOWFUNCTION_H #include "GeneralConfig.h"#define BESSELI_K_LENGTH10#define FLATTOPWIN_A00.215578995#defineFLATTOPWIN_A10.41663158#defineFLATTOPWIN_A20.277263158#defineFLATTOPWIN_A30.083578947#defineFLATTOPWIN_A40.006947368#define NUTTALL_A00.3635819#defi

3、neNUTTALL_A10.4891775#defineNUTTALL_A20.1365995#define0.35875#define0.0106411NUTTALL_A3 #define BLACKMANHARRIS_A0BLACKMANHARRIS_A10.48829#defineBLACKMANHARRIS_A20.14128#defineBLACKMANHARRIS_A30.01168dspErrorStatustaylorWin(dspUint_16 N, dspUint_16 nbar, dspDouble sll,triangularWin(dspUint_16tukeyWin

4、(dspUint_16 N,dspDouble *w);dspErrorStatusN, dspDouble *w);dspErrorStatus dspDouble r, dspDouble *w);dspErrorStatus bartlettWin(dspUint_16 N, dspDouble *w);dspErrorStatus bartLettHannWin(dspUint_16 N, dspDouble *w);dspErrorStatus blackManWin(dspUint_16 N, dspDouble *w);dspErrorStatus blackManHarrisW

5、in(dspUint_16 N, dspDouble*w);dspErrorStatus bohmanWin(dspUint_16 N,dspDouble *w);dspErrorStatuschebyshevWin(dspUint_16N, dspDouble r, dspDouble *w);dspErrorStatus flatTopWin(dspUint_16 N, dspDouble *w);dspErrorStatus gaussianWin(dspUint_16 N, dspDouble alpha, dspDouble *w);dspErrorStatus hammingWin

6、(dspUint_16 N,hannWin(dspUint_16 N,kaiserWin(dspUint_16 N,dspDouble *w);dspErrorStatus dspDouble *w);dspErrorStatus dspDouble beta, dspDouble *w);dspErrorStatus nuttalWin(dspUint_16 N, dspDouble *w);dspErrorStatus parzenWin(dspUint_16 N, dspDouble *w);dspErrorStatus rectangularWin(dspUint_16 N, dspD

7、ouble *w);#endifWindowFunction.h/* *fileWindowFunction.c *authorVincent Cui *e-mailwhcui1987 *version0.3 *data31-Oct-2014 *brief各种窗函数的 C 语言实现*/#include "WindowFunction.h"#include "GeneralConfig.h"#include "MathReplenish.h"#include "math.h"#include <stdl

8、ib.h>#include <stdio.h>#include <string.h>/*函数名： taylorWin * 说明：计算泰勒窗。泰勒加权函数输入：*输出：*返回：*调用：prod()连乘函数*其它:用过以后，需要手动释放掉 *w 的内存空间 *调用示例： ret = taylorWin(99, 4, 40, &w); 注意此处的 40 是 正数 表示 -40dB */dspErrorStatus taylorWin(dspUint_16 N, dspUint_16 nbar, dspDouble

9、sll, dspDouble*w)dspDouble A; dspDouble *retDspDouble;dspDouble *sf; dspDouble *result; dspDouble alpha,beta,theta; dspUint_16 i,j;/*A = R cosh(PI, A) = R*/ A =(dspDouble)acosh(pow(dspDouble)10.0,(dspDouble)sll/20.0) /PI; A = A * A;/* 开出存放系数的空间 */retDspDouble = (dspDouble *)malloc(sizeof(dspDouble)

10、*(nbar - 1);if(retDspDouble = NULL)returnDSP_ERROR; sf = retDspDouble;/* 开出存放系数的空间 */ retDspDouble = (dspDouble*)malloc(sizeof(dspDouble) * N);if(retDspDouble =NULL)return DSP_ERROR; result =retDspDouble;alpha = prod(1, 1, (nbar - 1); alpha *= alpha;beta= (dspDouble)nbar / sqrt( A + pow(nbar - 0.5),

11、 2) );for(i =1; i <= (nbar - 1); i+)*(sf + i - 1) =prod(1,1,(nbar -1 + i) * prod(1,1,(nbar -1 - i);theta =1;for(j = 1; j <= (nbar - 1); j+)theta *= 1 - (dspDouble)(i * i) / ( beta * beta *( A + (j - 0.5) * (j - 0.5) );*(sf + i - 1) =alpha * (dspDouble)theta / (*(sf + i - 1);/*奇数阶 */ if

12、(N % 2) = 1)for(i = 0; i< N; i+)alpha = 0;for(j = 1; j <= (nbar - 1); j+)alpha += (输出： *返回： * 调用：(sf + j - 1) * cos( 2 * PI * j *(dspDouble)(i - (N-1)/2)/N );*(result + i) = 1 + 2 * alpha; /*偶数阶*/else for(i = 0; i < N; i+)alpha = 0;for(j = 1; j <= (nbar- 1); j+)alpha += (

13、*(sf + j -1) * cos( PI * j * (dspDouble)(2 * (i - (N/2) + 1) /N );*(result + i) = 1 + 2 *alpha; *w = result;free(sf);return DSP_SUCESS;/* * 函数名： triangularWin * 说明：计算三角窗函数* 输入：*其它：用过以后，需要手动释放掉 *w 的内存空间 *调用示例： ret =*/dspErrorStatustriangularWin(99, &w);dspDoubletriangularWin(dspUint_16 N, dsp

14、Double *w) *ptr; dspUint_16 i;ptr = (dspDouble *)malloc(N * sizeof(dspDouble);if(ptr = NULL)return DSP_ERROR; /* 阶数为奇 */if(N % 2) = 1)for(i = 0; i < (N -1)/2); i+)*(ptr + i) = 2 *(dspDouble)(i + 1) / (N + 1);for(i = (N -*(ptr + i) = 2 * /* 阶数为偶1)/2); i < N; i+) (dspDouble)(N - i) / (N

15、+ 1);if(result = NULL)return DSP_ERROR;*/elsefor(i = 0; i < (N/2); i+)N;*(ptr + i) = (i + i + 1) * (dspDouble)1 /for(i = (N/2); i < N; i+)*(ptr + i) = *(ptr + N - 1 - i);*w = ptr;return DSP_SUCESS;tukey 窗函数 * 输入： */* * 函数名： tukeyWin * 说明：计算输出： *返回： linSpace() *调用：*其它：用过以后，需要手动释放掉 *w 的内

16、存空间 *调用示例： ret =index; dspDouble*x,*result,*retPtr;dspDoublealpha;tukeyWin(99, 0.5, &w); */dspErrorStatus tukeyWin(dspUint_16 N, dspDouble r, dspDouble*w)dspErrorStatus retErrorStatus; dspUint_16returnretErrorStatus = linSpace(0, 1, N, &x);if(retErrorStatus = DSP_ERROR)DSP_ERROR;resu

17、lt = (dspDouble *)malloc(N * sizeof(dspDouble);/*r <= 0 就是矩形窗 */ if(r <= 0)returnretErrorStatus = rectangularWin(N, &retPtr);if(retErrorStatus = DSP_ERROR)DSP_ERROR;/*将数据拷出来以后，释放调用的窗函数的空间 */memcpy(result, retPtr, ( N *sizeof(dspDouble);free(retPtr); /*r >=1 就是汉宁窗 */else

18、if(r >= 1)(dspDouble)r/2)/r)/2;returnretErrorStatus = hannWin(N, &retPtr);if(retErrorStatus = DSP_ERROR)DSP_ERROR;/*将数据拷出来以后，释放调用的窗函数的空间 */memcpy(result, retPtr, ( N *sizeof(dspDouble);free(retPtr); elsefor(index = 0; index < N; index+)alpha = *(x + index);if(alpha <(r/2

19、)*(result + index) =(dspDouble)(1 + cos( 2 * PI * (dspDouble)(alpha -(dspDouble)r/2)/r)/2;elseif(alpha >= (r/2) && (alpha <(1 - r/2)*(result + index) = 1;else*(result + index) =(dspDouble)(1 + cos( 2 * PI * (dspDouble)(alpha - 1 +free(x);*w = result;return DSP_SUCESS;/*

20、 * 函数名：bartlettWin * 说明：计算 bartlettWin 窗函数 * 输入： * 输出：* 返回： * 调用：*其它：用过以后，需要手动释放掉 *w的内存空间 *调用示例： ret =bartlettWin(99, &w); */dspErrorStatusdspDoublebartlettWin(dspUint_16 N, dspDouble *w) *ret; dspUint_16 n;ret = (dspDouble *)malloc(N * sizeof(dspDouble);if(ret = NULL)return DSP_ERROR;*(ret

21、+ n)for(n = 0; n < ( N - 1 ) / 2; n+)= 2 * (dspDouble)n / (N - 1);*(ret + n)for(n = ( N - 1 ) / 2; n < N; n+)= 2 - 2 * (dspDouble)n / ( N - 1 );*w = ret;return DSP_SUCESS;/* * 函数名： bartLettHannWin * 说明：计算bartLettHannWin*其它：用过以窗函数 * 输入： * 输出： * 返回： * 调用：后，需要手动释放掉 *w 的内存空间 *调用示例： ret= ba

22、rtLettHannWin(99, &w); */dspErrorStatus bartLettHannWin(dspUint_16 N, dspDouble*w)dspUint_16 n; dspDouble *ret;ret = (dspDouble *)malloc(N * sizeof(dspDouble);if(ret = NULL)return DSP_ERROR;/* 奇*/if( N % 2 ) = 1)for(n = 0; n < N; n+)*(ret + n) = 0.62 - 0.48 *myAbs( ( (dspDouble)n / (

23、N - 1 ) ) - 0.5 ) + 0.38 * cos( 2 * PI *( (dspDouble)n / ( N - 1 ) ) - 0.5 ) );for(n =0; n < (N-1)/2; n+)*(ret + n) =*(ret + N - 1 - n);/* 偶*/elsefor(n = 0; n < N; n+)*(ret + n) = 0.62 - 0.48 *myAbs( ( (dspDouble)n / ( N - 1 ) ) - 0.5 ) + 0.38 * cos( 2 * PI *for(n =( (dspDouble)n / ( N

24、 - 1 ) ) - 0.5 ) );0; n < N/2; n+)*(ret + n) = *(ret +N -1 - n);*w = ret;return DSP_SUCESS;/* * 函数名： blackManWin * 说明：计算 blackManWin 窗函数*输入： *输出： * 返回： *调用： *其它：用过以后，需要手动释放掉 *w 的内存空间 *调用示例： ret =blackManWin(99, &w); */dspErrorStatusblackManWin(dspUint_16 N, dspDouble *w)dspUint_16n; d

25、spDouble *ret; ret = (dspDouble *)malloc(N *sizeof(dspDouble); if(ret = NULL)returnDSP_ERROR;for(n = 0; n < N; n+)*(ret + n) = 0.42 -0.5 * cos(2 * PI * (dspDouble)n / ( N - 1 ) + 0.08 * cos( 4 * PI* ( dspDouble )n / ( N - 1 ) );*w = ret;return DSP_SUCESS;/* * 函数名： blackManHarrisWin * 说明：计算bla

26、ckManHarrisWin 窗函数 * 输入： *输出： * 返回：* 调用：*其它：用过以后，需要手动释放掉*w 的内存空间dspDouble *w)if(ret= NULL)return DSP_ERROR;*(ret + n) =for(n = 0; n < N; n+)BLACKMANHARRIS_A0 - BLACKMANHARRIS_A1 调用示例： ret = blackManHarrisWin(99, &w);minimum 4-term Blackman-harris window - From Matlab */dspErrorStatus

27、blackManHarrisWin(dspUint_16 N,dspUint_16 n; dspDouble *ret;cos( 2 * PI * (dspDouble)n / (N) ) + BLACKMANHARRIS_A2 * cos(4 * PI * (dspDouble)n/(N) ) - BLACKMANHARRIS_A3 * cos(6 * PI * (dspDouble)n/(N) ); *w = ret;return DSP_SUCESS;/* * 函数名： bohmanWin *说明：计算 bohmanWin 窗函数 *输入： * 输出： * 返回：* 调用： * 其它：用

28、过以后，需要手动释放掉 *w 的内存空间调用示例： ret =bohmanWin(99, &w);*/dspErrorStatusbohmanWin(dspUint_16 N, dspDouble *w)dspUint_16 n;dspDouble *ret; dspDouble x;ret = (dspDouble*)malloc(N * sizeof(dspDouble);if(ret = NULL)return DSP_ERROR;for(n = 0; n < N; n+)-1 +(dspDouble)2 / ( N - 1 ) ;/* 取绝对值 */&am

29、p;gt;= 0 ? x : ( x * ( -1 ) );*(ret + n) =cos( PI * x) + (dspDouble)(1 / PI) * sin( PI * x);*w = ret;return DSP_SUCESS;/* * 函数名： chebyshevWin * 说明：计算 chebyshevWin 窗函数 * 输入： * 输出： * 返回： * 调用： * 其它：用过以后，需要手动释放掉 *w 的内存空间 *调用示例： ret =chebyshevWin(99,100, &w); */dspErrorStatus chebyshevWin(dspUin

30、t_16 N, dspDouble r, dspDouble*w)dspUint_16 n,index; dspDouble *ret;dspDouble x, alpha, beta, theta, gama;ret = (dspDouble *)malloc(N * sizeof(dspDouble);if(ret = NULL)return DSP_ERROR;/*10八(/20)*/theta = pow(dspDouble)10,(dspDouble)(myAbs(r)/20);beta =pow(cosh(acosh(theta)/(N - 1),2);alpha = 1 - (d

31、spDouble)1/ beta;if(N % 2) = 1)/*计算一半的区间 */for( n = 1; n < ( N + 1 ) / 2; n+ )< n; index+)gama = 1;for(index = 1; indexx = index *(dspDouble)( N - 1 - 2 * n + index) /( n - index ) * (n + 1-index);gama = gama * alpha * x +1;*(ret + n) = (N - 1) * alpha *gama;theta = *( ret + (N - 1)/2

32、);*ret = 1;for(n = 0; n < ( N + 1 ) / 2; n+ )*(ret + n) = (dspDouble)(*(ret + n) /theta;/* 填充另一半 */for(; n < N; n+)*(ret + n) = retN - n - 1;else/* 计算一半的区间 */for( n = 1; n< ( N+ 1 ) / 2; n+ )gama = 1;for(index = 1; index < n; index+)x = index * (dspDouble)( N - 1 - 2 * n

33、+index) /( n - index ) * (n + 1 -index);gama= gama * alpha * x + 1;*(ret + n)= (N - 1) * alpha * gama;theta = *( ret + (N/2) - 1);*ret = 1;theta;for(n = 0; n < ( N + 1 ) / 2; n+ )*(ret + n) = (dspDouble)(*(ret + n) /* 填充另一半 */for(; n < N; n+)*(ret + n) = retN - n - 1;*w = ret;return DS

34、P_SUCESS;/* * 函数名： flatTopWin * 说明：计算 flatTopWin 窗函数 * 输入： * 输出： * 返回： * 调用： * 其它：用过以后，需要手动释放掉 *w 的内存空间 *调用示例： ret =flatTopWin(99, &w); */dspErrorStatusflatTopWin(dspUint_16 N, dspDouble *w)dspUint_16 n;dspDouble *ret; ret = (dspDouble *)malloc(N *sizeof(dspDouble); if(ret = NULL)returnDSP_E

35、RROR;for(n = 0; n < N; n+)*(ret + n) =FLATTOPWIN_A0 - FLATTOPWIN_A1 * cos(2 * PI *(dspDouble)n / (N - 1) +FLATTOPWIN_A2 * cos(4 * PI * (dspDouble)n / (N - 1) -FLATTOPWIN_A3 * cos(6 * PI * (dspDouble)n / (N - 1) +FLATTOPWIN_A4 * cos(8 * PI * (dspDouble)n / (N -1); *w = ret;return DSP_SUCESS;/*

36、 * 函数名： gaussianWin *说明：计算 gaussianWin 窗函数 * 输入： *输出： *返回： * 调用：*其它：用过以后，需要手动释放掉*w 的内存空间 *调用示例： ret = gaussianWin(99,2.5, &w);*/dspErrorStatusgaussianWin(dspUint_16 N, dspDoublealpha, dspDouble *w)dspUint_16 n; dspDouble k,beta, theta; dspDouble *ret;ret = (dspDouble *)malloc(N * sizeof(dspD

37、ouble);if(ret = NULL)return DSP_ERROR;for(n=0; n < N; n+) if(N % 2) = 1)k = n - (N - 1)/2;beta = 2 * alphaelse* (dspDouble)k / (N - 1);k = n - (N)/2;beta = 2 * alpha *theta =(dspDouble)k / (N - 1);pow(beta, 2);*(ret + n) = exp(-1) * (dspDouble)theta /2); *w = ret;return DSP_SUCESS;/* * 函数名：ha

38、mmingWin *说明：计算 hammingWin 窗函数 *输入：*输出：*返回：* 调用： * 其它：用过以后，需要手动释放掉 *w 的内存空间调用示例： ret =hammingWin(99, &w);*/dspErrorStatushammingWin(dspUint_16 N, dspDouble *w)dspUint_16n; dspDouble *ret; ret = (dspDouble *)malloc(N *sizeof(dspDouble); if(ret = NULL)returnDSP_ERROR;for(n = 0; n < N; n

39、+)*(ret + n) = 0.54 -0.46 * cos (2 * PI *( dspDouble )n / ( N - 1 ) );*w = ret;return DSP_SUCESS;/* * 函数名： hannWin * 说明：计算 hannWin 窗函数 * 输入：*输出： *返回： * 调用： *其它：用过以后，需要手动释放掉 *w 的内存空间 *调用示例： ret = hannWin(99,&w); */dspErrorStatus hannWin(dspUint_16 N,dspDouble *w)dspUint_16 n; dspDouble *ret;i

40、f(retret = (dspDouble *)malloc(N * sizeof(dspDouble);= NULL)return DSP_ERROR;*(ret + n) = 0.5 *for(n = 0; n < N; n+)( 1 - cos( 2 * PI * (dspDouble)n / (N - 1);*w = ret;return DSP_SUCESS;/* * 函数名： kaiserWin * 说明：计算 kaiserWin 窗函数 * 输入：*输出：*返回：*调用：besseliO第一类修正贝塞尔函数*其它：用过以后，需要手动释放掉 *w 的内存空间 *调用示

41、例： ret = kaiserWin(99, 5, &w); */dspErrorStatus kaiserWin(dspUint_16 N, dspDouble beta, dspDouble*w)dspUint_16 n; dspDouble *ret; dspDoubletheta;ret = (dspDouble *)malloc(N * sizeof(dspDouble);if(ret = NULL)return DSP_ERROR;theta = beta *for(n = 0; n < N; n+)sqrt( 1 - pow( ( (2 * (dsp

42、Double)n/(N -1) - 1),2 ) );*(ret + n) = (dspDouble)besseli(0, theta,BESSELI_K_LENGTH) / besseli(0, beta,BESSELI_K_LENGTH);*w = ret;return DSP_SUCESS;/* * 函数名： nuttalWin * 说明：计算 nuttalWin 窗函数 * 输入：放掉 *w 的内存空间 *输出： *返回： * 调用： *其它：用过以后，需要手动释调用示例： ret = nuttalWin(99,&w); */dspErrorStatus nuttalWin(dspUint_16 N,dspDouble *w)dspUint_16 n; dspDouble *ret;ret = (dspDouble *)malloc(N * sizeof(dspDouble);if(ret = NULL)return DSP_ERROR;*(ret + n)for(n = 0; n < N; n+)