参考文档
主要参考资料来源数字信号处理 基于计算机的方法 4版（中文版）和Matlab窗函数【Matlab帮助中心】加窗法两部分
分类
滤波器类型分类
- 低通
- 高通
- 带通
- 带阻
不同滤波器类型对应不同理想单位冲击函数hd(n)
窗函数分类
C语言实现
理想单位冲击响应函数
根据数字信号处理 基于计算机的方法 4版（中文版）中关于理想单位冲击响应函数介绍（如下图）实现。
- 低通
// wc为截至频率，N为窗口长度
double *IdealLp(double wc, int N) //低通
{
double M = (N - 1) / 2.0;
double *hd = (double *)malloc(N * sizeof(double));
for (int n = 0; n < N; n++)
{
if (n == M)
{
hd[n] = wc / PI;
}
else
{
hd[n] = sin(wc * (n - M)) / (PI * (n - M));
}
}
return hd;
}
- 高通
// wc为截至频率，N为窗口长度
double *IdealHp(double wc, int N) //高通
{
double M = (N - 1) / 2.0;
double *hd = (double *)malloc(N * sizeof(double));
for (int n = 0; n < N; n++)
{
if (n == M)
{
hd[n] = 1.0 - wc / PI;
}
else
{
hd[n] = -sin(wc * (n - M)) / (PI * (n - M));
}
}
return hd;
}
- 带通
// wcl为低频截至频率，wch为高频截至频率，N为窗口长度
double *IdealBp(double wcl, double wch, int N) // 带通
{
double M = (N - 1) / 2.0;
double *hd = (double *)malloc(N * sizeof(double));
hd = IdealLp(wcl, N);
hd = subtract(IdealLp(wch, N), hd, N);
return hd;
}
带通理想单位冲击响应函数可由低频和高频截至频率的低通理想单位冲击响应函数相减得到。
其中subtract为封装的double数组运算库中的减法函数，IdealLp为低通理想单位冲击响应函数。
- 带阻
// wcl为低频截至频率，wch为高频截至频率，N为窗口长度
double *IdealBs(double wcl, double wch, int N) // 带阻
{
double M = (N - 1) / 2.0;
double *hd = (double *)malloc(N * sizeof(double));
for (int n = 0; n < N; n++)
{
if (n == M)
{
hd[n] = 1.0 - (wch - wcl) / PI;
}
else
{
hd[n] = (sin(wcl * (n - M)) - sin(wch * (n - M))) / (PI * (n - M));
}
}
return hd;
}
至此，低通、高通、带通和带阻的理想单位冲击响应函数已全部实现。
窗函数
矩形窗
//全为1的数组
double *RectangleWindow(int n)
{
double *result = (double *)malloc(n * sizeof(double));
for (size_t i = 0; i < n; i++)
{
result[i] = 1.0;
}
return result;
}
三角窗
double *TriangWindow(int n)
{
double *result = (double *)malloc(n * sizeof(double));
if (n % 2 == 0) // even
{
for (size_t i = 1; i <= n; i++)
{
if (i >= 1 && i <= n / 2)
{
result[i - 1] = (2.0 * i - 1.0) / n;
}
else
{
result[i - 1] = 2.0 - (2.0 * i - 1.0) / n;
}
}
}
else if (n % 2 == 1) // odd
{
for (size_t i = 1; i <= n; i++)
{
if (i >= 1 && i <= (n + 1) / 2)
{
result[i - 1] = 2.0 * i / (n + 1.0);
}
else
{
result[i - 1] = 2.0 - 2.0 * i / (n + 1.0);
}
}
}
return result;
}
汉宁窗
double *HannWindow(int n, double *w)
{
double *result = (double *)malloc(n * sizeof(double));
for (size_t i = 0; i < n; i++)
{
result[i] = 0.5 - 0.5 * cos(w[i]);
}
return result;
}
其中参数w提前计算作为公用参数
double *w = (double *)malloc(n * sizeof(double));
for (size_t i = 0; i < n; i++)
{
w[i] = 2 * PI * i / (n - 1);
}
汉明窗
double *HammingWindow(int n, double *w)
{
double *result = (double *)malloc(n * sizeof(double));
for (size_t i = 0; i < n; i++)
{
result[i] = 0.54 - 0.46 * cos(w[i]);
}
return result;
}
布拉克曼窗
double *HammingWindow(int n, double *w)
{
double *result = (double *)malloc(n * sizeof(double));
for (size_t i = 0; i < n; i++)
{
result[i] = 0.54 - 0.46 * cos(w[i]);
}
return result;
}
平顶窗
double *FlattopWindow(int n, double *w)
{
double *result = (double *)malloc(n * sizeof(double));
double fta0 = 0.215578995;
double fta1 = 0.41663158;
double fta2 = 0.277263158;
double fta3 = 0.083578947;
double fta4 = 0.006947368;
for (size_t i = 0; i < n; i++)
{
result[i] = fta0 -
fta1 * cos(w[i]) +
fta2 * cos(2 * w[i]) -
fta3 * cos(3 * w[i]) +
fta4 * cos(4 * w[i]);
}
return result;
}
高斯窗
double *GaussianWindow(int n)
{
double *result = (double *)malloc(n * sizeof(double));
double alpha = 2.5;
double width = (n - 1) / 2.0;
for (size_t i = 0; i < n; i++)
{
double lenth = i - width;
result[i] = width == 0 ? 0 : exp(-0.5 * pow(alpha * lenth / width, 2));
}
return result;
}
其中alpha可以添加到参数中，这里选择固定值2.5,alpha选择原则：
- alpha < 1：如果您需要更窄的主瓣和更宽的副瓣，可以选择较小的 alpha。这对于需要更好的频率分辨率的应用可能有用，但会导致更慢的副瓣下降。
- 1 <= alpha <= 2：这是一种平衡选择，通常是最常用的。它提供了相对宽一些的主瓣和适度快速的副瓣下降。
- alpha > 2：如果您需要非常宽的主瓣和非常快速的副瓣下降，可以选择较大的 alpha。这对于某些特定的应用可能有用，但可能会降低频率分辨率。
凯赛窗
double *KaiserWindow(int n)
{
double *result = (double *)malloc(n * sizeof(double));
double beta = 7;
// β满足条件，α>50， β=0.1102（ α-8.7）；
// 50>=α>=21,β=0.5842(α−21)^0.4+0.07886(α−21);
// α<21,β=0;
for (size_t i = 0; i < n; i++)
{
// double x = beta * sqrt(1 - pow(1 - 2.0 * i / (n - 1), 2)); // old
double x = beta * sqrt(1 - pow((i-(n-1)/2.0)/((n-1)/2.0), 2)); //new
result[i] = I0function(x) / I0function(beta);
}
return result;
}
其中beta可以作为参数传入，一般根据滤波器性能指标进行计算，这里选择固定值7*；
方法I0function为第一类零阶修正贝塞尔函数，具体实现：
double I0function(double x)
{
double i0 = 0;
for (size_t i = 0; i < 25; i++)
{
i0 += pow(x, 2 * i) / pow(4, i) / pow(Fact(i), 2);
}
return i0;
}
double Fact(int n) // factorial阶乘
{
if (n == 0)
{
return 1.0;
}
double y = (double)n;
for (double m = n - 1; m > 0; m--)
{
y *= m;
}
return y;
}
其他窗函数
double *BartlettWindow(int n)
{
double *result = (double *)malloc(n * sizeof(double));
for (size_t i = 0; i < n; i++)
{
if (i <= (n - 1) / 2)
{
result[i] = 2.0 * i / (n - 1.0);
}
else
{
result[i] = 2.0 - 2.0 * i / (n - 1.0);
}
}
return result;
}
巴特利-汉宁窗
double *BartlettHanningWindow(int n)
{
double *result = (double *)malloc(n * sizeof(double));
for (size_t i = 0; i < n; i++)
{
result[i] = 0.62 - 0.48 * fabs((1.0 * i / (n - 1.0)) - 0.5) + 0.38 * cos(2.0 * PI * ((1.0 * i / (n - 1.0)) - 0.5));
}
return result;
}
布莱克曼-哈里斯窗
//type=1为对称性，=2为周期性
double *BlackmanHarrisWindow(int n, int type)
{
double *result = (double *)malloc(n * sizeof(double));
double a0 = 0.35875;
double a1 = 0.48829;
double a2 = 0.14128;
double a3 = 0.01168;
double m = type == 0 ? (n - 1.0) : n;
for (size_t i = 0; i < n; i++)
{
result[i] = a0 -
a1 * cos(2.0 * PI * i / m) +
a2 * cos(2.0 * 2.0 * PI * i / m) -
a3 * cos(3.0 * 2.0 * PI * i / m);
}
return result;
}
是不是很眼熟，和平顶窗公式类似，但是系数不一样。
Parzen
double *ParzenWindow(int n)
{
double *result = (double *)malloc(n * sizeof(double));
int index = 0;
for (int i = -(n - 1) / 2; i <= (n - 1) / 2; i++)
{
if (abs(i) >= 0 && abs(i) <= (n - 1) / 4)
{
result[index++] = 1 -
6.0 * pow(fabs(i) / (n / 2.0), 2) +
6.0 * pow(fabs(i) / (n / 2.0), 3);
}
else
{
result[index++] = 2 * pow((1.0 - fabs(i) / (n / 2.0)), 3);
}
}
return result;
}
纳托尔的布莱克曼-哈里斯窗
double *NuttallWindow(int n, int type)
{
double *result = (double *)malloc(n * sizeof(double));
double a0 = 0.3635819;
double a1 = 0.4891775;
double a2 = 0.1365995;
double a3 = 0.0106411;
double m = type == 0 ? (n - 1.0) : n;
for (size_t i = 0; i < n; i++)
{
result[i] = a0 -
a1 * cos(2.0 * PI * i / m) +
a2 * cos(2.0 * 2.0 * PI * i / m) -
a3 * cos(3.0 * 2.0 * PI * i / m);
}
return result;
}
又一个类似于平顶窗的窗函数，同样的，它们系数不一样。
切比雪夫窗
double *ChebyshevWindow(int n, double ripple) // 切比雪夫窗函数，ripple为削弱系数，默认100dB
{
}
没写完暂时空着，以后再加上
滤波器系数
通过上面代码可以获得对应长度的理想单位冲击响应函数系数hd(n)和窗函数系数w(n)，最后将它们对应相乘，再做归一化。
double *ApplyWindow(double *hd, int windowType, int n)
{
double *w = GetWindowCoefficient(windowType, n);
double *h = multiply(hd, w, n);//对应相乘
double sum_h = 0.0;
for (size_t i = 0; i < n; i++)
{
sum_h += h[i];
}
h = multiply_num(h, 1.0 / sum_h, n);
return h;
}
测试
- matlab代码实现
function CreateFirByWindow()
disp('Hello, Test Beginning!')
N = input('Please enter order!\n'); % Order
windowType=input('Please enter window type!(Rectangle 0,Hann 1,Hamming 2,Blackman 3,Flattop 4,Kaiser 5,Gaussian 6,Triang 7,Bartlett 8,Bartlett-Hann 9,Blackman-Harris(Symmetric) 10,Blackman-Harris(Periodic) 11,Parzen 12,Nuttall(Symmetric) 13,Nuttall(Periodic) 14)\n'); %Window Type
filterResponseType=input('Please enter filter response type!(LowPass 0, HighPass 1,BandPass 2,BandStop 3)\n'); %Filter Response Type
wc1 = input('Please enter wc1!(0<wc1<1)\n');
% wc1 = 1.0e-323;
wc2=0.6;
if filterResponseType==2||filterResponseType==3
wc2 = input('Please enter wc2!(0<wc2<1)\n');
else
wc2=0.6;
end
win = CreateWindow(windowType,N);
type='low';
switch filterResponseType
case 0
type='low';
case 1
type='high';
case 2
type='bandpass';
case 3
type='stop';
end
%
if filterResponseType<2
w=fir1(N, wc1, type, win,"scale");
else
w=fir1(N,[wc1,wc2],type,win,"scale");
end
w=vpa(w,20);
disp(w.'); %打印归一化前系数
function win = CreateWindow(windowType,order)
win=zeros(1,order+1);
switch windowType
case {0}
win=rectwin(order+1);
case {1}
win=hann(order+1);
case {2}
win=hamming(order+1);
case {3}
win=blackman(order+1);
case {4}
win=flattopwin(order+1);
case {5}
beta=7; %kaiser参数默认为7
win=kaiser(order+1,beta);
case {6}
alpha=2.5;%gaussian参数默认为2.5
win=gausswin(order+1,alpha);
case {7}
win=triang(order+1);
case {8}
win=bartlett(order+1);
case {9}
win=barthannwin(order+1);
case {10}
win=blackmanharris(order+1,"symmetric");
case {11}
win=blackmanharris(order+1,"periodic");
case {12}
win=parzenwin(order+1);
case {13}
win=nuttallwin(order+1,"symmetric");
case {14}
win =nuttallwin(order+1,"periodic");
case {15}
win=chebwin(order+1,100);
end
通过不同阶数（个位数，十位数，百位数，千位数）、奇偶性（阶数分为奇数和偶数）和极限值测试
阶数和奇偶性测试结果：
小数10位后有一定差别，总体结果基本一致
极限值测试结果
Matlab:
下限截至频率：1.0e-323,上限截至频率0.9999999999999999(16位小数)
C:
下限截至频率：1.0e-309，上限截至频率大于0.9999999999999999(大于16位小数)
基本通过，测试完成
5402
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
