[Matlab]FIR滤波器设计：(基本窗函数FIR滤波器设计)

​ IIR滤波器主要设计方法先设计一个模拟低通滤波器，然后把它转化为形式上的数字滤波器。但对于FIR滤波器来说，设计方法的关键要求之一就是保证线性相位条件。而IIR滤波器的设计方法中只对幅值特性进行设计，因此无法保证相位。所以FIR滤波器的设计需要采用完全不同的方法。FIR滤波器的设计方法主要有窗函数法、频率采样法、切比雪夫逼近法等。

​ 由于理想滤波器在边界频率处不连续，故时域信号$h_d(n)$一定是无限时宽的，也是非因果的。所以理想低通滤波器是无法实现的。如果实现一个具有理想线性相位特性的滤波器，其幅值特性只能采用逼近理想幅频特性的方法实现。如果对时域信号$h_d(n)$进行截取，并保证截取过程中序列保持对称，而且截取长度为N,则对称点为$\alpha =\frac{1}{2}\ast (N-1)$。截取后序列为$h(n)$，侧$h(n)$可用下式子表示：
$$h(n)=h_d(n)\ast w(n)$$
式中，$w(n)$为截取函数，又称为穿函数。如果窗函数为矩形序列，则称之为矩形窗。窗函数有多种形式，为保证加窗后系统的线性相位特性，必须保证加窗后序列关于$\alpha =\frac{1}{2}\ast (N-1)$点对称。常用的窗函数有矩形窗、汉宁窗、汉明窗、布莱克曼窗、凯塞窗。窗函数设计法的基本思想是用一个长度为N的序列$h(n)$代替$h_d(n)$作为实际设计的滤波器的单位脉冲响应。这种设计法成为窗函数设计法。显然在保证$h(n)$对称性的前提下，窗函数长度N越长，则$h(n)$越接近$h_d(n)$。但是误差是肯定存在的，这种误差成为截断误差。

​ 要确定如何设计一个FIR滤波器，首先得对加窗后的理想滤波器特性变化进行分析，并研究减少由截断引起的误差的途径，从而提出更好的滤波器设计方案。对于调整窗口长度可以有效地控制过渡带的宽度，但减少带内波动以及加大阻带衰减没有作用。所以必须挑选最为合适的窗函数对理想滤波器进行截取。下面简单介绍几种窗函数。一个实际的滤波器的单位脉冲响应可表示为：
$$h(n)=h_d(n)\ast w(n)$$
几种常见窗函数（只给出了窗函数的定义和幅度特性）：
$$W(e^{j\ast w})=W(w)(e^{-j\alpha w})$$

矩形窗FIR滤波器设计:

矩形窗的窗函数为：
$$w_R(n)=R_N(n)$$
幅度函数为：
$$R_N(w)=\frac{sin(wN/2)}{sin(w/2)}$$
它的主瓣宽度为$4\pi /N$,第一瓣比主瓣地13dB.

在Matlab中，实现矩形窗的函数为boxcar和recttwin ,其调用格式如下：

w=boxcar(N);
w=recttwin(N);
%%%显示窗函数的GUI工具
n = 60;
w = rectwin(n);
wvtool(w)%显示窗函数的GUI工具
%还提供了显示窗函数的GUI工具，如wvtool可以显示用来显示窗的形状和频域图形，wintool可以打开窗设计和分析工具，如运行
wvtool(hamming(64),hann(64),gausswin(64))
%%可以对比汉明窗、汉宁窗和高斯窗

其中，N是窗函数的长度，返回值w是一个N阶的向量，它的元素有窗函数的值组成。其中w=boxcar(N)等价于w=ones(N,1)。

案例分析：

利用矩形窗设计FIR带阻滤波器，主要参数如下:

给定抽样频率为$\Omega s=2\ast pi\ast 1.5\ast 10^4(rad/sec)$,

通带截至频率为$\Omega p1=2\ast pi\ast 0.75\ast 10^3(rad/sec),\Omega p2=2\ast pi\ast 6\ast 10^3(rad/sec)$
阻带截至频率为$\Omega st1=2\ast pi\ast 2.25\ast 10^3(rad/sec),\Omega st2=2\ast pi\ast 1.5\ast 10^3(rad/sec)$
阻带衰减${\delta }_2>=50dB$。

%%%%调用子程序1：
function hd=ideal_bs(Wcl,Wch,m); 
alpha=(m-1)/2; 
n=[0:1:(m-1)];
m=n-alpha+eps; 
hd=[sin(m*pi)+sin(Wcl*m)-sin(Wch*m)]./(pi*m)
%%%%调用子程序2：
function[db,mag,pha,w]=freqz_m2(b,a)
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))'; w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
%%%%运行MATLAB源代码如下：
clear all;
Wph=3*pi*6.25/15;%通带频率
Wpl=3*pi/15;
Wsl=3*pi*2.5/15;%阻带频率
Wsh=3*pi*4.75/15;
tr_width=min((Wsl-Wpl),(Wph-Wsh));%%过渡带带宽
%过渡带宽度
N=ceil(4*pi/tr_width);					%滤波器长度
n=0:1:N-1;
Wcl=(Wsl+Wpl)/2;						%理想滤波器的截止频率
Wch=(Wsh+Wph)/2;
hd=ideal_bs(Wcl,Wch,N);				%理想滤波器的单位冲击响应
w_ham=(boxcar(N))';
string=['矩形窗','N=',num2str(N)];
h=hd.*w_ham;						%截取取得实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2(h,[1]);
%计算实际滤波器的幅度响应
delta_w=2*pi/1000;
subplot(241);
stem(n,hd);
title('理想脉冲响应hd(n)')
axis([-1,N,-0.5,0.8]);
xlabel('n');ylabel('hd(n)');
grid on
subplot(242);
stem(n,w_ham);
axis([-1,N,0,1.1]);
xlabel('n');ylabel('w(n)');
text(1.5,1.3,string);
grid on
subplot(243);
stem(n,h);title('实际脉冲响应h(n)');
axis([0,N,-1.4,1.4]);
xlabel('n');ylabel('h(n)');
grid on
subplot(244);
plot(w,pha);title('相频特性');
axis([0,3.15,-4,4]);
xlabel('频率（rad）');ylabel('相位（Φ）');
grid on
subplot(245);
plot(w/pi,db);title('幅度特性（dB）');
axis([0,1,-80,10]);
xlabel('频率（pi）');ylabel('分贝数');
grid on
subplot(246);
plot(w,mag);title('频率特性')
axis([0,3,0,2]);
xlabel('频率（rad）');ylabel('幅值');
grid on
fs=15000;
t=(0:100)/fs;
x=cos(2*pi*t*750)+cos(2*pi*t*3000)+cos(2*pi*t*6100);
q=filter(h,1,x);
[a,f1]=freqz(x);
f1=f1/pi*fs/2;
[b,f2]=freqz(q);
f2=f2/pi*fs/2;
subplot(247);
plot(f1,abs(a));
title('输入波形频谱图');
xlabel('频率');ylabel('幅度')
grid on
subplot(248);
plot(f2,abs(b));
title('输出波形频谱图');
xlabel('频率');ylabel('幅度')
grid on

汉宁窗FIR滤波器设计:

汉宁窗（hanning window）又称升余弦窗，窗函数为：
$$w_{Hn}(n)=0.5[1-cos(\frac{2\pi \ast n}{N-1})]\ast R_N(n)$$
幅值函数为：
$$W_{Hn}(w)=0.5W_R(w)+0.25[W_R(w-\frac{2\pi }{N-1})+W_R(w+\frac{2\pi }{N-1})]$$
汉宁窗幅度函数由3部分相加而成，其结果是使主瓣集中了更多能量，而旁瓣3部分相加时相互抵消而变小，其代价是主瓣宽度增加到$8\pi /N$。第一瓣比主瓣低31dB,阻带衰减加大。

在Matlab中，实现汉宁窗的函数为hanning和barthannwin ,其调用格式如下：

w=hanning(N)
w=barthannwin(N)

案例1：绘制50个点的汉宁窗。

N=49;n=1:N;
wdhn=hanning(N);	
figure(3);
stem(n,wdhn,'.');
grid on
axis([0,N,0,1.1]);
title('50点汉宁窗');
ylabel('W(n)');
xlabel('n');
title('50点汉宁窗');

案例2：已知连续信号为$x(t)=cos(2\pi f_{1}t)+0.15cos(2\pi f_{2}t)$，其中$f_1=100Hz,f_2=150Hz$。若抽样频率$f_{s}am=600Hz$d对信号进行抽样，利用不同宽度N的矩形截断该序列，N取40，观察不同的窗对普分析结果的影响。

N=40; 
L=512;
f1=100;f2=150;fs=600;
ws=2*pi*fs;
t=(0:N-1)*(1/fs);
x=cos(2*pi*f1*t)+0.25*sin (2*pi*f2*t);
wh=boxcar(N)';
x=x.*wh;
subplot(221);stem(t,x);
title('加矩形窗时域图');
xlabel('n');ylabel('h(n)')
grid on
W=fft(x,L);
f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi); 
subplot(222);
plot(f,abs(fftshift(W)))
title('加矩形窗频域图');
xlabel('频率');ylabel('幅度')
grid on
x=cos(2*pi*f1*t)+0.15*cos(2*pi*f2*t); 
wh=hanning(N)';
x=x.*wh;
subplot(223);stem(t,x);
title('加汉宁窗时域图');
xlabel('n');ylabel('h(n)')
grid on
W=fft(x,L);
f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);
subplot(224);
plot(f,abs(fftshift(W)))
title('加汉宁窗频域图');
xlabel('频率');ylabel('幅度')
grid on

用汉宁窗对谐波信号进行分析：

clear; 
% 原始数据：直流:0V; 基波：49.5Hz,100V,10deg; HR2:0.5V,40deg;
hr0=0;f1=50.1; 
hr(1)=25*sqrt(2);deg(1)=10; 
hr(2)=0;deg(2)=0; 
hr(3)=1.755*sqrt(2);deg(3)=40; 
hr(4)=0;deg(4)=0; 
hr(5)=0.885*sqrt(2);deg(5)=70; 
hr(6)=0;deg(6)=0; 
hr(7)=1.125;deg(7)=110; 
M=7;f=[1:M]*f1;							%设定频率 
% 采样 
fs=10000; 
N=2048;									% 约10个周期 
T=1/fs; 
n=[0:N-1];t=n*T; 
x=zeros(size(t)); 
for k=1:M 
    x=x+hr(k)*cos(2*pi*f(k)*t+deg(k)*pi/180); 
end 
%分析： 
w=0.5-0.5*cos(2*pi*n/N);
Xk=fft(x.*w); 
amp=abs(Xk(1:N/2))/N*2;						%幅频 
pha=angle(Xk(1:N/2))/pi*180;					%相频 
for k=1:N/2 
    if(amp(k)<0.01) pha(k)=0;  %当谐波<10mV时，其相位=0 
    end 
    if(pha(k)<0) pha(k)=pha(k)+360;%调整到0-360度 
    end 
end 
fmin=fs/N; 
xaxis=fmin*n(1:N/2);
%横坐标为Hz 
kx=round([1:M]*50/fmin);
%各次谐波对应的下标(从0开始) 
for m=1:M 
    km(m)=searchpeaks(amp,kx(m)+1);			%km为谱峰(从1开始) 
    if(amp(km(m)+1)<amp(km(m)-1)) 
        km(m)=km(m)-1; 
    end 
    beta(m)=amp(km(m)+1)./amp(km(m)); 
    delta(m)=(2*beta(m)-1)./(1+beta(m)); 
end 
fx=(km-1+delta)*fmin;							%估计频率 
hrx=amp(km)*2.*pi.*delta.*(1-delta.*delta)./sin(pi*delta);
degx=pha(km)-delta.*180/N*(N-1);				%估计相位 
degx=mod(degx,360);							%调整到0-360度 
efx=(fx-f)./f*100;								%频率误差 
ehr=(hrx-hr)./hr*100;							%幅度误差 
edeg=(degx-deg);								%相位误差 
% 结果输出： 
subplot(2,2,1);
%画出采样序列 
plot(t,x); 
hold on; 
plot(t,x.*w,'r');
%加窗波形 
hold off; 
xlabel('x(k)'); 
title('原信号和加窗信号 '); 
subplot(2,2,2);
%画出FFT分析结果 
stem(xaxis,amp,'.r'); 
xlabel('频率'); 
title('幅频结果'); 
subplot(2,2,4); 
stem(xaxis,pha,'.r'); 
xlabel('角频率'); 
title('相频结果'); 
subplot(2,2,3); 
stem(ehr); 
title('幅度误差(%)'); 
%文本输出 
fid=fopen('result.txt','w'); 
fprintf(fid,'原始数据：f1=%6.1fHz, N=%.f,  fs=%.f \r\n\r\n',f1,N,fs); 
fprintf(fid,'谐波次数      1      2      3      4      5      6     7\r\n'); 
fprintf(fid,'设定频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',f); 
fprintf(fid,'估计频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',fx); 
fprintf(fid,'误差(%%)  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n',efx); 
fprintf(fid,'设定幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',hr); 
fprintf(fid,'估计幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',hrx); 
fprintf(fid,'误差(%%)  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n',ehr); 
fprintf(fid,'设定相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n',deg); 
fprintf(fid,'估计相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n',degx); 
fprintf(fid,'误差(度)  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n\r\n',edeg); 
%其他数据 
fprintf(fid,'谱峰位置理论值：\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',[1:M]*f1/fmin); 
fprintf(fid,'谱峰位置估计值：\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',km-1+delta); 
fprintf(fid,'误差（%%）\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',((km-1+delta)-[1:M]*f1/fmin)./([1:M]*f1/fmin)*100); 
fprintf(fid,'delta     ：\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',delta); 
fclose(fid);

%运行过程中的调用子程序：
function index1=searchpeaks(x,index) 
%在数组中寻找最大值对应的下标 
%x为数组，index 为给定的下标(index不能取最前或最后2个下标)，在前后2个数中（共5个数）查找最大值和紧邻的次最大值 
% indexmax 返回两个谱峰位置中的前一个谱峰对应的下标 
index1=index-2; 
for k=-1:2 
    if(x(index+k)>x(index1)) 
        index1=index+k; 
    end 
end 
if x(index1-1)>x(index1+1) 
    index1=index1-1; 
end

#result.txt输出结果
原始数据：f1=  50.1Hz, N=2048,  fs=10000 

谐波次数      1      2      3      4      5      6     7
设定频率 50.100 100.200 150.300 200.400 250.500 300.600 350.700
估计频率 50.100 78.819 150.302 181.252 250.499 279.138 350.701
误差(%)  -0.000 -21.338  0.001 -9.555 -0.000 -7.140  0.000

设定幅值 35.355  0.000  2.482  0.000  1.252  0.000  1.125
估计幅值 35.356  0.046  2.482  0.002  1.252  0.002  1.125
误差(%)   0.001    Inf  0.009    Inf  0.004    Inf  0.003

设定相位  10.00   0.00  40.00   0.00  70.00   0.00 110.00
估计相位  10.03  31.67  39.97 338.35  70.06 329.86 110.05
误差(度)    0.03  31.67  -0.03 338.35   0.06 329.86   0.05

谱峰位置理论值：
 10.2605 20.5210 30.7814 41.0419 51.3024 61.5629 71.8234
谱峰位置估计值：
 10.2605 16.1421 30.7818 37.1203 51.3022 57.1675 71.8235
误差（%）
 -0.0002 -21.3385 0.0012 -9.5551 -0.0004 -7.1396 0.0002
delta     ：
 0.2605 0.1421 0.7818 0.1203 0.3022 0.1675 0.8235

汉明窗FIR滤波器设计:

汉宁窗（hanming window）又称改进升余弦窗，窗函数为：
$$w_{Hn}(n)=[0.54-0.46\ast cos(\frac{2\pi \ast n}{N-1})]\ast R_N(n)$$
幅值函数为：
$$W_{Hn}(w)=0.54W_R(w)+0.23[W_R(w-\frac{2\pi }{N-1})+0.23W_R(w+\frac{2\pi }{N-1})]$$
汉明窗主瓣宽度与汉宁窗相同，$8\pi /N，99.96\%$的能量集中在主瓣，第一瓣比主瓣低41dB。

在Matlab中，实现汉宁窗的函数为hanming ,其调用格式如下：

w=hanming(N)

案例分析1：设计一个汉明窗低通滤波器：

%语音信号设计一个汉明窗低通滤波器：
[x,FS,bits]=wavread('C:\Windows\Media\Windows Ringout');
x=x(:,1);
figure(1);
subplot(211);plot(x);
title('语音信号时域波形图')
xlabel('n');ylabel('h(n)')
grid on
y=fft(x,1000);
f=(FS/1000)*[1:1000];  
subplot(212);
plot(f(1:300),abs(y(1:300)));
title('语音信号频谱图');
xlabel('频率');ylabel('幅度')
grid on
%产生噪声信号并加到语音信号
t=0:length(x)-1;
zs0=0.05*cos(2*pi*10000*t/1024);
zs=[zeros(0,20000),zs0];
figure(2);
subplot(211)
plot(zs)
title('噪声信号波形');
xlabel('n');ylabel('h(n)')
grid on
zs1=fft(zs,1200);
subplot(212)
plot(f(1:600),abs(zs1(1:600)));
title('噪声信号频谱');
xlabel('频率');ylabel('幅度')
grid on
x1=x+zs';
%sound(x1,FS,bits);
y1=fft(x1,1200);
figure(3);
subplot(211);plot(x1);
title('加入噪声后的信号波形');
xlabel('n');ylabel('h(n)')
grid on
subplot(212);
plot(f(1:600),abs(y1(1:600)));
title('加入噪声后的信号频谱');
xlabel('频率');ylabel('幅度')
grid on
%滤波
fp=7500;
fc=8500; 
wp=2*pi*fp/FS;
ws=2*pi*fc/FS;
Bt=ws-wp; 
N0=ceil(6.2*pi/Bt);     
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;         
hn=fir1(N-1,wc,hamming(N)); 
X=conv(hn,x);           
X1=fft(X,1200);
figure(4);
subplot(211);
plot(X);
title('滤波后的信号波形');
xlabel('n');ylabel('h(n)')
grid on
subplot(212);
plot(f(1:600),abs(X1(1:600))); 
title('滤波后的信号频谱')
xlabel('频率');ylabel('幅度')
grid on

案例分析2：已知连续信号为$x_a(t)=cos(100\pi t)+sin(100\pi t)+cos(50\pi t)$，用DFT分析其中$x_a(t)$的频谱结构，选择不同的截取长度$T_p$。观察存在的截断效应，试用加窗的方法减少谱间干扰。

clear;close all
fs=400;T=1/fs;						%采样频率和采样间隔
Tp=0.04;N=Tp*fs;					%采样点数N
N1=[N,4*N,8*N];					%设定三种截取长度 
for m=1:3
    n=1:N1(m);
    xn=cos(100*pi*n*T)+ sin(200*pi*n*T)+ cos(50*pi*n*T);
    Xk=fft(xn,4096);
fk=[0:4095]/4096/T;
subplot(3,2,2*m-1);plot(fk,abs(Xk)/max(abs(Xk)));
if m==1 title('矩形窗截取');
end
end
%hamming窗截断
for m=1:3
    n=1:N1(m);
    wn=hamming(N1(m));
    xn=cos(200*pi*n*T)+ sin(100*pi*n*T)+ cos(50*pi*n*T).*wn';
    Xk=fft(xn,4096);
fk=[0:4095]/4096/T;
subplot(3,2,2*m);plot(fk,abs(Xk)/max(abs(Xk)));
if m==1 title('hamming窗截取');
end
end

比较矩形窗与汉明窗的普分析结果可见，矩形窗比用汉明窗分辨率高（泄露小），但是谱间干扰大。因此汉明窗是以牺牲分辨率来换取谱间干扰降低。

布莱克曼FIR滤波器设计:

布莱克曼（Blackman window）的窗函数为：
$$w_{Bl}(n)=[0.42-0.5\ast cos(\frac{2\pi \ast n}{N-1})+0.08\ast cos(\frac{4\pi \ast n}{N-1})]\ast R_N(n)$$
幅值函数为：
$$W_{Hn}(w)=0.42W_R(w)+0.25[W_R(w-\frac{2\pi }{N-1})+W_R(w+\frac{2\pi }{N-1})]$$

$$+0.04[W_R(w-\frac{4\pi }{N-1}+W_R(w+\frac{4\pi }{N-1})]$$

布莱克曼窗幅度函数由5部分相加而成，5部分相加的结果使得旁瓣得到进一步抵消，阻带衰减加大而过渡带加大到$12\pi /N$。

在Matlab中，实现布莱克曼窗的函数为blackman ,其调用格式如下：

w=blackman(N);

:案例：用窗函数法设计数字带通滤波器。下阻带边缘：$W_{s1}=0.3\pi ，A_s=65dB$，下通带边缘:$W_{p1}=0.4\pi ，R_p=1dB$,上通带边缘：$W_{p2}=0.6\pi ，R_p=1dB$，上阻带边缘:$W_{s2}=0.7\pi ，R_p=65dB$。根据窗函数最小阻带衰减的特性，以及参照窗函数的基本参数表，选择布莱克曼窗可以达到75dB的最小阻带衰减，其过渡带为$11\pi /N$。

clear all;
wp1=0.4*pi;
wp2=0.6*pi;
ws1=0.3*pi;
ws2=0.7*pi;
As=65;
tr_width=min((wp1-ws1),(ws2-wp2)); 					%过渡带宽度 
M=ceil(11*pi/tr_width)+1							%滤波器长度
n=[0:1:M-1];
wc1=(ws1+wp1)/2;									%理想带通滤波器的下截止频率
wc2=(ws2+wp2)/2;									%理想带通滤波器的上截止频率
hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);
w_bla=(blackman(M))';								%布莱克曼窗
h=hd.*w_bla;
%截取得到实际的单位脉冲响应
[db,mag,pha,grd,w]=freqz_m(h,[1]);
%计算实际滤波器的幅度响应
delta_w=2*pi/1000;
Rp=-min(db(wp1/delta_w+1:1:wp2/delta_w))
%实际通带纹波
As=-round(max(db(ws2/delta_w+1:1:501)))
As=75
subplot(2,2,1);
stem(n,hd);
title('理想单位脉冲响应hd(n)')
axis([0 M-1 -0.4 0.5]);
xlabel('n');
ylabel('hd(n)')
grid on;
subplot(2,2,2);
stem(n,w_bla);
title('布莱克曼窗w(n)')
axis([0 M-1 0 1.1]);
xlabel('n');
ylabel('w(n)')
grid on;
subplot(2,2,3);
stem(n,h);
title('实际单位脉冲响应hd(n)')
axis([0 M-1 -0.4 0.5]);
xlabel('n');
ylabel('h(n)')
grid on;
subplot(2,2,4);
plot(w/pi,db);
axis([0 1 -150 10]);
title('幅度响应(dB)');
grid on;
xlabel('频率单位:pi');
ylabel('分贝数')

%调用小程序设计1：
function [db,mag,pha,grd,w] = freqz_m(b,a);
% Modified version of freqz subroutine
% ------------------------------------
% [db,mag,pha,grd,w] = freqz_m(b,a);
%  db = Relative magnitude in dB computed over 0 to pi radians
% mag = absolute magnitude computed over 0 to pi radians 
% pha = Phase response in radians over 0 to pi radians
% grd = Group delay over 0 to pi radians
%   w = 501 frequency samples between 0 to pi radians
%   b = numerator polynomial of H(z)   (for FIR: b=h)
%   a = denominator polynomial of H(z) (for FIR: a=[1])
%
[H,w] = freqz(b,a,1000,'whole');
    H = (H(1:1:501))'; w = (w(1:1:501))';
  mag = abs(H);
   db = 20*log10((mag+eps)/max(mag));
  pha = angle(H);
%  pha = unwrap(angle(H));
  grd = grpdelay(b,a,w);
%  grd = diff(pha);
%  grd = [grd(1) grd];
%  grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];
%  grd = median(grd)*500/pi;

%调用小程序设计2：
function hd=ideal_lp(wc,M);
%计算理想低通滤波器的脉冲响应
%[hd]=ideal_lp(wc,M)
%hd=理想脉冲响应0到M-1
%wc=截止频率
% M=理想滤波器的长度
alpha=(M-1)/2;
n=[0:1:(M-1)];
m=n-alpha+eps;
%加上一个很小的值eps避免除以0的错误情况出现
hd=sin(wc*m)./(pi*m);

凯塞窗FIR滤波器设计:

凯塞-贝塞窗（Kaiser-Basel window)的窗函数为：
$$w_k(n)=\frac{I_0(\beta )}{I_0(\alpha )}\ast R_N(n)$$
式中，$\beta =\alpha \sqrt{(1-(\frac{2n}{N-1})-1{)}^2}$。 $I_0(x)$是零阶第一类修正贝塞函数，可用下面级数计算：
$$I_0(x)=1+\sum _{k=1}^{+\infty }（\frac{1}{k!}(\frac{x}{2}{)}^k{）}^2$$
​ $I_0(x)$q取15-25项就可以满足精度要求。通常$\alpha$用以控制窗的形状，$\alpha$加大，主瓣加宽，旁瓣减小，典型数据$4<\alpha <9$。当$\alpha =5.44$s时，窗函数接近汉明窗；当$\alpha =7.865$s时，窗函数接近于布莱克曼窗。其幅值函数为：
$$W_k(w)=w_k(0)+2\sum _{n=1}^{\frac{(N-1)}{2}}(w_k(n)cos(wn))$$
在Matlab中，实现汉宁窗的函数为kaiser,其调用格式如下：

w=kaiser(N);

在Matlab中设计标准响应FIR滤波器可使用fir1函数。fir1函数以经典方法实现加窗性相位FIR滤波器的设计，它可以设计出标准的低通、高通、带通、带阻滤波器。fir1函数用法为：

b = fir1(n,Wn,‘ftype’,wimdow)

各个参数的含义如下：

• b -滤波器系数，n-滤波器阶数。
• Wn -截止频率，$0<=W_n<=1$,$W_n=1$对应于采样频率的一半。当设计带通和带阻滤波器时，$W_n=[W_1,W_2],W_1<w<W_2$
• ftype -当指定ftype时，可设计高通和带阻滤波器。ftype=hight时，设计高通FIR滤波器；ftype=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数。
• window–窗函数。窗函数的长度应等于FIR滤波器系数的个数，即阶数n+1。

案例分析：利用凯塞窗函数设计一个带通滤波器，上截止频率2500Hz,下截止频率1000Hz，过渡带宽200Hz,通带纹波允许差0.1，带阻纹波不大于允差0.02dB，通带幅值为1。

Fs=8000;N=216;
fcuts=[1000 1200 2300 2500];
mags=[0 1 0];
devs=[0.02 0.1 0.02];
[n,Wn,beta,ftype]=kaiserord(fcuts,mags,devs,Fs);
n=n+rem(n,2);
hh=fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
[H,f]=freqz(hh,1,N,Fs);
plot(f,abs(H));
xlabel('频率 (Hz)');
ylabel('幅值|H(f)|');
grid on;

窗函数设计法：

根据前面几节的分析：设计一个FIR低通滤波器通常按照下面的步骤执行：

1. 根据滤波器设计要求，确定滤波器的过渡带宽和阻带衰减要求，选择合适窗函数的类型并进行估计窗函数的宽度N。
2. 根据所求的理想滤波器求出单位脉冲响应$h_d(n)$
3. 根据求得的h(n)求出其频率响应。
4. 根据频率响应验证是否满足技术指标。
5. 若不满足指标要求，则应调整窗函数类型或者长度，然后重复(1),(2),(3)(4)步，直到满足要求为止。

注意：matlab中数据通常是以列向量形式存储的，所以两个向量相乘必须进行转置。计算滤波器的单位脉冲响应h(n),根据窗函数设计理论$h(n)=h_d(n)\ast w(n)$,在matlab中用语句hn=hd*wd实现h(n).

窗函数设计法程序设计如下：

function [h]=usefir1(mode,n,fp,fs,window,r,sample)
% mode:模式(1--高通; 2--低通; 3--带通; 4--带阻)
% n:阶数, 加窗的点数为阶数加1
% fp:高通和低通时指示截止频率, 带通和带阻时指示下限频率
% fs:带通和带阻时指示上限频率
% window:加窗(1--矩形窗; 2--三角窗; 3--巴特窗; 4--汉明窗; 
%5--汉宁窗; 6--布莱克曼窗; 7--凯泽窗; 8--契比雪夫窗)
% r代表加chebyshev窗的r值和加kaiser窗时的beta值
% sample:采样率
% h:返回设计好的FIR滤波器系数
if window==1 w=boxcar(n+1);
end
if window==2 w=triang(n+1);end
if window==3 w=bartlett(n+1);end
if window==4 w=hamming(n+1);end
if window==5 w=hanning(n+1);end
if window==6 w=blackman(n+1);end
if window==7 w=kaiser(n+1,r);end
if window==8 w=chebwin(n+1,r);
end
wp=2*fp/sample;
ws=2*fs/sample;
if mode==1 h=fir1(n,wp,'high',w);
end
if mode==2 h=fir1(n,wp,'low',w);
end
if mode==3 h=fir1(n,[wp,ws],w);
end
if mode==4 h=fir1(n,[wp,ws],'stop',w);
end
m=0:n;
subplot(131);
plot(m,h);grid on;	
title('冲激响应');
axis([0 n 1.1*min(h) 1.1*max(h)]);
ylabel('h(n)');xlabel('n');
freq_response=freqz(h,1);
magnitude=20*log10(abs(freq_response));
m=0:511; f=m*sample/(2*511);
subplot(132);
plot(f,magnitude);grid on;
title('幅频特性');
axis([0 sample/2 1.1*min(magnitude) 1.1*max(magnitude)]);
ylabel('f幅值');xlabel('频率');
phase=angle(freq_response);
subplot(133);plot(f,phase);grid on;
title('相频特性');
axis([0 sample/2 1.1*min(phase) 1.1*max(phase)]);
ylabel('相位');xlabel('频率');

案例分析：假设需要设计一个40阶的带通FIR滤波器，采用汉明窗，采样频率为10kHz,两个截止频率分别为2kHz和3kHz,则需要在Matlab的命令行窗口输入：

h=usefir1(3,60,2000,3000,4,2,10000);