samp10_2.m

MATLAB7.x数字信号处理 光盘内容

%Samp10_2
Fs=100;dt=1/Fs;   %给定模拟输入信号采样间隔
f1=6;f2=10;f3=19;   %模拟输入信号的频率成分
N=500;             %数据点数
t=[0:N-1]*dt;   %模拟输入信号的时间序列
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); %模拟输入信号
X=fft(x);   %得到输入信号的Fouier变换
%设计宽带Butterworth模拟滤波器
wp=[1 40]*2*pi;ws=[0.2 42.5]*2*pi; %通带和阻带截止频率
Rp=1;Rs=20;   %通带波纹和阻带衰减
[Order,Wn]=buttord(wp,ws,Rp,Rs,'s');   %求解Butterworth滤波器的最小阶数
w=linspace(0,Fs/2,N/2)*2*pi;   %设置绘制幅频特性的频率
[b1,a1]=butter(Order,Wn,'s');     %设计Butterworth带通滤波器
H=freqs(b1,a1,w);   %计算带通滤波器的频率响应
magH=abs(H);phaH=unwrap(angle(H));  %求得振幅谱和相位谱
figure(1)
subplot(3,1,1),plot(w/(2*pi),20*log10(magH));  %绘制幅频特性
xlabel('频率/Hz')
ylabel('振幅/dB');
ylim([-100 0]);xlim([0,50]); grid on
title('宽带模拟带通滤波器');

%模拟输出
H=freqs(b1,a1,w);   %滤波器的幅频响应
H1=zeros(1,N);
for ii=1:N/2   %构建与Fourier变换得到前后数据共轭的形式
    H1(ii)=H(ii);    
    H1(N-ii)=conj(H1(ii));
end
X=fft(x);   %将原信号进行Fourier变换
Y=zeros(1,N);
for ii=1:N
    Y(ii)=X(ii).*H1(ii);   %按（10-4）式得到的输出
end
y=real(ifft(Y));     %将频率域的数据通过Fourier变换转换到时间域
subplot(3,1,2),plot(t,x),title('输入信号')  %绘制输入信号
subplot(3,1,3),plot(t,y)  %绘制输出信号
title('宽带滤波器的模拟输出信号') 
xlabel('时间/s')

%恢复地面运动
XX=zeros(1,N);
for ii=1:N
    if(abs(H1(ii))>1.0e-1)   %为防止位于阻带内数据不稳定
        XX(ii)=Y(ii)./H1(ii);   %恢复地面运动公式（10-7）
     end
end
xx=real(ifft(XX));  %转换到时间域
figure(2)
plot(t,x,t,xx,':r')
legend('原始信号','恢复信号',1)  %给出图例
title('原始信号和恢复信号')
xlabel('时间/s')
%窄带滤波器的设计
wp=[5 7]*2*pi;ws=[3.5 7.5]*2*pi;  %窄带滤波器的通带和阻带边界频率
Rp=1;Rs=20;   %窄带滤波器的通带波纹和阻带衰减
[Order,Wn]=buttord(wp,ws,Rp,Rs,'s'); %计算窄带滤波器的最小阶数和截止频率  
w=linspace(0,Fs/2,N/2)*2*pi;   %设置绘制窄带滤波器幅频特性的频率
[b2,a2]=butter(Order,Wn,'s');   %设计Butterworth窄带滤波器
H=freqs(b2,a2,w);    %计算Butterworth窄带滤波器的频率响应
magH=abs(H);phaH=unwrap(angle(H));  %计算振幅谱和相位谱
figure(3)
subplot(3,1,1),plot(w/(2*pi),20*log10(magH));  %绘制幅频特性
xlabel('频率/Hz'),ylabel('振幅/dB'); ylim([-100 0]);xlim([0,50])
title('窄带模拟滤波器');
grid on
%模拟输出
H=freqs(b2,a2,w);
H2=zeros(1,N);
for ii=1:N/2
    H2(ii)=H(ii);
    H2(N-ii)=conj(H2(ii));
end
X=fft(x);
Y1=zeros(1,N);
for ii=1:N
    Y1(ii)=X(ii).*H2(ii);
end
y1=real(ifft(Y1));
subplot(3,1,2),plot(t,x),title('输入信号')  %绘制输入信号
subplot(3,1,3),plot(t,y1)  %绘制输出信号
title('窄带滤波器的模拟输出') 
xlabel('时间/s')
%运用宽频带仪器的输出仿真得到窄带仪器上的输出
figure(4)
XY=zeros(1,N);
for ii=1:N
    if(abs(H1(ii))>1.0e-2)
        XY(ii)=Y(ii)*H2(ii)/H1(ii);  %公式（10-6）
    end
end
xy=real(ifft(XY));   %转换到时间域
plot(t,y1,t,xy,':r')
legend('窄带输出','仿真输出',1)
xlabel('时间/s')