%《MATLAB在電子信息課程中的應用》第六章例6.10程序q610
% 調幅信號通過窄帶濾波器
% 電子工業出版社出版  陳懷琛 吳大正 高西全合著 2001年10月

clear
N=2^16;tf=2*pi;
t=linspace(0,tf,N);dt=tf/N;		% 信號周期為2*pi,分成1000份
w1=[99,100,101];			% 輸入信號的三個頻率分量
U=[0.5,1,0.5];			% 三個頻率分量對應的向量（虛部為零）
b=[2,0];a=[1,2,10001];	% 濾波器分子分母系數向量
u1=U*cos(w1'*t+angle(U')*ones(1,N));	% 輸入信號的時間曲線
U1=fft([u1,u1],2*N)*dt;
w=[0:2*N-1]*2*pi/2/N/dt;
H=freqs(b,a,w);	% 求濾波器在三個頻點上的頻率響應，也可用H=freqs(b,a,w);
% 畫出濾波器的頻率響應曲線，只用三個頻點，圖形不好看，讀者可修改程序得到其完整的頻率響應
figure(1)
subplot(2,1,1),plot(w,abs(H)),grid on		% 幅度
v1=axis;axis([0 400 v1(3) v1(4)])
subplot(2,1,2),plot(w,angle(H)),grid on		% 相位
v2=axis;axis([0 400 v2(3) v2(4)])
U2r=U1(1:N).*H(1:N);
U2=[U2r,fliplr(conj(U2r))];
u2=ifft(U2,2*N)/dt;								% 求和
% 巧妙地利用元素群運算和矩陣運算相結合可把四條語句合成一條語句如下
% u2=abs(U.*H)*cos(w'*t+angle((U.*H).')*ones(1,1001));
% 注意對復數矩陣(U.*H)，(U.*H)'為其共軛轉置，(U.*H).'為轉置而不共軛
figure(2)							% 畫出原信號和濾波后信號的波形作比較
subplot(2,1,1),plot(t,u1)
subplot(2,1,2),plot(t,u2(1:N))