目錄

• 源碼
• • WindowFunction.c
  • WindowFunction.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 whcui1987@163.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 whcui1987@163.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.5，24Hz, 相位為30°的方波信號進行FFT，有大小為2.5的直流偏置
1024點FFT

FFT代碼見適用于單片機的FFT快速傅里葉變換算法，51單片機都能用

無窗

漢寧窗

平頂窗

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的常见窗函数的C语言实现及其形状，适用于单片机、DSP作FFT运算的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。