常见窗函数的C语言实现及其形状，适用于单片机、DSP作FFT运算

您所在的位置：网站首页 › 巴特利家具 › 常见窗函数的C语言实现及其形状，适用于单片机、DSP作FFT运算

常见窗函数的C语言实现及其形状，适用于单片机、DSP作FFT运算

2023-08-13 08:55| 来源: 网络整理| 查看: 265

目录源码WindowFunction.cWindowFunction.h 使用形状三角窗巴特利特窗巴特利特-汉宁窗布莱克曼窗布莱克曼-哈里斯窗博曼窗切比雪夫窗平顶窗高斯窗海明窗汉宁窗纳托尔窗Parzen窗矩形窗（模拟）效果无窗汉宁窗平顶窗

平台：Windows 10 20H2 Visual Studio 2015 Python 3.8.12 (default, Oct 12 2021, 03:01:40) [MSC v.1916 64 bit (AMD64)] :: Anaconda, Inc. on win32

原作见窗函数的C语言实现 —— Vincent.Cui

可配合C语言实现的FFT与IFFT源代码，不依赖特定平台使用

原代码大量使用了动态内存分配，考虑到部分单片机的限制，我把它们又改回了数组传参的形式。

由于缺少besseli、prod和linSpace函数，有三个窗函数暂时被我用条件编译注释掉了。

源码 WindowFunction.c /* *file WindowFunction.c *author Vincent Cui *e-mail [email protected] *version 0.3 *data 31-Oct-2014 *brief 各种窗函数的C语言实现 */ #include "WindowFunction.h" #include #include #if prod_Flag /*函数名：taylorWin *说明：计算泰勒窗。泰勒加权函数 *输入： *输出： *返回： *调用：prod()连乘函数 *其它：用过以后，需要手动释放掉*w的内存空间 * 调用示例：ret = taylorWin(99, 4, 40, &w); 注意此处的40是正数表示-40dB */ dspErrorStatus taylorWin(dspUint_16 N, dspUint_16 nbar, dspDouble 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) * (nbar - 1)); if (retDspDouble == NULL) return DSP_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), 2)); for (i = 1; i theta *= 1 - (dspDouble)(i * i) / (beta * beta * (A + (j - 0.5) * (j - 0.5))); } *(sf + i - 1) = alpha * (dspDouble)theta / (*(sf + i - 1)); } /*奇数阶*/ if ((N % 2) == 1) { for (i = 0; i alpha += (*(sf + j - 1)) * cos(2 * PI * j * (dspDouble)(i - ((N - 1) / 2)) / N); } *(result + i) = 1 + 2 * alpha; } } /*偶数阶*/ else { for (i = 0; i alpha += (*(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; } #endif /* *函数名：triangularWin *说明：计算三角窗函数 *输入： *输出： *返回： *调用： *调用示例：ret = triangularWin(99, w); */ dspErrorStatus triangularWin(uint16_t N, double w[]) { uint16_t i; /*阶数为奇*/ if ((N % 2) == 1) { for (i = 0; i w[i] = 2 * (double)(N - i) / (N + 1); } } /*阶数为偶*/ else { for (i = 0; i w[i] = w[N - 1 - i]; } } return DSP_SUCESS; } #if linSpace_Flag /* *函数名：tukeyWin *说明：计算tukey窗函数 *输入： *输出： *返回：linSpace() *调用： *其它：用过以后，需要手动释放掉*w的内存空间 * 调用示例：ret = tukeyWin(99, 0.5, &w); */ dspErrorStatus tukeyWin(dspUint_16 N, dspDouble r, dspDouble **w) { dspErrorStatus retErrorStatus; dspUint_16 index; dspDouble *x, *result, *retPtr; dspDouble alpha; retErrorStatus = linSpace(0, 1, N, &x); if (retErrorStatus == DSP_ERROR) return DSP_ERROR; result = (dspDouble *)malloc(N * sizeof(dspDouble)); if (result == NULL) return DSP_ERROR; /*r retErrorStatus = hannWin(N, &retPtr); if (retErrorStatus == DSP_ERROR) return DSP_ERROR; /*将数据拷出来以后，释放调用的窗函数的空间*/ memcpy(result, retPtr, (N * sizeof(dspDouble))); free(retPtr); } else { for (index = 0; index *(result + index) = (dspDouble)(1 + cos(2 * PI * (dspDouble)(alpha - (dspDouble)r / 2) / r)) / 2; } else if ((alpha >= (r / 2)) && (alpha *(result + index) = (dspDouble)(1 + cos(2 * PI * (dspDouble)(alpha - 1 + (dspDouble)r / 2) / r)) / 2; } } } free(x); *w = result; return DSP_SUCESS; } #endif /* *函数名：bartlettWin *说明：计算bartlettWin窗函数 *输入： *输出： *返回： *调用： *调用示例：ret = bartlettWin(99, w); */ dspErrorStatus bartlettWin(uint16_t N, double w[]) { uint16_t n; for (n = 0; n w[n] = 2 - 2 * (double)n / ((N - 1)); } return DSP_SUCESS; } /* *函数名：bartLettHannWin *说明：计算bartLettHannWin窗函数 *输入： *输出： *返回： *调用： *调用示例：ret = bartLettHannWin(99, w); */ dspErrorStatus bartLettHannWin(uint16_t N, double w[]) { uint16_t n; /*奇*/ if ((N % 2) == 1) { for (n = 0; n w[n] = w[N - 1 - n]; } } /*偶*/ else { for (n = 0; n w[n] = w[N - 1 - n]; } } return DSP_SUCESS; } /* *函数名：blackManWin *说明：计算blackManWin窗函数 *输入： *输出： *返回： *调用： *调用示例：ret = blackManWin(99, w); */ dspErrorStatus blackManWin(uint16_t N, double w[]) { uint16_t n; for (n = 0; n uint16_t n; for (n = 0; n uint16_t n; double x; for (n = 0; n uint16_t n, index; double x, alpha, beta, theta, gama; /*10^(r/20)*/ theta = pow((double)10, (double)(fabs(r) / 20)); beta = pow(cosh(acosh(theta) / (N - 1)), 2); alpha = 1 - (double)1 / beta; if ((N % 2) == 1) { /*计算一半的区间*/ for (n = 1; n x = index * (double)(N - 1 - 2 * n + index) / ((n - index) * (n + 1 - index)); gama = gama * alpha * x + 1; } w[n] = (N - 1) * alpha * gama; } theta = w[(N - 1) / 2]; w[0] = 1; for (n = 0; n w[n] = w[N - n - 1]; } } else { /*计算一半的区间*/ for (n = 1; n x = index * (double)(N - 1 - 2 * n + index) / ((n - index) * (n + 1 - index)); gama = gama * alpha * x + 1; } w[n] = (N - 1) * alpha * gama; } theta = w[(N / 2) - 1]; w[0] = 1; for (n = 0; n w[n] = w[N - n - 1]; } } return DSP_SUCESS; } /* *函数名：flatTopWin *说明：计算flatTopWin窗函数 *输入： *输出： *返回： *调用： *调用示例：ret = flatTopWin(99, w); */ dspErrorStatus flatTopWin(uint16_t N, double w[]) { uint16_t n; for (n = 0; n uint16_t n; double k, beta, theta; for (n = 0; n k = n - (N - 1) / 2; beta = 2 * alpha * (double)k / (N - 1); } else { k = n - (N) / 2; beta = 2 * alpha * (double)k / (N - 1); } theta = pow(beta, 2); w[n] = exp((-1) * (double)theta / 2); } return DSP_SUCESS; } /* *函数名：hammingWin *说明：计算hammingWin窗函数 *输入： *输出： *返回： *调用： *调用示例：ret = hammingWin(99, w); */ dspErrorStatus hammingWin(uint16_t N, double w[]) { uint16_t n; for (n = 0; n uint16_t n; for (n = 0; n dspUint_16 n; dspDouble *ret; dspDouble theta; ret = (dspDouble *)malloc(N * sizeof(dspDouble)); if (ret == NULL) return DSP_ERROR; for (n = 0; n uint16_t n; for (n = 0; n uint16_t n; double alpha, k; if ((N % 2) == 1) { for (n = 0; n w[n] = 1 - 6 * pow(alpha, 2) + 6 * pow(alpha, 3); } else { w[n] = 2 * pow((1 - alpha), 3); } } } else { for (n = 0; n w[n] = 1 - 6 * pow(alpha, 2) + 6 * pow(alpha, 3); } else { w[n] = 2 * pow((1 - alpha), 3); } } } return DSP_SUCESS; } /* *函数名：rectangularWin *说明：计算rectangularWin窗函数 *输入： *输出： *返回： *调用： *调用示例：ret = rectangularWin(99, w); */ dspErrorStatus rectangularWin(uint16_t N, double w[]) { uint16_t n; for (n = 0; n DSP_ERROR = 0, DSP_SUCESS, }dspErrorStatus; dspErrorStatus triangularWin(uint16_t N, double w[]); dspErrorStatus bartlettWin(uint16_t N, double w[]); dspErrorStatus bartLettHannWin(uint16_t N, double w[]); dspErrorStatus blackManWin(uint16_t N, double w[]); dspErrorStatus blackManHarrisWin(uint16_t N, double w[]); dspErrorStatus bohmanWin(uint16_t N, double w[]); dspErrorStatus chebyshevWin(uint16_t N, double r, double w[]); dspErrorStatus flatTopWin(uint16_t N, double w[]); dspErrorStatus gaussianWin(uint16_t N, double alpha, double w[]); dspErrorStatus hammingWin(uint16_t N, double w[]); dspErrorStatus hannWin(uint16_t N, double w[]); dspErrorStatus nuttalWin(uint16_t N, double w[]); dspErrorStatus parzenWin(uint16_t N, double w[]); dspErrorStatus rectangularWin(uint16_t N, double w[]); #if besseli_Flag dspErrorStatus kaiserWin(uint16_t N, double beta, double w[]); #endif #if prod_Flag dspErrorStatus taylorWin(uint16_t N, uint16_t nbar, double sll, double w[]); #endif #if linSpace_Flag dspErrorStatus tukeyWin(uint16_t N, double r, double w[]); #endif #endif 使用

FFT_N为存放时域数值的数组大小，一般与所采用的FFT点数一致

double Window[FFT_N] = {0}; bartLettHannWin(FFT_N, Window);

调用后Window[]内便存入了窗函数的系数，再将这些系数与存放时域数值的数组元素一一相乘应该就行。

形状

以下均为1024点生成的窗函数形状，数据由VS2015产生，图像由 python3 绘制。

三角窗

dspErrorStatus triangularWin(uint16_t N, double w[]); 在这里插入图片描述