当前位置：首页 > news >正文

C++实现窗函数

news 2026/3/27 3:47:29

1、h

#ifndef WIN4C_H__
#define WIN4C_H__
#pragma once#define PI  3.141592653589793#define EPS  1e-8#include <vector>std::vector<double> rectangle(int);
std::vector<double> boxcar(int); //same as matlab, rectangle window.
std::vector<double> bartlett(int);
std::vector<double> hanning(int, bool symmetric=true);
std::vector<double> hamming(int);
std::vector<double> blackman(int);
std::vector<double> kaiser(int, double alpha=0.5);
std::vector<double> gausswin(int, double alpha=2.5);
std::vector<double> chebwin(int, double r = 100);
std::vector<double> triang(int);
std::vector<double> tukeywin(int, double r = 0.5);
double I0(double alpha);/*NOTE: Use the HANN function to get a Hanning window which has thefirst and last zero-weighted samples.
*/
std::vector<double> hann(int, bool symmetric = true);#endif

2、cpp

#include "win4c.h"
#include "fftw3.h"std::vector<double> boxcar(int n)
{return rectangle(n);
}std::vector<double> rectangle(int n)
{std::vector<double> win(n, 1);return win;
}std::vector<double> bartlett(int n)
{std::vector<double> win(n);for (int i = 0; i < (n + 1) / 2; ++i){win[i] = 2.0 * i / (n - 1.0);win[n - 1 - i] = win[i];}return win;
}std::vector<double> hanning(int n, bool symmetric)
{std::vector<double> win(n);if (symmetric){for (int i = 0; i < (n + 1) / 2; ++i){win[i] = 0.5 - 0.5 * cos(2 * PI * (i + 1) / (n + 1));win[n - 1 - i] = win[i];}}else{win.push_back(0);std::vector<double> tmp = hanning(n - 1, true);win.insert(win.begin() + 1, tmp.begin(), tmp.end());}return win;
}std::vector<double> hamming(int n)
{std::vector<double> win(n);for (int i = 0; i < (n + 1) / 2; ++i){win[i] = (0.54 - 0.46 * cos (2 * PI * i / (n - 1.0)));win[n - 1 - i] = win[i];}return win;
}std::vector<double> blackman(int n)
{std::vector<double> win(n);for (int i = 0; i < (n + 1) / 2; ++i){win[i] = (0.42 - 0.50 * cos(2 * PI * i / (n - 1.0)) + 0.08 * cos(2 * 2 * PI * i / (n - 1.0)));win[n - 1 - i] = win[i];}return win;
}std::vector<double> kaiser(int n, double alpha)
{std::vector<double> win(n);for (int i = 0; i < (n + 1) / 2; ++i){double beta = 2 * alpha * sqrt(i*(n - i - 1.0)) / (n - 1.0);win[i] = I0(beta) / I0(alpha);win[n - 1 - i] = win[i];}return win;
}std::vector<double> gausswin(int n, double alpha)
{std::vector<double> win(n);double center = (n - 1) / 2.0;for (int i = 0; i < (n + 1) / 2; ++i){double tmp = alpha * (i - center) / center;win[i] = exp(-0.5 * tmp * tmp);win[n - 1 - i] = win[i];}return win;
}std::vector<double> chebwin(int M, double r)
{std::vector<double> win(M);double order = M - 1.0;double beta = cosh(1.0 / order * acosh(pow(10, (abs(r) / 20.))));double *x = new double[M];for (int i = 0; i < M; i++){double temp = cos(PI * i / M)*beta;if (temp > 1){temp = cosh(order * acosh(temp));}            else if (temp < -1){temp = (2 * (M % 2) - 1) * cosh(order * acosh(-temp));}else{temp = cos(order * acos(temp));}            x[i] = temp;}fftw_complex *fftIn = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * M);fftw_complex *fftOut = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * M);fftw_plan plan = fftw_plan_dft_1d(M, fftIn, fftOut, FFTW_FORWARD, FFTW_MEASURE);double maxvalue(-1e8);if (M % 2)//M Ϊ����
    {for (int j = 0; j < M; j++){fftIn[j][0] = x[j];fftIn[j][1] = 0;}fftw_execute(plan);int n = (M + 1) / 2; win[n - 1] = fftOut[0][0];maxvalue = fftOut[0][0];for (int j = n - 1; j >= 1; j--){win[n - 1 - j] = fftOut[j][0];win[j + n - 1] = fftOut[j][0];if (maxvalue < fftOut[j][0]){maxvalue = fftOut[j][0];}                }}else{for (int j = 0; j < M; j++){fftIn[j][0] = x[j] * cos(PI / M * j);fftIn[j][1] = x[j] * sin(PI / M * j);}fftw_execute(plan);int n = M / 2 + 1;for (int j = n - 1; j >= 1; j--){win[n - 1 - j] = fftOut[j][0];win[j + n - 2] = fftOut[j][0];if (maxvalue < fftOut[j][0]) {maxvalue = fftOut[j][0];}                }}for (int j = 0; j < M; j++){win[j] = win[j] / maxvalue;}delete [] x;x = nullptr;fftw_destroy_plan(plan);fftw_free(fftIn);fftw_free(fftOut);return win;
}std::vector<double> triang(int n)
{std::vector<double> win(n);for (int i = 0; i < (n + 1) / 2; ++i){win[i] = (2.0 * ( i + 1) - 1 + n % 2) / (n + n % 2);win[n - 1 - i] = win[i];}return win;
}std::vector<double> hann(int n, bool symmetric)
{std::vector<double> win(n);win.push_back(0);std::vector<double> tmp = hanning(n - 2, symmetric);win.insert(win.begin() + 1, tmp.begin(), tmp.end());win.push_back(0);return win;
}std::vector<double> tukeywin(int n, double r)
{if (r <= 0){return boxcar(n);}if (r >= 1){return hann(n);}std::vector<double> win(n);std::vector<double> t(n);double step = 1.0 / (n - 1);for (int i = 0; i < n - 1; i++){t[i] = i * step;}t[n - 1] = 1;double per = r / 2.0;int tl = floor(per * (n - 1)) + 1;int th = n - tl + 1;//Window is defined in three sections : taper, constant, taper//    w = [((1 + cos(pi / per * (t(1:tl) - per))) / 2);  ones(th - tl - 1, 1); ((1 + cos(pi / per * (t(th:end) - 1 + per))) / 2)];for (int i = 0; i < tl; i++){win[i] = (1 + cos(PI / per * (t[i] - per))) / 2.0;}for (int i = tl; i < th - 1; i++){win[i] = 1;}for (int i = th - 1; i < n; i++){win[i] = (1 + cos(PI / per * (t[i] - 1 + per))) / 2.0;}return win;
}
double I0(double alpha)
{double dNew;double K = alpha / 2.0;const int MAXTERM = 25 + 1;double J = 1.0, dOld = 1.0;bool converge = false;// Use series expansion definition of Bessel.for (int i = 1; i < MAXTERM; ++i){J *= K / i;dNew = dOld + J * J;if ((dNew - dOld) < EPS){converge = true;break;}dOld = dNew;}if (!converge)return 0;return dNew;
}