实验一：低通采样定理和内插与抽取实现

1、实验一：低通采样定理和内插与抽取实现一 实验目的1. 连续信号和系统的表示方法，以及坊真方法。2. 用 MATLAB 实现连续信号采用与重构的方法，3. 采样信号的插值和抽取等重采样实现方法。4. 用时域采样信号重构连续时域信号的原理和方法。5. 用 MATLAB 绘图函数表示信号的基本方法，实验数据的可视化表示。二 原理1 、时域抽样定理令连续信号 xa(t)的傅里叶变换为 Xa(j ), 抽样脉冲序列 p(t)傅里叶变换为 P(j ),抽样后的信号 x(t)的傅里叶变换为 X(j )若采用均匀抽样 , 抽样周期 Ts, 抽样频率为 s=2 fs, 由前面分析可知 :抽样的过程可以通过抽样脉

2、冲序列 p(t) 与连续信号 xa(t)相乘来完成 , 即满足 : x(t)=xa(t) p(t) ，又周期信号 f(t) 傅里叶变换 为：故可以推得 p(t)的傅里叶变换为 :其中:根据卷积定理可知 :得到抽样信号 x(t) 的傅里叶变换为 :其表明 :信号在时域被抽样后 ,他的频谱 X(j ) 是连续信号频谱 X(j )的形状 以抽样频率 为间隔周期重复而得到 ,在重复过程中幅度被 p(t)的傅里叶级数 Pn 加权。因为Pn只是 n的函数,所以 X(j )在重复的过程中不会使其形状发生变化。假定信号 x(t)的频谱限制在 - m+ m 的范围内, 若以间隔 Ts对 xa(t)进行 抽样,可

3、知抽样信号 X(t)的频谱 X(j )是以 s 为周期重复。显然 ,若在抽样的 过程中 s<2 m,则 X(j ) 将发生频谱混叠现象 ,只有在抽样的过程中满足 s>=2 m 条件,X(j )才不会产生频谱的混叠 ,接收端完全可以由 x(t) 恢复原 连续信号 xa(t),这就是低通信号抽样定理的核心内容。2、信号的重建从频域看 ,设信号最高频率不超过折叠频率Xa(j )=Xa(j ) | |< s/2Xa(j )=0 | |> s/2 则理想取样后的频谱就不会产生混叠 ,故有 :让取样信号 x(t) 通过一带宽等于折叠频率的理想低通滤波器H(j )=T | |<

4、 s/2H(j )=0 | |> s/2 滤波器只允许通过基带频谱 ,即原信号频谱 ,故:Y(j )=X(j )H(j )=Xa(j )因此在滤波器的输出得到了恢复的原模拟信号 :y(t)=xa(t)从时域上看 ,上述理想的低通滤波器的脉冲响应为根据卷积公式可求得理想低通滤波器的输出为由上式显然可得则:上式表明只要满足取样频率高于两倍信号最高频率 ,连续时间函数 xa(t)就可 用他的取样值 xa(nT)来表达而不损失任何信息， 这时只要把每一个取样瞬时值与 内插函数式相乘求和即可得出 xa(t),在每一取样点上 ,由于只有该取样值所对应的 内插函数式不为零 ,所以各个取样点上的信号值不

5、变。内容1. 连续时间信号线性滤波实现设计一个 Butterworth 模拟带通滤波器，设计指标为：通带频率： 1000- 2000Hz，两侧过渡带宽 500Hz，通带波纹 1dB，阻带衰减 100dB。假设一个信号， 其中 f1=100Hz,f2=1500Hz,f3=2900Hz 。信号的采样频率为 10000Hz。试将原信号 与通过该滤波器的模拟信号进行比较。参考程序如下：wp=1000 2000*2*pi;ws=500 2500*2*pi;Rp=1;Rs=100; %滤波器设计参数，对于给定 Hz 应乘以 2N,Wn=buttord(wp,ws,Rp,Rs,'s'); %

6、求得滤波器的最小阶数和截止频率w=linspace(1,3000,1000)*2*pi;%设置绘制频率响应的频率点b,a=butter(N,Wn,'s');%设计模拟 Butterworth 滤波器H=freqs(b,a,w);%计算给定频率点的复数频率响应magH=abs(H);phaH=unwrap(angle(H); % 计算幅频响应和相频响应 plot(w/(2*pi),20*log10(magH);%以频率为横坐标绘制幅频响应xlabel(' 频率 /Hz');ylabel(' 振幅 /dB'); title('Butterwo

7、rth 模拟带通滤波器 ');hold on;plot(1000 1000,ylim,'r');plot(2000 2000,ylim,'r');% 绘带边界 grid on figure(2) dt=1/10000;%模拟信号采样间隔f1=100;f2=1500;f3=2900;% 输入信号的三个频率成分 t=0:dt:0.04;%给定模拟时间段x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); % 输入信号 H=tf(b,a); %滤波器在 MATLAB系统中的表示 y,t1=lsim(H,

8、x,t);%模拟输出subplot(2,1,1),plot(t,x),title(' 输入信号 ') % 绘出输入信号 subplot(2,1,2),plot(t1,y)%绘制输出信号title(' 输出信号 '),xlabel(' 时间 /s')2. 理想矩形滤波器的时域表示clear; clc; fh=100; for I=1:4 k=I; fs=k*2*fh; N=10*k; n=-N:N; dt=1/fs; T=N*dt; t=-T:dt:T;%h=2*fh/fs*sinc(n/k);%h=sinc(n/k); subplot(2,2,I

9、) plot(t,h); hold on ； stem(t,h); hold on; plot(t,zeros(length(t),'linewidth',3); title('fs/2fh=k,k=',num2str(k),'fontsize',28); axis('off') endfs/2fh=k,k=2fs/2fh=k,k=4fs/2fh=k,k=1 fs/2fh=k,k=33. 连续时间信号的采样和重建 、分别用 150HZ及 300HZ对信号采样源信号为：fa=5*sin(2*pi*40*t1)+1.8*sin(4*p

10、i*40*t1)+0.8*sin(5*pi*40*t1) ，用 150Hz 的频率 对 f(t) 进行采样，其采样图如图 1 所示；用 300Hz 的频率对 f(t) 进行采样，其采样图如图 2 所示。程序如下： fs1=150;t1=-0.1:1/fs1:0.1; fa=5*sin(2*pi*40*t1)+1.8*sin(4*pi*40*t1)+0.8*sin(5*pi*40*t1); figure(1);plot(t1,fa),xlabel('fs1=150Hz 时 ,fa 采样时域图 '); hold off;fs2=300;t2=-0.1:1/fs2:0.1;fb=5*

11、sin(2*pi*40*t2)+1.8*sin(4*pi*40*t2)+0.8*sin(5*pi*40*t2);时,fb 采样时域图 ');0.1figure(2);plot(t2,fb),xlabel('fs2=300Hz-8-0.1-0.08-0.06-0.04-0.020 0.020.040.060.08fs1=150Hz时 ,fa 采 样 时 域 图8-0.1 -0.08 -0.06 -0.04 -0.020 0.02 0.04 0.06 0.08 0.1fs2=300Hz时 ,fb 采 样 时 域 图图1 150HZ采样频率对信号采样图图 2 300HZ采样频率对信号

12、采样图 、对信号进行快速离散傅里叶变换将两个采样信号进行快速离散傅里叶变换 (FFT), 用 150Hz的频率对 f(t) 进行 采样，其采样后快速傅立叶变换频谱图图 3 所示；用 300Hz 的频率对 f(t) 进行 采样，其采样后快速傅立叶变换频谱图图 4 所示。程序如下：f=40;fs=150;N=300;k=0:N-1;t=-0.1:1/fs:0.1;w1=150*k/N;fa=5*sin(2*pi*f*t)+1.8*sin(4*pi*f*t)+0.8*sin(5*pi*f*t);xfa=fft(fa,N);xf1=abs(xfa);时,fa 经fft 后频谱图 . 单位： Hz

13、9;);figure(1);plot(w1,xf1),xlabel('fs=150Hz f=40;fs=300;N=300;k=0:N-1;t=-0.1:1/fs:0.1;w2=300*k/N fb=5*sin(2*pi*f*t)+1.8*sin(4*pi*f*t)+0.8*sin(5*pi*f*t);xfb=fft(fb,N);xf2=abs(xfb);figure(2);plot(w2,xf2),xlabel('fs=300Hz 时,fb 经fft 后频谱图 .单位： Hz ');80fs=150Hz 时 ,fa 经 fft 后 频 谱 图 . 单 位 ： Hz图

14、3 150HZ采样后经 FFT后频谱图图4 300HZ采样后经 FFT后频谱图 、信号的重建我们可以通过利用内插法把原信号从采样信号中恢复出来，观察信号在满 足怎样的采样条件下能够恢复为原信号，图 5和图 6分别为恢复后的原信号 程序如下：Wm=180*pi;Wc=Wm;fs1=150;Ws=2*pi*fs1; n=-800:800;nTs1=n/fs1; fa=5.1*sin(2*pi*40*nTs1)+1.8*sin(4*pi*40*nTs1)+0.8*sin(5*pi*40*nTs1); Dt=1/fs1;t1=-0.1:Dt:0.1;fa1=fa/fs1*Wc/pi*sinc(Wc/p

15、i)*(ones(length(nTs1),1)*t1-nTs1'*ones(1,length( t1);figure(1);plot(t1,fa1);axis(-0.1 0.1 -8 8); xlabel('fs=150Hz 时 ,fa 利用内插由样本重建原信号图 .');Wm=180*pi;Wc=Wm;fs2=300;Ws=2*pi*fs2; n=-800:800;nTs2=n/fs2;fb=5.1*sin(2*pi*40*nTs2)+1.8*sin(4*pi*40*nTs2)+0.8*sin(5*pi*40*nTs2); Dt=1/fs2;t1=-0.1:Dt:0

16、.1;fb1=fb/fs2*Wc/pi*sinc(Wc/pi)*(ones(length(nTs2),1)*t1-nTs2'*ones(1,length( t1);figure(2);plot(t1,fb1);axis(-0.1 0.1 -8 8); xlabel('fs=300Hz 时 ,fb 利用内插由样本重建原信号图 .'); grid;8图5 150HZ 采样后的信号的重建信号8图6 300HZ 采样后的信号的重建信号4. 采样信号的抽取和插值t = 0:.00025:1;% Time vectorx = sin(2*pi*30*t) + sin(2*pi*60

17、*t); figure(1) subplot(211) stem(t(1:120),x(1:120);hold on y = decimate(x,4);% View the original and decimated signals: stem(x(1:120), axis(0 120 -2 2) % Original signaltitle( 'Original Signal' )subplot(212)stem(y(1:30)title( 'Decimated Signal' )% Decimated signalDecimated Signal2%Fi

18、le_C3:UpSampl.m%该程序仿真通过零阶保持内插对信号进行上抽样，分析上抽样对信号频谱的影响clear allclcf0=0.06;% 信号数字频率N=256;dt=1;% 抽样时间t=0:N-1*dt;sig=sin(2*pi*f0*t);I=4; % 内插因子N1=N*I;Addsig(1:N1)=0;for k=1:Nfor i=1:IdLIAddsig(I*k-I+i)=sig(k); % 对原始数字信号补零内插， 抽样时间变为 endendfigure(1)subplot(2,1,1)t2=0:N1-1;stem(t,sig);title('Original'

19、;)xlabel('Sample time /n')axis(0 19 -1.2 1.2)subplot(2,1,2)stem(Addsig);axis(0 80 -1.2 1.2)xlabel('Sample time /n')title('Signal output of Zero-order hold interpolator')Am=fft(sig);AddAm=fft(Addsig);f1=0:2/N:2-1/N;f2=0:2/N1:2-1/N1;figure(2)subplot(2,1,1) plot(f1,20*log10(abs(

20、Am);axis(0 1 -20 60);ylabel('Amplitude(dB)')xlabel('omega/ pi')title('Amplitude Spectrum')subplot(2,1,2)plot(f2,20*log10(abs(AddAm)axis(0 1 -20 60);ylabel('Amplitude(dB)')xlabel('omega/ pi')title('Amplitude Spectrum of Zero-order hold interpolation signal&

21、#39;)b,a=butter(20,1/I);% 设计截止频率为 pi/l，阶数为 10 的低通巴特沃思滤波器y=filter(b,a,Addsig);% 对补零后的信号进行低通滤波，完成上抽样过程UpAm=fft(y);figure(3) subplot(2,1,1) plot(y);axis(0 1023 -1.2 1.2) title('output signal out of Interpolator') xlabel('Sample time /n') subplot(2,1,2)plot(f2,20*log10(abs(UpAm) axis(0 1

22、 -30 50);ylabel('Amplitude(dB)') xlabel('omega/ pi')%File_C3:UpSampl.m%该程序仿真通过零阶保持内插对信号进行上抽样，分析上抽样对信号频谱的影响 clear allclcf0=0.06; % 信号数字频率N=256;dt=1; % 抽样时间t=0:N-1*dt; sig=sin(2*pi*f0*t);I=4;% 内插因子N1=N*I;Addsig(1:N1)=0;for k=1:Nfor i=1:IdLIAddsig(I*k-I+i)=sig(k); % 对原始数字信号补零内插， 抽样时间变为

23、endendfigure(1) subplot(2,1,1) t2=0:N1-1; stem(t,sig); title('Original') xlabel('Sample time /n') axis(0 19 -1.2 1.2) subplot(2,1,2) stem(Addsig);axis(0 80 -1.2 1.2)xlabel('Sample time /n')title('Signal output of Zero-order hold interpolator') Am=fft(sig);AddAm=fft(Addsig); f1=0:2/N:2-1/N;f2=0:2/N1:2-1/N1;figure(2) subplot(2,1,1) plot(f1,20*log10(abs(Am);axis(0 1 -20 60); ylabel('Amplitude(dB)') xlabel('omega/ pi') title('Amplitude