FIR数字滤波器设计(窗函数法) C语言实现
背景介绍：
理想滤波器在物理上是不可实现的，其单位脉冲响应是无限长、非因果的。窗函数法，就是从时域出发，用有限长、因果的单位脉冲响应h(n)去逼近理想滤波器的无限长、非因果的单位脉冲响应的方法。窗函数法又叫傅里叶级数法。
更多背景资料，请看数字信号处理（李永全），P175。
方法简介：
设N-1阶FIR数字滤波器的单位冲击响应为h(n),则传递函数H(z)为：
窗函数法的设计步骤如下：
1.根据给定的理想频率响应Hd(e^jw),利用傅里叶反变换，求出单位冲击响应hd(n):
2.将hd(n)乘以窗函数w(n),得到所要求的FIR滤波器系数h(n):
3.求卷积：
使用说明
- 子函数语句：
void firwin(int n, int band, int wn, int fs, double h[], double kaiser=0.0, double fln=0.0, double fhn=0.0);
- 形参说明
n：滤波器的阶数
band：滤波器的类型,取值1-4,分别为低通、带通、带阻、高通滤波器
wn：窗函数的类型，取值1-7，分别对应矩形窗、图基窗、三角窗、汉宁窗、海明窗、布拉克曼窗和凯塞窗
fs：采样频率
h：存放滤波器的系数
kaiser：beta值
fln：带通下边界频率
fhn：带通上边界频率
源代码
void firwin(int n, int band, int wn, int fs, double h[], double kaiser, double fln, double fhn)
{
int i;
int n2;
int mid;
double s;
double pi;
double wc1;
double wc2;
double beta;
double delay;
beta = kaiser;
pi = 4.0 * atan(1.0); //pi=PI;
if ((n % 2) == 0)/*如果阶数n是偶数*/
{
n2 = (n / 2) - 1;/**/
mid = 1;//
}
else
{
n2 = n / 2;//n是奇数,则窗口长度为偶数
mid = 0;
}
delay = n / 2.0;
wc1 = 2 * pi * fln;
wc2 = 2 * pi * fhn;
switch (band)
{
case 1:
{
for (i=0; i<=n2; ++i)
{
s = i - delay;
h[i] = (sin(wc1 * s / fs) / (pi * s)) * window(wn, n+1, i, beta);//低通,窗口长度=阶数+1，故为n+1
h[n - i] = h[i];
}
if (mid == 1)
{
h[n / 2] = wc1 / pi;//n为偶数时，修正中间值系数
}
break;
}
case 2:
{
for (i=0; i<=n2; i++)
{
s = i - delay;
h[i] = (sin(wc2 * s / fs) - sin(wc1 * s / fs)) / (pi * s);//带通
h[i] = h[i] * window(wn, n+1, i, beta);
h[n-i] = h[i];
}
if (mid == 1)
{
h[n / 2] = (wc2 - wc1) / pi;
}
break;
}
case 3:
{
for (i=0; i<=n2; ++i)
{
s = i - delay;
h[i] = (sin(wc1 * s / fs) + sin(pi * s) - sin(wc2 * s / fs)) / (pi * s);//带阻
h[i] = h[i] * window(wn, n+1, i, beta);
h[n - i] = h[i];
}
if (mid == 1)
{
h[n / 2] = (wc1 + pi - wc2) / pi;
}
break;
}
case 4:
{
for (i=0; i<=n2; i++)
{
s = i - delay;
h[i] = (sin(pi * s) - sin(wc1 * s / fs)) / (pi * s);//高通
h[i] = h[i] * window(wn, n+1, i, beta);
h[n-i] = h[i];
}
if (mid == 1)
{
h[n / 2] = 1.0 - wc1 / pi;
}
break;
}
}
}
//n：窗口长度 type：选择窗函数的类型 beta：生成凯塞窗的系数
static double window(int type, int n, int i, double beta)
{
int k;
double pi;
double w;
pi = 4.0 * atan(1.0); //pi=PI;
w = 1.0;
switch (type)
{
case 1:
{
w = 1.0; //矩形窗
break;
}
case 2:
{
k = (n - 2) / 10;
if (i <= k)
{
w = 0.5 * (1.0 - cos(i * pi / (k + 1))); //图基窗
}
if (i > n-k-2)
{
w = 0.5 * (1.0 - cos((n - i - 1) * pi / (k + 1)));
}
break;
}
case 3:
{
w = 1.0 - fabs(1.0 - 2 * i / (n - 1.0));//三角窗
break;
}
case 4:
{
w = 0.5 * (1.0 - cos( 2 * i * pi / (n - 1)));//汉宁窗
break;
}
case 5:
{
w = 0.54 - 0.46 * cos(2 * i * pi / (n - 1));//海明窗
break;
}
case 6:
{
w = 0.42 - 0.5 * cos(2 * i * pi / (n - 1)) + 0.08 * cos(4 * i * pi / (n - 1));//布莱克曼窗
break;
}
case 7:
{
w = kaiser(i, n, beta);//凯塞窗
break;
}
}
return(w);
}
static double kaiser(int i, int n, double beta)
{
double a;
double w;
double a2;
double b1;
double b2;
double beta1;
b1 = bessel0(beta);
a = 2.0 * i / (double)(n - 1) - 1.0;
a2 = a * a;
beta1 = beta * sqrt(1.0 - a2);
b2 = bessel0(beta1);
w = b2 / b1;
return(w);
}
static double bessel0(double x)
{
int i;
double d;
double y;
double d2;
double sum;
y = x / 2.0;
d = 1.0;
sum = 1.0;
for (i=1; i<=25; i++)
{
d = d * y / i;
d2 = d * d;
sum = sum + d2;
if (d2 < sum*(1.0e-8))
{
break;
}
}
return(sum);
}
得到系数之后，与输入信号求卷积即可！
