登录建设厅网站的是企业锁吗,外贸公司英文,760关键词排名查询,北京中高端网站建设公司目录源码WindowFunction.cWindowFunction.h使用形状三角窗巴特利特窗巴特利特-汉宁窗布莱克曼窗布莱克曼-哈里斯窗博曼窗切比雪夫窗平顶窗高斯窗海明窗汉宁窗纳托尔窗Parzen窗矩形窗#xff08;模拟#xff09;效果无窗汉宁窗平顶窗平台#xff1a;Windows 10 20H2 Visual …
目录源码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 whcui1987163.com
*version 0.3
*data 31-Oct-2014
*brief 各种窗函数的C语言实现
*/#include WindowFunction.h
#include math.h
#include stdlib.h#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 (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 ((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 (*(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 ((N - 1) / 2); i){w[i] 2 * (double)(i 1) / (N 1);}for (i ((N - 1) / 2); i N; i){w[i] 2 * (double)(N - i) / (N 1);}}/*阶数为偶*/else{for (i 0; i (N / 2); i){w[i] (i i 1) * (double)1 / N;}for (i (N / 2); i N; 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 0 就是矩形窗*/if (r 0){retErrorStatus rectangularWin(N, retPtr);if (retErrorStatus DSP_ERROR)return DSP_ERROR;/*将数据拷出来以后释放调用的窗函数的空间*/memcpy(result, retPtr, (N * sizeof(dspDouble)));free(retPtr);}/*r 1 就是汉宁窗*/else if (r 1){retErrorStatus hannWin(N, retPtr);if (retErrorStatus DSP_ERROR)return DSP_ERROR;/*将数据拷出来以后释放调用的窗函数的空间*/memcpy(result, retPtr, (N * sizeof(dspDouble)));free(retPtr);}else{for (index 0; index N; index){alpha *(x index);if (alpha (r / 2)){*(result index) (dspDouble)(1 cos(2 * PI * (dspDouble)(alpha - (dspDouble)r / 2) / r)) / 2;}else if ((alpha (r / 2)) (alpha (1 - r / 2))){*(result index) 1;}else{*(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 (N - 1) / 2; n){w[n] 2 * (double)n / (N - 1);}for (n (N - 1) / 2; n N; 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 N; n){w[n] 0.62 - 0.48 * fabs(((double)n / (N - 1)) - 0.5) 0.38 * cos(2 * PI * (((double)n / (N - 1)) - 0.5));}for (n 0; n (N - 1) / 2; n){w[n] w[N - 1 - n];}}/*偶*/else{for (n 0; n N; n){w[n] 0.62 - 0.48 * fabs(((double)n / (N - 1)) - 0.5) 0.38 * cos(2 * PI * (((double)n / (N - 1)) - 0.5));}for (n 0; n N / 2; 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 N; n){w[n] 0.42 - 0.5 * cos(2 * PI * (double)n / (N - 1)) 0.08 * cos(4 * PI * (double)n / (N - 1));}return DSP_SUCESS;
}/*
*函数名blackManHarrisWin
*说明计算blackManHarrisWin窗函数
*输入
*输出
*返回
*调用
*调用示例ret blackManHarrisWin(99, w);
* minimum 4-term Blackman-harris window -- From Matlab
*/
dspErrorStatus blackManHarrisWin(uint16_t N, double w[])
{uint16_t n;for (n 0; n N; n){w[n] BLACKMANHARRIS_A0 - BLACKMANHARRIS_A1 * cos(2 * PI * (double)n / (N)) \BLACKMANHARRIS_A2 * cos(4 * PI * (double)n / (N)) - \BLACKMANHARRIS_A3 * cos(6 * PI * (double)n / (N));}return DSP_SUCESS;
}/*
*函数名bohmanWin
*说明计算bohmanWin窗函数
*输入
*输出
*返回
*调用
*调用示例ret bohmanWin(99, w);
*/
dspErrorStatus bohmanWin(uint16_t N, double w[])
{uint16_t n;double x;for (n 0; n N; n){x -1 n * (double)2 / (N - 1);/*取绝对值*/x x 0 ? x : (x * (-1));w[n] (1 - x) * cos(PI * x) (double)(1 / PI) * sin(PI * x);}return DSP_SUCESS;
}/*
*函数名chebyshevWin
*说明计算chebyshevWin窗函数
*输入
*输出
*返回
*调用
*调用示例ret chebyshevWin(99,100, w);
*/
dspErrorStatus chebyshevWin(uint16_t N, double r, double w[])
{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 (N 1) / 2; n){gama 1;for (index 1; index n; index){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 (N 1) / 2; n){w[n] (double)(w[n]) / theta;}/*填充另一半*/for (; n N; n){w[n] w[N - n - 1];}}else{/*计算一半的区间*/for (n 1; n (N 1) / 2; n){gama 1;for (index 1; index n; index){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 (N 1) / 2; n){w[n] (double)(w[n]) / theta;}/*填充另一半*/for (; n N; 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 N; n){w[n] FLATTOPWIN_A0 - FLATTOPWIN_A1 * cos(2 * PI * (double)n / (N - 1)) \FLATTOPWIN_A2 * cos(4 * PI * (double)n / (N - 1)) - \FLATTOPWIN_A3 * cos(6 * PI * (double)n / (N - 1)) \FLATTOPWIN_A4 * cos(8 * PI * (double)n / (N - 1));}return DSP_SUCESS;
}/*
*函数名gaussianWin
*说明计算gaussianWin窗函数
*输入
*输出
*返回
*调用
*调用示例ret gaussianWin(99,2.5, w);
*/
dspErrorStatus gaussianWin(uint16_t N, double alpha, double w[])
{uint16_t n;double k, beta, theta;for (n 0; n N; n){if ((N % 2) 1){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 N; n){w[n] 0.54 - 0.46 * cos(2 * PI * (double)n / (N - 1));}return DSP_SUCESS;
}/*
*函数名hannWin
*说明计算hannWin窗函数
*输入
*输出
*返回
*调用
*调用示例ret hannWin(99, w);
*/
dspErrorStatus hannWin(uint16_t N, double w[])
{uint16_t n;for (n 0; n N; n){w[n] 0.5 * (1 - cos(2 * PI * (double)n / (N - 1)));}return DSP_SUCESS;
}#if besseli_Flag
/*
*函数名kaiserWin
*说明计算kaiserWin窗函数
*输入
*输出
*返回
*调用besseli()第一类修正贝塞尔函数
*其它用过以后需要手动释放掉*w的内存空间
* 调用示例ret kaiserWin(99, 5, w);
*/
dspErrorStatus kaiserWin(dspUint_16 N, dspDouble beta, dspDouble **w)
{dspUint_16 n;dspDouble *ret;dspDouble theta;ret (dspDouble *)malloc(N * sizeof(dspDouble));if (ret NULL)return DSP_ERROR;for (n 0; n N; n){theta beta * sqrt(1 - pow(((2 * (dspDouble)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;
}
#endif/*
*函数名nuttalWin
*说明计算nuttalWin窗函数
*输入
*输出
*返回
*调用
*调用示例ret nuttalWin(99, w);
*/
dspErrorStatus nuttalWin(uint16_t N, double w[])
{uint16_t n;for (n 0; n N; n){w[n] NUTTALL_A0 - NUTTALL_A1 * cos(2 * PI * (double)n / (N - 1)) \NUTTALL_A2 * cos(4 * PI * (double)n / (N - 1)) - \NUTTALL_A3 * cos(6 * PI * (double)n / (N - 1));}return DSP_SUCESS;
}/*
*函数名parzenWin
*说明计算parzenWin窗函数
*输入
*输出
*返回
*调用
*调用示例ret parzenWin(99, w);
*/
dspErrorStatus parzenWin(uint16_t N, double w[])
{uint16_t n;double alpha, k;if ((N % 2) 1){for (n 0; n N; n){k n - (N - 1) / 2;alpha 2 * (double)fabs(k) / N;if (fabs(k) (N - 1) / 4){w[n] 1 - 6 * pow(alpha, 2) 6 * pow(alpha, 3);}else{w[n] 2 * pow((1 - alpha), 3);}}}else{for (n 0; n N; n){k n - (N - 1) / 2;alpha 2 * (double)fabs(k) / N;if (fabs(k) (double)(N - 1) / 4){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 N; n){w[n] 1;}return DSP_SUCESS;
}WindowFunction.h
/*
*file WindowFunction.h
*author Vincent Cui
*e-mail whcui1987163.com
*version 0.3
*data 31-Oct-2014
*brief 各种窗函数的C语言实现
*/#ifndef _WINDOWFUNCTION_H_
#define _WINDOWFUNCTION_H_#include stdint.h#define besseli_Flag 0 //缺少besseli函数
#define prod_Flag 0 //缺少prod函数
#define linSpace_Flag 0 //缺少linSpace函数#define BESSELI_K_LENGTH 10#define FLATTOPWIN_A0 0.215578995
#define FLATTOPWIN_A1 0.41663158
#define FLATTOPWIN_A2 0.277263158
#define FLATTOPWIN_A3 0.083578947
#define FLATTOPWIN_A4 0.006947368#define NUTTALL_A0 0.3635819
#define NUTTALL_A1 0.4891775
#define NUTTALL_A2 0.1365995
#define NUTTALL_A3 0.0106411#define BLACKMANHARRIS_A0 0.35875
#define BLACKMANHARRIS_A1 0.48829
#define BLACKMANHARRIS_A2 0.14128
#define BLACKMANHARRIS_A3 0.01168#define PI 3.14159265358979323846264338327950288419717 //定义圆周率值typedef enum
{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[]);
巴特利特窗
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[]); r 100 dB
平顶窗
dspErrorStatus flatTopWin(uint16_t N, double w[]);
高斯窗
dspErrorStatus gaussianWin(uint16_t N, double alpha, double w[]); alpha 2.5 alpha 8
海明窗
dspErrorStatus hammingWin(uint16_t N, double w[]);
汉宁窗
dspErrorStatus hannWin(uint16_t N, double w[]);
纳托尔窗
dspErrorStatus nuttalWin(uint16_t N, double w[]);
Parzen窗
dspErrorStatus parzenWin(uint16_t N, double w[]);
矩形窗
dspErrorStatus rectangularWin(uint16_t N, double w[]);
模拟效果
采样频率为100Hz 对一个振幅为2.524Hz, 相位为30°的方波信号进行FFT有大小为2.5的直流偏置 1024点FFT
FFT代码见适用于单片机的FFT快速傅里叶变换算法51单片机都能用
无窗 汉宁窗 平顶窗
