q610a.m

<matlab在电子信息中的应用>一书中个章节的配套源代码。

%《MATLAB在电子信息课程中的应用》第六章例6.10程序q610
% 调幅信号通过窄带滤波器:用fft及ifft函数来求
% 电子工业出版社出版  陈怀琛 吴大正 高西全合著 2001年10月

clear
N=2^10;tf=2*pi;
t=linspace(0,tf,N+1);
dt=tf/N;t=t(1:N);		% 信号周期为2*pi,分成N份
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,N)*dt;
w=[0:N-1]*2*pi/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/2).*H(1:N/2);
%U2=[U2r,fliplr(conj(U2r))];
U2=[U2r,0,fliplr(conj(U2r(2:N/2)))];    % 为了构成共扼对称的U2,这一句很关键，要仔细读
u2=ifft(U2,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)