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

实验一：低通采样定理和内插与抽取实现一 实验目的1. 连续信号和系统的表示方法，以及坊真方法。2. 用MATLAB实现连续信号采用与重构的方法，3. 采样信号的插值和抽取等重采样实现方法。4. 用时域采样信号重构连续时域信号的原理和方法。5. 用MATLAB绘图函数表示信号的基本方法，实验数据的可视化表示。二 原理1 、时域抽样定理令连续信号xa(t)的傅里叶变换为Xa(j),抽样脉冲序列p(t)傅里叶变换为P(j),抽样后的信号x(t)的傅里叶变换为X(j)若采用均匀抽样,抽样周期Ts,抽样频率为s=2fs,由前面分析可知:抽样的过程可以通过抽样脉冲序列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)进行抽样,可知抽样信号X(t)的频谱X(j)是以s为周期重复。显然,若在抽样的过程中s=2m条件,X(j)才不会产生频谱的混叠,接收端完全可以由x(t)恢复原连续信号xa(t),这就是低通信号抽样定理的核心内容。2、信号的重建从频域看,设信号最高频率不超过折叠频率: Xa(j)=Xa(j) |s/2则理想取样后的频谱就不会产生混叠,故有: 让取样信号x(t)通过一带宽等于折叠频率的理想低通滤波器:H(j)=T |s/2滤波器只允许通过基带频谱,即原信号频谱,故:Y(j)=X(j)H(j)=Xa(j)因此在滤波器的输出得到了恢复的原模拟信号:y(t)=xa(t)从时域上看,上述理想的低通滤波器的脉冲响应为:根据卷积公式可求得理想低通滤波器的输出为:由上式显然可得:则:上式表明只要满足取样频率高于两倍信号最高频率,连续时间函数xa(t)就可用他的取样值xa(nT)来表达而不损失任何信息，这时只要把每一个取样瞬时值与内插函数式相乘求和即可得出xa(t),在每一取样点上,由于只有该取样值所对应的内插函数式不为零,所以各个取样点上的信号值不变。三 内容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); %求得滤波器的最小阶数和截止频率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模拟带通滤波器);hold on;plot(1000 1000,ylim,r);plot(2000 2000,ylim,r);%绘通带边界grid onfigure(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)2.理想矩形滤波器的时域表示clear; clc;fh=100;for I=1:4k=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)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)end3. 连续时间信号的采样和重建、分别用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采样时域图);hold off;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采样时域图); 图1 150HZ采样频率对信号采样图图2 300HZ采样频率对信号采样图、对信号进行快速离散傅里叶变换将两个采样信号进行快速离散傅里叶变换(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 );图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/pi)*(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.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;图5 150HZ采样后的信号的重建信号图6 300HZ采样后的信号的重建信号4.采样信号的抽取和插值t = 0:.00025:1; % Time vectorx = sin(2*pi*30*t) + sin(2*pi*60*t);figure(1)subplot(211)stem(t(1:120),x(1:120);hold ony = 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) % Decimated signaltitle(Decimated Signal)%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:N for i=1:I Addsig(I*k-I+i)=sig(k); %对原始数字信号补零内插，抽样时间变为dLI 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 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)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 -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:N for i=1:I Addsig(I*k-I+i)=sig(k); %对原始数字信号补零内插，抽样时间变为dLI 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(Amp