%《MATLAB在電子信息課程中的應用》第七章例7.13程序q713
% IIR濾波器實現及FT的時域卷積定理驗證
% 電子工業出版社出版  陳懷琛 吳大正 高西全合著 2001年10月

clear;close all
%
%產生輸入信號x(n)
n=0:255;N=4096;
x=cos(0.04*pi*n)+cos(0.08*pi*n)+cos(0.4*pi*n);
w=randn(size(x));%產生正態零均值噪聲
x=x+0.3*w;
%
%求H(z)分子分母多項式系數向量B和A
b=[1,2,1];  %(1+z-1)2 的展開系數
B=0.0003738*conv(conv(b,b),b);%嵌套調用卷積函數conv
a1=[1,-1.2686,0.7051];
a2=[1,-1.0106,0.3583];
a3=[1,-0.9044,0.2155];
A=conv(conv(a1,a2),a3);
%
%對x(n)濾波
y=filter(B,A,x);
%繪圖
X=fft(x,N);
Y=fft(y,N);
subplot(3,2,1);stem(x,'.')
axis([0,max(n)/4,min(x),max(x)]);line([0,max(n)],[0,0])
title('輸入信號x(n)');xlabel('n');ylabel('x(n)')
subplot(3,2,5);stem(y,'.')
axis([0,max(n)/4,min(y),max(y)]);line([0,max(n)],[0,0])
title('輸出信號 y(n)');xlabel('n');ylabel('y(n)')
k=0:N-1;f=2*k/N;
subplot(3,2,2);plot(f,abs(X))
title('輸入信號x(n)的幅頻曲線');xlabel('ω/π');ylabel('|FT[x(n)]|')
axis([0,0.5,0,max(abs(X))]);
subplot(3,2,6);plot(f,abs(Y))
title('輸出信號y(n)的幅頻曲線');xlabel('ω/π');ylabel('|FT[y(n)]|')
axis([0,0.5,0,max(abs(Y))]);
[h,f]=freqz(B,A,N,'whole');
figure(2)
subplot(3,2,1);plot(f/pi,abs(h))
title('濾波器幅頻響應曲線');xlabel('ω/π');ylabel('H幅度')
axis([0,0.5,0,max(abs(h))]);
Ym=h'.*X;
subplot(3,2,2);plot(f/pi,abs(Ym))
title('|FT[x(n)]FT[h(n)]|');xlabel('ω/π');ylabel('Ym幅度')
axis([0,0.5,0,max(abs(Ym))]);