离散时间傅里叶变换 DTFT.pdf

数字信号处理第3章DTFT（离散时间傅立叶变换）授课老师：胡双红联系电话：[phone]计算机与通信工程学院第三章 离散时间傅里叶分析?DTFT 离散时间傅里叶变换?DTFT性质?LTI系统的频域表示?采样和模拟信号重建3.1 DTFT 离散时间傅里叶变换?定义定义?计算计算?无限长序列的无限长序列的DFT?有限长序列的有限长序列的DFT?特性：特性：周期性和对称性周期性和对称性?MATLAB实现实现?无限长序列无限长序列?有限长序列有限长序列?验证实序列的周期性验证实序列的周期性3.1 离散时间傅里叶变换（离散时间傅里叶变换（DTFT)?定义：如果定义：如果x(n)绝对可加，则其绝对可加，则其DTFT存在且为存在且为?X(ejw)的的IDTFT为为?其中其中称为数字频率，单位：弧度（称为数字频率，单位：弧度（rad)=nnjjenxnxFeX)()()(=deeXeXFnxnjjj)(21)()(1eg1:求求x(n)=(0.2)nu(n)的离散时间傅里叶变换。的离散时间傅里叶变换。?解：因为序列解：因为序列x(n)绝对可加，所以其傅里叶变换存在且为绝对可加，所以其傅里叶变换存在且为2.02.011)2.0()2.0()()(00=jjjnjnjnnjjeeeeeenxeX65432432234)()(jjjjjjnjjeeeeeeenxeX+=eg2:求有限长序列求有限长序列x(n)=4,3,2,1,2,3,4的离散时间傅里叶变换的离散时间傅里叶变换.?解：重要性质：重要性质：?周期性周期性?对称性对称性:对于对于实值实值序列序列x(n),?或者：或者：)()(2+=jjeXeXsincosjej+=)()(*jjeXeX=奇对称）（偶对称）（奇对称）偶对称）()()(|)(|)(|)(Im)(Im()(Re)(RejjjjjjjjeXeXeXeXeXeXeXeX=MATLAB实现：实现：?对无限长序列求离散时间傅里叶变换：对无限长序列求离散时间傅里叶变换：?先用解析方法得出序列离散时间傅里叶变换的表达式；先用解析方法得出序列离散时间傅里叶变换的表达式；?将数学表达式转化为将数学表达式转化为MATLAB语言；语言；?画出它的幅度和相位（或实部和虚部）画出它的幅度和相位（或实部和虚部）例例3：对例：对例1中的中的X(ejw)在在0，之间的之间的501个等分点上求值，并画出它的幅度、相位、实部和虚部个等分点上求值，并画出它的幅度、相位、实部和虚部?解题步骤：解题步骤：?A.计算计算DTFT表达式表达式?B.将表达式化为将表达式化为MATLAB语言语言X=exp(j*w)./(exp(j*w)-0.2*ones(1,501);?C.画出它的幅度、相位、实部和虚部画出它的幅度、相位、实部和虚部2.0)(=jjjeeeX-MATLAB脚本-w=0:500*pi/500;%0,pi轴等分成轴等分成501个点个点 X=exp(j*w)./(exp(j*w)-0.2*ones(1,501);magX=abs(X);angX=angle(X);%计算计算X的幅度和相位的幅度和相位 realX=real(X);imagX=imag(X);%计算计算X的实部和虚部的实部和虚部 subplot(2,2,1);plot(w/pi,magX);grid xlabel(frequency in pi units);ylabel(Magnitude);title(Manitude Part)subplot(2,2,3);plot(w/pi,angX);grid;xlabel(frequncy in pi units);ylabel(Radians);title(Angle Part)subplot(2,2,2);plot(w/pi,realX);grid xlabel(frequency in pi units);ylabel(Real Part);title(Real Part)subplot(2,2,4);plot(w/pi,imagX);grid xlabel(frequency in pi units);ylabel(Imaginary);title(Imaginary Part)有限长序列的离散时间傅里叶变换有限长序列的离散时间傅里叶变换?直接用直接用MATLAB进行计算。进行计算。?主要步骤：主要步骤：将公式先为转化矩阵乘法，再转化为将公式先为转化矩阵乘法，再转化为MATLAB语言语言设x(n)在n1nn2之间有N个样本，在wk=(/M)k，k=0,1,M点上对X(ejw)求值，则离散时间傅里叶变换x(n)和X(ejwk)被安排为列向量x和X时有X(M+1)?1=W(M+1)?Nx N?1其中我们将k和n)安排成行向量k和n，则有W=exp(-j(/M)kTn)MATLAB中，x(n)和X(ejwk)分别被安排为行向量，则MATLAB:X=x*(exp(-j*pi/M).(n*k)MknxeeXnnnknMjjk,1,0,)()(21)/(=)1(1)1(1)exp()()(+=MNTNTMTknMjxX()MknnneWknMj,1,0,;21/?=板书演示矩阵乘展开形式-MATLAB脚本-n=0:6;x=4,3,2,1,2,3,4;k=0:500;w=(pi/500)*k;X=x*(exp(-j*pi/500).(n*k);magX=abs(X);angX=angle(X);realX=real(X);imagX=imag(X);subplot(2,2,1);plot(k/500,magX);grid xlabel(w);ylabel(幅度);title(Magnitude Part)subplot(2,2,3);plot(k/500,angX/pi);grid xlabel(w);ylabel(相位(rad);title(Angel Part)subplot(2,2,2);plot(k/500,realX);grid xlabel(w);ylabel(实部);title(Real Part)subplot(2,2,4);plot(k/500,imagX);grid xlabel(w);ylabel(虚部);title(Imaginary Part)例4：对例2中的序列x(n)的DTFT，在0，之间的501个等分频率上进行数值求解产生离散时间信号x(n)定义频谱计算的频率间隔DTFT分别计算幅度，相位，实部和虚部 n=0:10;x=(0.9*exp(j*pi/3).n;k=-200:200;w=(pi/100)*k;X=x*(exp(-j*pi/100).(n*k);magX=abs(X);angX=angle(X);subplot(2,1,1);plot(w/pi,magX);grid xlabel(w frequency in units of pi);ylabel(|X|);title(Magnitude Part)subplot(2,1,2);plot(w/pi,angX);grid xlabel(w frequency in units of pi);ylabel(radians/pi);title(Angle Part)例5：设x(n)=0.9exp(j/3)n，0n10,求X(ejw)并研究它的周期性解：因为x(n)是复值序列，故只满足周期性质我们在-2，-2之间的401个频率点上求出并画出它的DTFT以观察它的周期性-MATLAB脚本-产生离散时间信号x(n)定义频谱计算的频率间隔DTFT分别计算幅度和相位例6：设x(n)=(-0.9)n，0n10,求X(ejw)并研究它的周期性解：由DTFT性质得，实序列的DTFT是w的周期函数且共轭对称在-2,2之间的401个频率点上求出并画出它的DTFT以验证周期性和对称性-MATLAB脚本-n=-10:10;x=(-0.9).n;k=-200:200;w=(pi/100)*k;X=x*(exp(-j*pi/100).(n*k);magX=abs(X);angX=angle(X);subplot(2,1,1);plot(w/pi,magX);grid xlabel(frequency in units of pi);ylabel(|X|);title(Magnitude Part)subplot(2,1,2);plot(w/pi,angX/pi);grid xlabel(frequency in units of pi);ylabel(radians/pi);title(Angle Part)产生离散时间信号x(n)定义频谱计算的频率间隔DTFT分别计算幅度和相位3.2 DTFT性质?九条性质九条性质?解析式解析式?数学推导数学推导?MATLAB验证验证?图形法图形法?数值法数值法?1.线性?2.时移?3.频移?4.共轭?5.反转?6.实序列对称性)()()()(2121nxbFnxaFnbxnaxF+=+kjjeeXknxF=)()()()()(00=jnjeXenxF)()(*jeXnxF=)()(jeXnxF=)(Im)()(Re)()()()(jojeoeeXnxFeXnxFnxnxnx=+=则3.2 DTFT性质所有性质板书推导证明过程所有性质板书推导证明过程?7.卷积?8.相乘?9.能量定义能量谱密度)()()()()(*)(212121jjeXeXnxFnxFnxnxF=deXeXnxnxFjj)()(21)()(2121=0222|)(|)(|21|)(|deXdeXnxEjjx2|)(|)(jxeX=用MATLAB验证DTFT的性质（以频移性和时移性为例）例7：用x(n)=cos(n/2)，0n100和y(n)=e jn/4x(n)验证频移性(P46)图形法解：利用MATLAB分别作出x(n)和y(n)的DTFT的幅度与相位-MATLAB脚本-n=0:100;x=cos(pi*n/2);k=-100:100;w=(pi/100)*k;%产生 的频率w X=x*(exp(-j*pi/100).(n*k);%计算x的DTFT y=exp(j*pi*n/4).*x;Y=y*(exp(-j*pi/100).(n*k);%计算y的DTFT%图形验证 subplot(2,2,1);plot(w/pi,abs(X);grid;axis(-1,1,0,60);xlabel(frequency in pi units);ylabel(|X|);title(Magnitude of X);subplot(2,2,2);plot(w/pi,angle(X)/pi);grid;axis(-1,1,-1,1);xlabel(frequency in pi units);ylabel(radiands/pi);title(Angle of X);subplot(2,2,3);plot(w/pi,abs(Y);grid;axis(-1,1,0,60);xlabel(frequency in pi units);ylabel(|Y|);title(Magnitude of Y);subplot(2,2,4);plot(w/pi,angle(Y)/pi);grid;axis(-1,1,-1,1);xlabel(frequency in pi units);ylabel(radiands/pi);title(Angle of Y);由图可见，由图可见，Y(ejw)相对相对X(ejw)确实在幅度和相位上有确实在幅度和相位上有/4的频移的频移验证时移性：数值法例8：设x(n)是在0n10内，0，1间均匀分布的随机序列，并令y(n)=x(n-2)，用MATLAB验证时移性。(P47)-MATLAB脚本-x=rand(1,11);n=0:10;%产生随机序列x(n)k=0:500;w=(pi/500)*k;%定义频率变量w的范围 X=x*(exp(-j*pi/500).(n*k);%计算x(n)的DTFT X(ejw)Y_check=(exp(-j*2).w).*X;%用性质计算x(n-2)的DTFT y,m=sigshift(x,n,2)；%产生x(n)的移位序列y(m)Y=y*(exp(-j*pi/500).(m*k);%用定义求y(m)的DTFT Y(ejw)error=max(abs(Y-Y_check)%求两者之间的差的最大值error=1.7196e-014在MATLAB的有限精度内，这两个数组是一致的3.3 LTI系统的频域表示?频率响应定义频率响应定义?系统对各种输入的响应系统对各种输入的响应?解析表达式解析表达式?仿真验证仿真验证?频率响应函数求法频率响应函数求法?解析法解析法?MATLAB数值求解数值求解?脉冲响应的离散时间傅里叶变换=njjenheH)()(频率响应定义频率响应定义对复指数ejwn的响应h(n)knjkkeAknjjkkkeeHA)(H(ejw)njenx0)(=njjeeHny00)()(=板书推导等式来由板书推导等式来由?对正弦序列的响应?对任意序列的响应H(ejw)+kkkknA)cos(+kjkkjkkkeHneHA)(cos(|)(|H(ejw)(jeX)()()(jjjeXeHeY=解：解：1）利用上面的结论可得：）利用上面的结论可得：2）在稳态下，输入为）在稳态下，输入为x(n)=cos(0.05n)u(n)，w0=0.05,0=0，系统对应频率响应为，系统对应频率响应为结论：该信号在输出端被放大结论：该信号在输出端被放大4.0928倍，移位了倍，移位了3.42个样本个样本()jjeeH=8.011()5377.005.005.00928.48.011jjjeeeH=()()()42.305.0cos0928.45377.005.0cos0928.4=nnyss例例9：LTI系统由下面差分方程描述：系统由下面差分方程描述：y(n)-0.8y(n-1)=x(n)1）求）求H(ejw)2)对输入对输入x(n)=cos(0.05n)u(n)计算并画出稳态响应计算并画出稳态响应yss(n)2）用）用MATLAB画出激励和响应的图形，验证结果画出激励和响应的图形，验证结果-MATLAB脚本脚本-b=1;a=1,-0.8;%输入系数矩阵输入系数矩阵 n=0:100;x=cos(0.05*pi*n);%输入激励序列输入激励序列 y=filter(b,a,x);%计算响应序列计算响应序列 subplot(2,1,1);stem(n,x);grid%画出激励枝干图画出激励枝干图 xlabel(n);ylabel(x(n);title(Input sequence);subplot(2,1,2);stem(n,y);grid%画出响应枝干图画出响应枝干图 xlabel(n);ylabel(y(n);title(Output sequence);?若若LTI系统用差分方程表示时，如何求其频率响应系统用差分方程表示时，如何求其频率响应H(ejw)?当激励当激励x(n)=ejwn时，时，y(n)=H(ejw)ejwn，代入上式并整理得：，代入上式并整理得：=+MmmNllmnxblnyany01)()()(=+=NlljlMmmjmjeaebeH101)(=+MmmnjmNllnjjlnjjebeeHaeeH0)(1)()()(频率响应函数的求法频率响应函数的求法MATLAB实现定义：bm，al，m=0，1，M，l=0,1，2，N和wkH=(b*exp(-j*m*w)./(a*exp(-j*l*w)()2225.05.0121jjjjnjeeeeeH+=例例10：一个三阶低通滤波器由下面差分方程描述：一个三阶低通滤波器由下面差分方程描述y(n)+0.5 y(n-1)+0.25 y(n-2)=x(n)+2 x(n-1)+x(n-2)求出并画出其幅度和相位响应，并验证它是低通滤波器求出并画出其幅度和相位响应，并验证它是低通滤波器解：解：1）数学方法）数学方法由前面分析得：当由前面分析得：当x(n)=ejwn时，时，y(n)=H(ejw)ejwn，代入方程得，代入方程得H(ejw)ejwn+0.5 H(ejw)ejw(n-1)+0.25H(ejw)e jw(n-2)=ejwn+2 e jw(n-1)+e jw(n-2)消去公共因式消去公共因式ejwn并重新整理，得：并重新整理，得：MATLAB解法：-MATLAB脚本-b=1,2,1;a=1,0.5,0.25;%系数向量系数向量 m=0:length(b)-1;l=0:length(a)-1;%时延向量时延向量 K=500;k=0:K;w=pi*k/K;%0，pi以以pi/500为间隔为间隔 num=b*exp(-j*m*w);%计算系统函数分子计算系统函数分子 den=a*exp(-j*l*w);%计算系统函数分母计算系统函数分母 H=num./den;%求出系统函数求出系统函数 magH=abs(H);angH=angle(H);%计算系统函数的模和相角计算系统函数的模和相角 subplot(2,1,1);plot(w/pi,magH);grid%绘制 幅频特性图，加网格 xlabel(w(frequency in pi units);ylabel(|H|);title(MagnitudeResponse)subplot(2,1,2);plot(w/pi,angH/pi);grid xlabel(w(frequency in pi units);ylabel(Phase in pi Radians);title(Phase Response)从幅度响应的图中可见，这确实是个低通滤波器3.4 采样和模拟信号重建采样和模拟信号重建?采样采样?采样定理采样定理?MATLAB仿真仿真?重建重建?理想恢复理想恢复?D/A转换转换?MATLAB仿真仿真?混叠公式：混叠公式：?若有若有x(n)=xa(nTs)则则?即采样信号的频谱等于原信号的频谱在整个频域以（即采样信号的频谱等于原信号的频谱在整个频域以（2/Ts）为周期扩展。）为周期扩展。?如果每个周期互不重叠，则有可能从如果每个周期互不重叠，则有可能从X(ejw)中恢复出原信号的傅里叶变换中恢复出原信号的傅里叶变换Xa(j)，从而恢复出原信号，从而恢复出原信号xa(t)。)2(1)(lTTjXTeXslsasj=采样：此公式通信原理中有介绍，板书推导此公式通信原理中有介绍，板书推导2）带限信号）带限信号?定义：若有一有限频率定义：若有一有限频率0存在，使得存在，使得Xa(j)=0,(|0),该信号是带限的该信号是带限的?若若 0Ts，则，则sssssjTTTTjXTeX2 F0，则信号，则信号xa(t)能从其样本能从其样本x(n)中恢复而不产生混叠。中恢复而不产生混叠。?因此因此 Fs=2F0称为称为 Nyquist 频率。频率。?方法：用一个数组对一个模拟信号进行仿真。方法：用一个数组对一个模拟信号进行仿真。?步骤：步骤：?定义一个足够小的时间间隔定义一个足够小的时间间隔t;则则xG(m)=xa(mt)?用用x G(m)的傅里叶变换来近似代替的傅里叶变换来近似代替x a(t)的傅里叶变换的傅里叶变换?用用MATLAB语言实现。语言实现。Xa=xa*exp(-j*t*w)*Dt=mtmjGmtmjGaemxttemxjX)()()(MATLAB实现例11：设xa(t)=e-1000|t|，求出并画出它的傅里叶变换-MATLAB脚本-%产生模拟信号 Dt=0.00005;t=-0.005:Dt:0.005;xa=exp(-1000*abs(t);%计算连续时间傅里叶变换 Wmax=2*pi*2000;K=500;k=0:K;W=k*Wmax/K;Xa=xa*exp(-j*t*W)*Dt;%在0,4000rad/s上计算Xa(j)Xa=real(Xa);W=-fliplr(W),W(2:501);Xa=fliplr(Xa),Xa(2:501);%在-4000,0rad/s上复制Xa(j)subplot(2,1,1);plot(t*1000,xa);xlabel(t/ms);ylabel(xa(t);title(Analog Signal)subplot(2,1,2);plot(W/(2*pi*1000),Xa*1000);xlabel(f/kHz);ylabel(Xa(jW)*1000);title(Continuous-time Fourier Transform)例例12：为研究在频域数量上的采样效果，对上例中的：为研究在频域数量上的采样效果，对上例中的xa(t)用两种不同的采样频率采样用两种不同的采样频率采样1）在）在Fs=5000样本样本/s对对x a(t)采样得到采样得到x1(n)，求出并画出，求出并画出X1(ejw).2）在）在Fs=1000样本样本/s对对x a(t)采样得到采样得到x2(n)，求出并画出，求出并画出X2(ejw).解：因为解：因为x a(t)的带宽为的带宽为2kHz，Nyquist速率为速率为4000样本样本/s，1）Fs2F0，所以混叠不存在，所以混叠不存在%Discrete-time Signal Ts=0.0002;n=-25:25;x=exp(-1000*abs(n*Ts);%Discrete-time Fourier Tansform K=500;k=0:K;w=pi*k/K;X=x*exp(-j*n*w);X=real(X);w=-fliplr(w),w(2:K+1);X=fliplr(X),X(2:K+1);subplot(2,1,1);plot(t*1000,xa);xlabel(t/ms);ylabel(x1(n)title(Discrete Signal);hold on stem(n*Ts*1000,x);gtext(Ts=0.2msec);hold off subplot(2,1,2);plot(w/pi,X);title(Discrete-time Fourier Transform)xlabel(w(frequency in pi units);ylabel(X1(w);2）Fs=1000%Discrete-time Signal Ts=0.001;n=-5:5;x=exp(-1000*abs(n*Ts);%Discrete-time Fourier Tansform K=500;k=0:K;w=pi*k/K;X=x*exp(-j*n*w);X=real(X);w=-fliplr(w),w(2:K+1);X=fliplr(X),X(2:K+1);subplot(2,1,1);plot(t*1000,xa);xlabel(t/ms);ylabel(x1(n);title(Discrete Signal);hold on stem(n*Ts*1000,x);gtext(Ts=1msec);hold off subplot(2,1,2);plot(w/pi,X);title(Discrete-time Fourier Transform)xlabel(w(frequency in pi units);ylabel(X1(w);?一、理想抽样一、理想抽样?步骤：步骤：?1）将样本转换为一个加权的冲激串；）将样本转换为一个加权的冲激串；?2）将冲激串经由一个带宽限制到）将冲激串经由一个带宽限制到-Fs/2，Fs/2的理想模拟低通滤波器过滤。的理想模拟低通滤波器过滤。?其中其中+=)()1()()0()()1()()(ssnsTtxtxTtxnTtnx冲激串转换理想低通滤波器x(n)xa(n)=nssanTtFcnxtx)(sin)()(xxxcsin)(sin=重建：重建：二、D/A转换?零阶保持（零阶保持（ZOH）内插：）内插：?样本值在样本间隔内一直保持到下一样本被接收为止样本值在样本间隔内一直保持到下一样本被接收为止?用用stairs函数绘图实现函数绘图实现?一阶保持（一阶保持（FOH）内插：）内插：?相邻样本用直线连接相邻样本用直线连接?用用plot函数绘图实现函数绘图实现?三次样条内插：三次样条内插：?相邻样本用三阶多项式描述相邻样本用三阶多项式描述?用用spline函数实现数据求解并输出函数实现数据求解并输出?1）sinc(x)函数：产生函数：产生sin(x)/x函数函数?若在一个很细的栅格若在一个很细的栅格t上内插上内插xa(t),有这就能作为一个矩阵向量乘法运算实现有这就能作为一个矩阵向量乘法运算实现?n=n1:n2;t=t1:t2;Fs=1/Ts;nTs=n*Ts;%Ts为抽样间隔为抽样间隔?xa=x*sinc(Fs*(ones(length(n),1)*t-nTs*ones(1,length(t);?1）利用）利用spline(nTs,x,t)函数实现三次样条内插。函数实现三次样条内插。?2）作图：）作图：?用用stairs(nTs,x)函数画出函数画出ZOH波形波形?用用plot(nTs,x)函数画出函数画出FOH波形波形2121)(sin)()(ttmtnTtmFcnxtmxnnnssa=，MATLAB实现例例13：用：用sinc函数对例函数对例12中的两个样本中的两个样本x1(n)和和x2(n)进行重建，并对结果进行讨论进行重建，并对结果进行讨论解：在解：在-0.005t0.005内以栅格间隔内以栅格间隔0.00005s对对x1(n)和和x2(n)重建重建1）对）对x1(n)进行重建，进行重建，%Discrete-time Signal x1(n)Ts=0.0002;n=-25:25;nTs=n*Ts;x=exp(-1000*abs(nTs);%Analog Signal reconstructionDt=0.00005;t=-0.005:Dt:0.005;Fs=1/Ts;xa=x*sinc(Fs*(ones(length(n),1)*t-nTs*ones(1,length(t);%check error=max(abs(xa-exp(-1000*abs(t)error=0.0363plot(t*1000,xa)；hold on；stem(nTs*1000,x)xlabel(t/ms);ylabel(xa(t);title(the Reconstruction Signal of x1(n)由图可见，重建由图可见，重建成功成功2）对x2(n)进行重建-MATLAB脚本脚本-%Discrete-time Signal x2(n)Ts=0.001;n=-5:5;nTs=n*Ts;x=exp(-1000*abs(nTs);%Analog Signal reconstructionDt=0.00005;t=-0.005:Dt:0.005;Fs=1/Ts;xa=x*sinc(Fs*(ones(length(n),1)*t-nTs*ones(1,length(t);%checkerror=max(abs(xa-exp(-1000*abs(t)error=0.1852plot(t*1000,xa);hold on;stem(nTs*1000,x)xlabel(t/ms);ylabel(xa(t)title(the Reconstruction Signal of x2(n);hold off由图和误差结果可见，由于混叠现象的产生，重建信号由图和误差结果可见，由于混叠现象的产生，重建信号严重失真严重失真例例14：分别利用：分别利用ZOH和和FOH内插画出内插画出x1(n)的重建信号，并对结果进行讨论的重建信号，并对结果进行讨论-MATLAB脚本脚本-%Discrete-time Signal x1(n)Ts=0.0002;n=-25:25;nTs=n*Ts;x=exp(-1000*abs(nTs);%Plotssubplot(2,1,1);stairs(nTs*1000,x);xlabel(t/ms);ylabel(xa(t)title(Reconstructed Signal from x1(n)using ZOH);hold on;stem(nTs*1000,x);hold off;subplot(2,1,2);plot(nTs*1000,x);xlabel(t/ms);ylabel(xa(t)title(Reconstructed Signal from x1(n)using FOH);hold on;stem(nTs*1000,x);hold off;FOH重建效果较好，但在重建效果较好，但在t=0处重现效果不正确处重现效果不正确ZOH重建较为粗糙，必须对模拟信号进一步处理重建较为粗糙，必须对模拟信号进一步处理例例15：用：用spline函数对函数对x1(n)和和x2(n)进行重建，并对结果进行讨论进行重建，并对结果进行讨论解：解：1）对）对x1(n)进行重建进行重建%Discrete-time Signal x1(n)Ts=0.0002;n=-25:25;nTs=n*Ts;x=exp(-1000*abs(nTs);%Analog Signal reconstructionDt=0.00005;t=-0.005:Dt:0.005;xa=spline(nTs,x,t);%checkerror=max(abs(xa-exp(-1000*abs(t)error=0.0317subplot(2,1,1);plot(t*1000,xa);xlabel(t/ms);ylabel(xa(t);title(the Reconstruction Signal of x1(n)hold on;stem(nTs*1000,x);hold off2）对）对x2(n)进行重建进行重建-MATLAB脚本脚本-%Discrete-time Signal x2(n)Ts=0.001;n=-5:5;nTs=n*Ts;x=exp(-1000*abs(nTs);%Analog Signal reconstructionDt=0.00005;t=-0.005:Dt:0.005;xa=spline(nTs,x,t);%checkerror=max(abs(xa-exp(-1000*abs(t)error=0.1679subplot(2,1,2);plot(t*1000,xa);xlabel(t/ms);ylabel(xa(t);title(the Reconstruction Signal of x1(n)hold on;stem(nTs*1000,x);hold off将最大误差和重建图形与将最大误差和重建图形与sinc函数相比较，可见效果较好频域内的混叠现象使得时域波形失真严重，但较函数相比较，可见效果较好频域内的混叠现象使得时域波形失真严重，但较sinc函数的波形好函数的波形好考虑一模拟信号考虑一模拟信号xa(t)=sin(20t)，0t1，在，在Ts=0.01，0.05和和0.1s间隔对它采样得到间隔对它采样得到x(n)1）对每一）对每一Ts画出画出x(n)2）采用）采用sinc内插（用内插（用t=0.001)从样本从样本x(n)重建模拟信号重建模拟信号ya(t)，并从图中求出在，并从图中求出在ya(t)中的频率中的频率3）采用三次样条内插从样本）采用三次样条内插从样本x(n)重建模拟信号重建模拟信号ya(t)，并从你的图中求出在，并从你的图中求出在ya(t)中的频率中的频率解：1）-MATLAB脚本-Dt=0.001;t=0:Dt:1;xa=sin(20*pi*t);%产生模拟信号Ts1=0.01;n1=0:100;x1=sin(20*pi*n1*Ts1);%产生Ts=0.01s的离散时间信号Ts2=0.05;n2=0:20;x2=sin(20*pi*n2*Ts2+pi/2);%Ts=0.05s的离散时间信号Ts3=0.1;n3=0:10;x3=sin(20*pi*n3*Ts3);%产生Ts=0.1s的离散时间信号subplot(4,1,1);plot(t,xa);%画模拟信号subplot(4,1,2);stem(n1*Ts1,x1);%画Ts=0.01s的离散时间信号subplot(4,1,3);stem(n2*Ts2+0.025,x2);%画Ts=0.05s的离散时间信号subplot(4,1,4);stem(n3*Ts3,x3);%画Ts=0.1s的离散时间信号注意临界状态时必须采1才能还原2）-MATLAB脚本-t=0:0.001:1;Fs1=1/Ts1;nTs1=n1*Ts1;ya1=x1*sinc(Fs1*(ones(length(n1),1)*t-nTs1*ones(1,length(t);Fs2=1/Ts2;nTs2=n2*Ts2;ya2=x2*sinc(Fs2*(ones(length(n2),1)*t-nTs2*ones(1,length(t);Fs3=1/Ts3;nTs3=n3*Ts3;ya3=x3*sinc(Fs3*(ones(length(n3),1)*t-nTs3*ones(1,length(t);subplot(3,1,1);plot(t,ya1);subplot(3,1,2);plot(t,ya2);subplot(3,1,3);plot(t,ya3);3）-MATLAB脚本-ya1=spline(nTs1,x1,t);ya2=spline(nTs2,x2,t);ya3=spline(nTs3,x3,t);subplot(3,1,1);plot(t,ya1);subplot(3,1,2);plot(t,ya2);subplot(3,1,3);plot(t,ya3);