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('信號濾波后頻域圖');