背景介紹：
理想濾波器在物理上是不可實現的，其單位脈沖響應是無限長、非因果的。窗函數法，就是從時域出發，用有限長、因果的單位脈沖響應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);
}
得到系數之后，與輸入信號求卷積即可！