%....... FIR數字濾波器設計..........%
clc,clear
% wp=0.2*pi;Rp=0.25dB; ws=0.3*pi;As=50dB;
% 查第237頁表7.1可知用Hamming窗、Blackman窗均可(最小阻帶衰減>=50dB)，
% 但Hamming窗具有較小的過渡帶(6.6pi/M)，故選擇Hamming窗。
wp=0.2*pi;ws=0.3*pi;
tr_width=ws-wp;
M=ceil(6.6*pi/tr_width)+1   %利用Hamming窗的過渡帶=(6.6pi/M),求窗寬M
n=0:1:M-1;
wc=(ws+wp)/2;       %求通帶截止頻率wc
hd=ideal_lp(wc,M);  %求理想低通沖激響應hd(n)
w_ham=(hamming(M))';
h=hd.*w_ham;        %求FIR濾波器低通沖激響應h(n)
[H,w]=freqz(h,1,1000,'whole');
db=20*log10((abs(H)+eps)/max(abs(H)));  %求FIR濾波器頻響的dB值
delta_w=2*pi/1000;                      %將2pi等分1000份
Rp1=-min(db(1:1:(wp/delta_w+1)))        %求Passband Ripple
As1=-max(db(ws/delta_w+1:1:501))        %求Min Stop attenuation
%plots
figure(3);clf;
subplot(131); stem(n,hd); title('Ideal Impulse Reponse  hd(n)');
axis([0,M-1,-0.1,0.3]);  ylabel('hd(n)');
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','m', 'LineWidth',2);
subplot(132); stem(n,w_ham); title('Hamming Window  w(n)');
axis([0,M-1,0,1.3]);  ylabel('w(n)');
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','r', 'LineWidth',2);
subplot(133); stem(n,h); axis([0,M-1,-0.1,0.3]);
ylabel('h(n)');title('Actual Impulse  Response  h(n)');
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','b', 'LineWidth',2);
figure(1);clf;
subplot(221); stem(n,hd); title('Ideal Impulse Reponse');
axis([0,M-1,-0.12,0.3]);  ylabel('hd(n)');
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','m', 'LineWidth',1);
subplot(222); stem(n,w_ham); title('Hamming Window');
axis([0,M-1,0,1.3]);  ylabel('w(n)');
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','g', 'LineWidth',1);
subplot(223); stem(n,h); axis([0,M-1,-0.1,0.3]);
ylabel('h(n)');title('Actual Impulse  Response');
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','r', 'LineWidth',1);
subplot(224); plot(w/pi,db);
axis([0,1,-100,10]);  ylabel('Decibels');title('Magnitude  Response in dB');
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);
set(gca,'YTickMode','manual','YTick',[-50 -0.25 0]);grid;
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','k', 'LineWidth',1);
figure(4)
subplot(121); plot(w/pi,db);
axis([0,1,-100,10]);  ylabel('Magnitude');title('Magnitude  Response in dB');
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);
set(gca,'YTickMode','manual','YTick',[-50  0]);grid;
subplot(122); plot(w/pi,angle(H));
axis([0,1,-3.5,3.5]);  ylabel('angle');title('Phase  Response in rad');
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','b', 'LineWidth',2);
% 用頻率抽樣法設計FIR濾波器,過渡帶內一個樣本T1，N=40。
% wp=0.2pi,ws=0.3pi,Rp=0.25dB,As=50dB
T1=0.35;
N=40,alpha=(N-1)/2;   %N為偶數
l=[0:1:N-1]; wl=(2*pi/N)*l;
Hrs=[ones(1,5),T1,zeros(1,29),T1,ones(1,4)]; %偶對稱
Hdr=[1 1 0 0 ];wdl=[0 0.2 0.3 1];  k1=0:(N/2-1);  k2=(N/2+1):N-1;%依據公式7-107
angH=[-alpha*(2*pi)/N*k1,0,alpha*(2*pi)/N*(N-k2)];  H=Hrs.*exp(j*angH);
h=real(ifft(H,N));[H,w]=freqz(h,1,1000,'whole');
%[db mag pha grd w]=freqz_m(h,1);
db=20*log10((abs(H)+eps)/max(abs(H)));%求FIR濾波器頻響的dB值
delta_w=2*pi/1000;                    %將2pi等分1000份
Rp2=-min(db(1:1:(wp/delta_w+1)))      %求Passband Ripple
As2=-max(db(ws/delta_w+1:1:501))      %求Min Stop attenuation
%plot
figure(5);clf;
subplot(221);plot(wl(1:21)/pi,Hrs(1:21),'o',wdl,Hdr,'linewidth',2);
title('Ideal Amplitude Response');
axis([0 1 -0.1 ,1.2]);  ylabel('Hr(k)');
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);
set(gca,'YTickMode','manual','YTick',[0 T1 1]); grid;
subplot(222); stem(l,h,'m');
title('Impulse Response');axis([-1,N,-0.1,0.3]);ylabel('h(n)');
subplot(223); plot(w/pi,abs(H),wl(1:31)/pi,Hrs(1:31),'o','linewidth',2);axis([0 1 -0.2,1.2]);
title('Amplitude Response');xlabel('frequency in pi units');ylabel('Hr(w)');grid;
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);
set(gca,'YTickMode','manual','YTick',[0 T1 1]);
subplot(224); plot(w/pi,db,'r');axis([0,1,-100,10]);grid;  ylabel('Decibels');
title('Magnitude Response');xlabel('frequency in pi units');
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);
set(gca,'YTickMode','manual','YTick',[-50,0])