FIR数字滤波器设计频率抽样法MATLAB仿真

http://blog.csdn.net/shenziheng1/article/details/53366644

https://wenku.baidu.com/view/ee4ec55f27d3240c8447ef73.html    实验报告，案例

http://www.docin.com/p-397194639.html   实验报告，案例

实际上就是由连续谱，连续取样N个点，相当于得到N个点的离散傅里叶变换值。通过反变换得到取real，得到时域信号。

频率抽样法是从频域出发，在频域直接设计，把给定的理想频率响应加以等间隔抽样，并以此作为实际FIR滤波器的频率响应。设所需滤波器的频率响应为。

现要求设计一个M阶的FIR滤波器h[k]，使得在M+1个抽样点上，FIR滤波器的频率响应与所需的频率响应相等，即

由设计的要求给定，h[k]通过设计来确定。如果M+1个方程是线性无关的，则可以通过求解M+1阶的线性方程来得出FIR滤波器的h[k]。对的一些特殊抽样方法，上述方程的解可以直接由IDFT得到。由于要求设计出的滤波器是实系数的线性相位FIR滤波器，所以的抽样值还需要满足线性相位滤波器的约束条件。

 I型和II型线性相位滤波器的，III型和IV型线性相位滤波器的。为了使设计出的滤波器具有线性相位， 在M+1个抽样点上的值应为：

下面分别讨论四种线性相位滤波器在抽样点上的值。

2.Ⅰ型线性相位FIR滤波器

 

 

I型（M为偶数，h[k]偶对称）线性相位FIR滤波器在M+1个抽样点值为：

 

上式表明I型线性相位FIR滤波器 在 的值可由 在 的值确定 。在的值确定后，对做M+1点的IDFT即可得到I型线性相位滤波器的h[k]。

 

 

3.Ⅱ型线性相位FIR滤波器

 

 

II型（M为奇数，h[k]偶对称）线性相位FIR滤波器在M+1个抽样点值为：

上式表明 II型线性相位FIR滤波器在的值可由在的值确定。

 

 

4.Ⅲ型线性相位FIR滤波器

 

 

III型（M为偶数，h[k]奇对称）线性相位FIR滤波器在M+1个抽样点值为：

上式表明 III型滤波器线性相位FIR滤波器在的值可由在的值确定。

 

 

5.Ⅳ型线性相位FIR滤波器

 

 

IV型（M为奇数，h[k]奇对称）线性相位FIR滤波器在M+1个抽样点值为：

上式表明 IV型线性相位FIR滤波器在的值可由在的值确定。

 

6.小结

对进行频率抽样，就是在z平面单位圆上的N个等间隔点上抽样出频率响应值。在单位圆上可以有两种抽样方式，第一种是第一个抽样点在w=0处，第二种是第一个抽样点在w=pi/M处，每种方式可分为M为偶数与M为奇数两种。

为了提高逼近质量，使逼近误差更小，也就是减小在通带边缘由于抽样点的徒然变化而引起的起伏变化（这种起伏振荡使阻带内最小衰减变小，例如从衰减30dB变小为衰减20dB）。和窗口法的平滑截断一样，这里是使理想频率响应的不连续点的边缘加上一些过渡的抽样点（在这些点上抽样的最佳值由计算机算出），从而增加过渡带，减小频带边缘的突变，也就是减小了起伏振荡，增大了阻带最小衰减。这些抽样点上的取值不同，效果也就不同。如果精心设计过渡带的抽样值，就有可能使它的游泳频带的博文减小，从而设计出较好的滤波器。一般过渡带取一、二、三点抽样值即可得到满意结果。

在理想低通滤波器的设计中，若不增加过渡点，阻带和通带之间的衰减约为-21dB，如果在通带和阻带之间增加一个采样点，阻带的最小衰减可以提高到-65dB，如果增加两个采样点，阻带的最小衰减可以提高到-75dB，如果增加3个采样点，阻带的最小衰减可以提高到-85dB至-95dB。
频率抽样法的优点是可以在频域直接设计，并且适合于最优化设计；缺点是抽样频率只能等于2pi/M的整数倍或等于2pi/M的整数倍上加上pi/M，因而不能确保截止频率Wc的自由取值。要想实现自由选择频率，则必须增加抽样点数M，但这种计算量加大。

 

 

M=63; Wp=0.5*pi;%所需频率采样点个数及通带截止频率  
m=0:(M+1)/2; Wm=2*pi*m./(M+1);%2*pi 划分为M+1个点  对称谱；  m指示了单边谱
mtr=floor(Wp*(M+1)/(2*pi))+2;%向负方向入floor(3.5)=3;floor(-3.2)=-4  
Ad=[Wm<=Wp];%截止频率之前都为1  
Ad(mtr)=0.38;%截止点对应的频点下标  
Hd=Ad.*exp(-1i*0.5*M*Wm);%构造频域采样向量H(k)  
Hd=[Hd conj(fliplr(Hd(2:(M+1)/2)))];  
%fliplr函数实现矩阵的左右翻转conj是求复数的共轭  
h=real(ifft(Hd));%h（n）=IDFT[H(k)]  
w=linspace(0,pi,1000);%用于产生0,pi之间的1000点行矢量  
H=freqz(h,[1],w);%滤波器的幅频特性图  
figure(1)  
plot(w/pi,20*log10(abs(H)));%参数分别是归一化频率与幅值  
xlabel('归一化角频率');ylabel('增益/分贝');title('滤波器的增益响应');  
axis([0 1 -50 0.5]);  
f1=100;f2=300;f3=700;fs=2000;%待滤波正弦信号频率及采样频率  
figure(2)  
subplot(211)  
t=0:1/fs:0.25;%定义时间范围和步长  
s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);%滤波前信号  
plot(t,s);%滤波前的信号图像  
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');  
subplot(212)  
Fs=fft(s,512); AFs=abs(Fs);%将信号变换到频域及信号频域图的幅值  
f=(0:255)*fs/512;%频率采样  
plot(f,AFs(1:256));%滤波前的信号频域图  
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');  
figure(3)  
sf=filter(h,1,s);%使用filter函数对信号进行滤波  
subplot(211)  
plot(t,sf)%滤波后的信号图像  
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');  
axis([0.2 0.25 -2 2]);%限定图像坐标范围  
subplot(212)  
Fsf=fft(sf,512); AFsf=abs(Fsf);%滤波后的信号频域图及信号频域图的幅值  
f=(0:255)*fs/512;%频率采样  
plot(f,AFsf(1:256))%滤波后的信号频域图  
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');