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

实验一：低通采样定理和内插与抽取实现实验目的连续信号和系统的表示方法，以及坊真方法。用MATLAB实现连续信号采用与重构的方法，采样信号的插值和抽取等重采样实现方法。用时域采样信号重构连续时域信号的原理和方法。用MATLAB绘图函数表示信号的基本方法，实验数据的可视化表示。原理1、时域抽样定理令连续信号xa(t)的傅里叶变换为Xa(jQ)，抽样脉冲序列p(t)傅里叶变换为P(jQ),抽样后的信号xA(t)的傅里叶变换为X气Q)若采用均匀抽样，抽样周期Ts,抽样频率为Qs=2兀fs,由前面分析可知：抽样的过程可以通过抽样脉冲序列p(t)与连续信号xa(t)相乘来完成，即满足:xNt)=xa(t)p(t)，又周期信号f(t)傅里叶变换为：=2兀E乩O'Q—E]故可以推得p(t)的傅里'叶变换为：00代」、..矿"*3/:1已二户勺出根据卷积定理可知■:-X(jQ)=嘉XJj切*P(jQ)得到抽样信号x(t)的傅里叶变换为：X"二Q-E]月=—oa其表明:信号在时域被抽样后，他的频谱X(jQ)是连续信号频谱X(jQ)的形状以抽样频率Q为间隔周期重复而得到,在重复过程中幅度被p(t)的傅里叶级数Pn加权。因为Pn只是n的函数，所以X(jQ)在重复的过程中不会使其形状发生变化。假定信号x(t)的频谱限制在-Qm〜+Qm的范围内，若以间隔Ts对xa(t)进行抽样,可知抽样信号XA(t)的频谱X气Q)是以Qs为周期重复。显然，若在抽样的过程中Qs<2Qm,则X气Q)将发生频谱混叠现象，只有在抽样的过程中满足Qs>=2Qm条件,X^CjQ)才不会产生频谱的混叠，接收端完全可以由xA(t)恢复原连续信号xa(t)，这就是低通信号抽样定理的核心内容。2、信号的重建从频域看,设信号最高频率不超过折叠频率：Xa(jQ)=Xa(jQ)IQl<Qs/2Xa(jQ)=0IQl>Qs/2则理想取样后的频谱就不会产生混叠,故有：00EJ印二+j◎一皿OJ让取样信号约(t)通过一带宽I等于折叠频率的理想低通滤波器：H(jQ)=TIQl<Qs/2H(jQ)=0IQl>Qs/2滤波器只允许通过基带频谱，即原信号频谱，故：Y(jQ)=X气Q)H(jQ)=Xa(jQ)因此在滤波器的输出得到了恢复的原模拟信号：y(t)=xa(t)从时域上看，上述理想的低通滤波器的脉冲响应为：.兀]s引11时ha=T-fH(jC|L\1。=—-—根据卷积公式可求得理想低通滤波器的输出为：卞001)-£W)hfl(I-nT)R=—8由上式显然可得：ha(t-nT)=sin(^/T)(t-nF)/(^/T)(t-nT)则:™silll-rT})(时=仃)兀=nf"上式表明只要满足取样频率高于两倍信号最高频率7,连续时间函数xa(t)就可用他的取样值xa(nT)来表达而不损失任何信息，这时只要把每一个取样瞬时值与内插函数式相乘求和即可得出xa(t)，在每一取样点上,由于只有该取样值所对应的内插函数式不为零，所以各个取样点上的信号值不变。内容设计连续时间信号线性滤波器分离信号组份已知信号x(t)=^99(m+1)cos(2兀(100m+50)t),试设计滤波器，分离出如下信号:m=0m=1,2,3……50m=51,52,53……100m=40,41,42..….60m=1,2,3......40,61,62......100据以下采样频率：(a)f=20000Hz(b)f=10000Hz(c)f=30000Hz求信号频谱及相应的滤波器。参考程序如下：设计一个Butterworth模拟带通滤波器，设计指标为：通带频率：1000-2000Hz，两侧过渡带宽500Hz，通带波纹1dB,阻带衰减100dB。假设一个信号，其中f1=100Hz,f2=1500Hz,f3=2900Hz。信号的采样频率为10000Hz。试将原信号与通过该滤波器的模拟信号进行比较。参考程序如下：wp=[10002000]*2*pi;ws=[5002500]*2*pi;Rp=1;Rs=100;%滤波器设计参数，对于给定Hz应乘以2[N,Wn]=buttord(wp,ws,Rp,Rs,'s');%求得滤波器的最小阶数和截止频率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('Butterworth模拟带通滤波器')；holdon;plot([10001000],ylim,'r');plot([20002000],ylim,'r');%绘带边界gridonfigure(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,x,t);%模拟输出subplot(2,1,1),plot(t,x),title('输入信号')%绘出输入信号subplot(2,1,2),plot(t1,y)%绘制输出信号title('输出信号')，xlabel('时间/s')连续时间信号的采样和重建已知信号X(t)=磬(m+1)cos(2兀(100m+50)t)，试以f广20000Hz采样频率对该m=0信号采样，并用插值公式重建该信号。参考程序：、分别用150HZ及300HZ对信号采样源信号为：fa=5*sin(2*pi*40*t1)+1.8*sin(4*pi*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采样时域图');holdoff;fs2=300;t2=-0.1:1/fs2:0.1;fb=5*sin(2*pi*40*t2)+1.8*sin(4*pi*40*t2)+0.8*sin(5*pi*40*t2);figure(2);plot(t2,fb),xlabel('fs2=300Hz时，fb采样时域图');-0.08-0.06-0.04-0.0200.020.040.060.080.1-8-0.1fs1=150Hz时，fa采样时域图图1150HZ采样频率对信号采样图86420-2-4-6-0.08-0.06-0.04-0.0200.020.040.060.08-0.08-0.06-0.04-0.0200.020.040.060.080.1-8-0.1fs1=150Hz时，fa采样时域图86420-2-4-6-0.08-0.06-0.04-0.0200.020.040.060.080.1fs2=300Hz时，fb采样时域图-8-0.1图2300HZ采样频率对信号采样图、对信号进行快速离散傅里叶变换将两个采样信号进行快速离散傅里叶变换(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);figure(1);plot(w1,xf1),xlabel('fs=150Hz时，fa经fft后频谱图.单位：Hz');

f=40;fs=300;N=300;k=0:N-1;t=-0.1:1/fs:0.1;w2=300*k/Nfb=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');807060504030201050100fs=150Hz时，fa经fft后频谱图.单位：Hz150图3150HZ采样后经FFT后频谱图15010050807060504030201050100fs=150Hz时，fa经fft后频谱图.单位：Hz1501501005050100150200fs=300Hz时，fb经fft后频谱图.单位：Hz250300、信号的重建我们可以通过利用内插法把原信号从采样信号中恢复出来，观察信号在满足怎样的采样条件下能够恢复为原信号，图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/pi)*(ones(length(nTs1),1)*t1-nTs1'*ones(1,length(t1)));figure(1);plot(t1,fa1);axis([-0.10.1-88]);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.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.10.1-88]);xlabel('fs=300Hz时,fb利用内插由样本重建原信号图.’)；grid;-8-0.1-0.08-0.06-0.04-0.0200.02-8-0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1fs=150Hz时，fa利用内插由样本重建原信号图8图6300HZ采样后的信号的重建信号8采样信号的抽取和插值⑴.已知信号X(t)=£9(m+1)cos(2兀(100m+50)t)，以f广20000Hz采样频率采样m=0后，设计抽取因子D和滤波器，分离出如下信号m=1,2,3......50信号组分。参考程序如下：求信号频谱及相应的滤波器。t=0:.00025:1;%Timevectorx=sin(2火pi*30火t)+sin(2火pi*60火t);figure(1)subplot(211)stem(t(1:120),x(1:120));holdony=decimate(x,4);%Viewtheoriginalanddecimatedsignals:stem(x(1:120)),axis([0120-22])%Originalsignaltitle('OriginalSignal')subplot(212)stem(y(1:30))%Decimatedsignaltitle('DecimatedSignal')⑵.已知信号七(t)=^(m+1)coS(10(m+31J),七(t)=尤(m+1)co&(10(h+37》，以m=0m=50fs=10000Hz采样频率采样后，设计插值因子I和滤波器，合成信号气(t)+%(t)。参考程序如下：%File_C3:UpSampl.m%该程序仿真通过零阶保持内插对信号进行上抽样，分析上抽样对信号频谱的影响clearallclcf0=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;fork=1:Nfori=1:IAddsig(I*k-I+i)=sig(k);%对原始数字信号补零内插，抽样时间变为dL『Iendendfigure(1)subplot(2,1,1)t2=[0:N1-1];stem(t，sig);title('Original')xlabel('Sampletime/n')axis([019-1.21.2])subplot(2,1,2)stem(Addsig);axis([080-1.21.2])xlabel('Sampletime/n')title('SignaloutputofZero-orderholdinterpolator')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([01-2060]);ylabel('Amplitude(dB)')xlabel('\omega/\pi')title('AmplitudeSpectrum')subplot(2,1,2)plot(f2,20*log10(abs(AddAm)))axis([01-2060]);ylabel('Amplitude(dB)')xlabel('\omega/\pi')title('AmplitudeSpectrumofZero-orderholdinterpolationsignal')[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([01023-1.21.2])title('outputsignaloutofInterpolator')xlabel('Sampletime/n')subplot(2,1,2)plot(f2,20*log10(abs(UpAm)))axis([01-3050]);ylabel('Amplitude(dB)')xlabel('\omega/\pi')%File_C3:UpSampl.m%该程序仿真通过零阶保持内插对信号进行上抽样，分析上抽样对信号频谱的影响clearallclcf0=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;fork=1:Nfori=1:IAddsig(I*k-I+i)=sig(k);%对原始数字信号补零内插，抽样时间变为dL『Iendendfigure(1)subplot(2,1,1)t2=[0:N1-1];stem(t,sig);title('Original')xlabel('Sampletime/n')axis([019-1.21.2])subplot(2,1,2)stem(Addsig);axis([080-1.21.2])xlabel('Sampletime/n')title('SignaloutputofZero-orderholdinterpolator')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([01-2060]);ylabel('Amplitude(dB)')xlabel('\omega/\pi')title('AmplitudeSpectrum')subplot(2,1,2)plot(f2,20*lo