MATLAb连续时间傅里叶变换.ppt

尹霄丽Electronic Engineer School of BUPTElectronic Engineer School of BUPT第八章第八章 连续时间傅里叶变换连续时间傅里叶变换连续时间傅里叶变换（连续时间傅里叶变换（CTFT）将连续时间将连续时间将连续时间将连续时间傅里叶级数傅里叶级数傅里叶级数傅里叶级数(CTFS)(CTFS)推广到既能对推广到既能对推广到既能对推广到既能对周期周期周期周期连连连连续时间信号，又能对续时间信号，又能对续时间信号，又能对续时间信号，又能对非周期非周期非周期非周期连续时间信号进行频谱分析。连续时间信号进行频谱分析。连续时间信号进行频谱分析。连续时间信号进行频谱分析。这是一种重要而强有力的方法，因为有很多信号当从时这是一种重要而强有力的方法，因为有很多信号当从时这是一种重要而强有力的方法，因为有很多信号当从时这是一种重要而强有力的方法，因为有很多信号当从时域来看时呈现出很复杂的结构，但域来看时呈现出很复杂的结构，但域来看时呈现出很复杂的结构，但域来看时呈现出很复杂的结构，但从频域来看却很简单从频域来看却很简单从频域来看却很简单从频域来看却很简单。另外，许多另外，许多另外，许多另外，许多LTILTI系统的特性行为在频域要比在时域容易系统的特性行为在频域要比在时域容易系统的特性行为在频域要比在时域容易系统的特性行为在频域要比在时域容易理解得多。为了更有效地应用频域方法，重要的是要将理解得多。为了更有效地应用频域方法，重要的是要将理解得多。为了更有效地应用频域方法，重要的是要将理解得多。为了更有效地应用频域方法，重要的是要将信号的时域特性是如何与它的频域特性信号的时域特性是如何与它的频域特性信号的时域特性是如何与它的频域特性信号的时域特性是如何与它的频域特性联系联系联系联系起来的建立起来的建立起来的建立起来的建立直观的认识直观的认识直观的认识直观的认识。频谱计算中的问题频谱计算中的问题连续连续离散（抽样，抽样间隔如何选取？）离散（抽样，抽样间隔如何选取？）无穷积分无穷积分有限长（截断）有限长（截断）8.1 8.1 连续时间傅里叶变换的数值近似连续时间傅里叶变换的数值近似傅里叶变换的近似表示傅里叶变换的近似表示8.2 8.2 连续时间信号的采样连续时间信号的采样8.3 8.3 理想抽样信号的傅里叶变换（利用卷积定理）理想抽样信号的傅里叶变换（利用卷积定理）冲激抽样信号的频谱冲激抽样信号的频谱说明说明抽样定理抽样定理8.4 DTFT8.4 DTFT的引出（利用时移性质）的引出（利用时移性质）DTFT:Discrete-time Fourier transformDTFT:Discrete-time Fourier transform为为研究离散研究离散研究离散研究离散时间时间系系系系统统的的的的频频率响率响率响率响应应作准作准作准作准备备，从抽，从抽，从抽，从抽样样信信信信号的傅里叶号的傅里叶号的傅里叶号的傅里叶变换变换引出：引出：引出：引出：离散离散离散离散时间时间信号信号信号信号的傅里叶的傅里叶的傅里叶的傅里叶变换变换DTFTDTFT就是抽就是抽就是抽就是抽样样信号的傅立信号的傅立信号的傅立信号的傅立叶叶叶叶变换变换。比较比较利用快速傅里叶变换计算频谱利用快速傅里叶变换计算频谱fftfft函数函数FFT Discrete Fourier transform.FFT(X)is the discrete Fourier transform(DFT)of vector X.For matrices,the FFT operation is applied to each column.For N-D arrays,the FFT operation operates on the first non-singleton dimension.FFT(X,N)is the N-point FFT,padded with zeros if X has less than N points and truncated if it has more.FFT实现的是实现的是DTFT的一个周期的抽样，实际的频谱近的一个周期的抽样，实际的频谱近似为似为fft函数的使用说明函数的使用说明补充说明补充说明例题例题解：解：画图（利用解析式）画图（利用解析式）%ss8_2.m and double_side_exp_spectrum.m%ss8_2.m and double_side_exp_spectrum.mTs=Ts=0.050.05;t=-5:Ts:5;t=-5:Ts:5;x=x=exp(-2*abs(t)exp(-2*abs(t);subplot(2,1,1);subplot(2,1,1);h=plot(t,x);set(h,linewidth,2);h=plot(t,x);set(h,linewidth,2);xlabel(t/s);ylabel(exp(-2|t|);xlabel(t/s);ylabel(exp(-2|t|);N=256;N=256;w=-pi/Ts+(0:N-1)/N*(2*pi/Ts)w=-pi/Ts+(0:N-1)/N*(2*pi/Ts);X=X=4./(w.*w+4)4./(w.*w+4);subplot(2,1,2);h=plot(w,X);subplot(2,1,2);h=plot(w,X);set(h,linewidth,2);set(h,linewidth,2);xlabel(xlabel(omegaomega rad/s);ylabel(X(jomega);rad/s);ylabel(X(jomega);抽样间隔如何选取？抽样间隔如何选取？(b)%exe4_2_bcde.m%exe4_2_bcde.mclear;clear;T=10;Ts=0.01;T=10;Ts=0.01;t=(-T/2):Tst=(-T/2):Ts:(T/2-Ts):(T/2-Ts);N N=length(t);=length(t);x=exp(-2*abs(t);x=exp(-2*abs(t);X=X=fftfft(x,N);(x,N);X=Ts*X=Ts*fftshiftfftshift(X);(X);w=-pi/Ts+(0:N-1)/N*(2*pi/Ts);w=-pi/Ts+(0:N-1)/N*(2*pi/Ts);SEMILOGY Semi-log scale plot.SEMILOGY Semi-log scale plot.SEMILOGY(.)is the same as PLOT(.),except a SEMILOGY(.)is the same as PLOT(.),except a logarithmic(base 10)scale is used for the Y-axis.logarithmic(base 10)scale is used for the Y-axis.abs_X=4./(4+w.*w)abs_X=4./(4+w.*w);subplot(2,1,1);h=subplot(2,1,1);h=semilogysemilogy(w,abs(X);(w,abs(X);set(h,linewidth,2);set(h,linewidth,2);xlabelxlabel(omega(omega rad/s);rad/s);ylabel(ylabel(log_1_0(|X(jomega)|)log_1_0(|X(jomega)|););hold onhold onsemilogy(w,abs_X,r:);semilogy(w,abs_X,r:);legend(fft,real);legend(fft,real);subplot(2,1,2);h=plot(w,subplot(2,1,2);h=plot(w,unwrapunwrap(angleangle(X);(X);set(h,linewidth,1);set(h,linewidth,1);xlabel(omega rad/s);ylabel(phi(omega);xlabel(omega rad/s);ylabel(phi(omega);8.5 8.5 连续时间傅里叶变换性质连续时间傅里叶变换性质目的：直观、深刻地理解傅里叶变换的性质；目的：直观、深刻地理解傅里叶变换的性质；主要内容：主要内容：奇偶虚实性；信号的幅度谱与相位谱奇偶虚实性；信号的幅度谱与相位谱尺度变换特性尺度变换特性频移性质和调制定理；频移性质和调制定理；抽样信号的重建抽样信号的重建方法方法sound函数函数SOUND Play vector as sound.SOUND(Y,FS)sends the signal in vector Y(with sample frequency FS)out to the speaker on platforms that support sound.Values in Y are assumed to be in the range-1.0=y 1);position=find(y21);position=find(y21);position=find(y21);2.2.将大于将大于将大于将大于1 1 1 1的部分置为的部分置为的部分置为的部分置为1 1 1 1：y2(position)=1;y2(position)=1;y2(position)=1;y2(position)=1;8.6 8.6 幅度调制和连续时间傅里叶变换幅度调制和连续时间傅里叶变换本地本地本地本地载载波波波波解调解调举例举例莫尔斯电报编码莫尔斯电报编码A .-A .-H .H .O -O -V -V -B -B -I .I .P .-.P .-.W .-W .-C -.-.C -.-.J .-J .-Q -.-Q -.-X -.-X -.-D -.D -.K -.-K -.-R .-.R .-.Y-.-Y-.-E .E .L .-.L .-.S S Z -.Z -.F .-.F .-.M -M -T -T -G -.G -.N -.N -.U .-U .-(a)%exe4_6_a.m%exe4_6_a.mclear;clear;load ctftmod.matload ctftmod.matZ=Z=dash dash dot dotdash dash dot dot;plot(t,Z,r);plot(t,Z,r);(b)freqs(bf,af);freqs(bf,af);freqsFREQS Laplace-transform(s-domain)frequency response.H=FREQS(B,A,W)returns the complex frequency response vector H of the filter B/A:given the numerator and denominator coefficients in vectors B and A.The frequency response is evaluated at the points specified in vector W(in rad/s).The magnitude and phase can be graphed by calling FREQS(B,A,W)with no output arguments.传输函数传输函数 nb-1 nb-2 nb-1 nb-2 B(s)b(1)s +b(2)s +.+b(nb)B(s)b(1)s +b(2)s +.+b(nb)H(s)=-=-H(s)=-=-na-1 na-2 na-1 na-2 A(s)a(1)s +a(2)s +.+a(na)A(s)a(1)s +a(2)s +.+a(na)B,A 矩阵的写法矩阵的写法例题例题运行结果运行结果其他用法其他用法H,W=FREQS(B,A)automatically picks a set of 200 frequencies W on which the frequency response is computed.FREQS(B,A,N)picks N frequencies.See also logspace,polyval,invfreqs,and freqz（离散系统）（离散系统）.(c)分析分析(d)(e)相干接收相干接收需要使用本地载波（接收端）需要使用本地载波（接收端）同步解调：本地载波与发送端载波同频同相同步解调：本地载波与发送端载波同频同相正交调制技术简介正交调制技术简介第一种情况：本地载波与调制载波同频同相第一种情况：本地载波与调制载波同频同相高频信号高频信号高频信号高频信号恢复出的原始恢复出的原始信号信号第二种情况：本地载波与调制载波同频不同第二种情况：本地载波与调制载波同频不同相相 只有高频信号，只有高频信号，经过低通滤波器经过低通滤波器后被滤除？后被滤除？第三种情况：本地载波与调制载波不同频第三种情况：本地载波与调制载波不同频差拍信号差拍信号高频信号高频信号y=x.*cos(2*pi*f1*t);D -.y=x.*sin(2*pi*f1*t);P .-.y=x.*cos(2*pi*f2*t);y=x.*sin(2*pi*f2*t);S .总结总结本地本地本地本地载载波波波波8.7 8.7 由欠采样引起的混叠由欠采样引起的混叠基本题基本题MATLABMATLAB实现实现%exe7_1_a.mT=1/8192;n=0:8191;t=n*T;f0=1000;x=sin(2*pi*f0*t);(b)(b)取前取前5050个样本：个样本：x(1:50)x(1:50)(c)(c)%exe7_1_c.m%exe7_1_c.mfs=8192;T=1/fs;f0=800;W=2*pi*f0*T;n=0:fs;x=sin(W*n);sound(x,fs);X=fft(x,56);stem(abs(X);(d)(d)提示：提示：通过修改通过修改exe7_1_c.mexe7_1_c.m中的数据来实现中的数据来实现深入题深入题%exe7_1_g.m%exe7_1_g.mfs=8192;fs=8192;T=1/fs;T=1/fs;n=0:fs*10;n=0:fs*10;t=n*T;t=n*T;%f0=3000/2/pi;%f0=3000/2/pi;%bate=2000;%bate=2000;f0=100;f0=100;bate=5000;bate=5000;x=x=sin(2*pi*f0*t+bate*t.*t/2)sin(2*pi*f0*t+bate*t.*t/2);sound(x);sound(x);specgram(x,8192);specgram(x,8192);SPECTROGRAMSPECTROGRAM Spectrogram using a Short-Time Fourier Transform(STFT,短时傅里叶变换短时傅里叶变换).S=SPECTROGRAM(X)returns the spectrogram of the signal specified by vector X in the matrix S.By default,X is divided into eight segments with 50%overlap,each segment is windowed with a Hamming window.The number of frequency points used to calculate the discrete Fourier transforms is equal to the maximum of 256 or the next power of two greater than the length of each segment of X.8.7 8.7 由样本重建信号由样本重建信号零零零零阶阶保持保持保持保持一一一一阶阶保持保持保持保持抽抽抽抽样样函数函数函数函数demo9sam_inversesam_inverse.m9sam_inversesam_inverse.mSa函数作为内插函数（理想化）函数作为内插函数（理想化）SaSa函数作为内插函数（理想化）函数作为内插函数（理想化）SaSa函数作为内插函数（理想化）函数作为内插函数（理想化）sincsinc函数内插函数内插sincsinc函数内插的函数内插的MATLABMATLAB实现实现分析：在各抽样值处插入一个分析：在各抽样值处插入一个sinc函数，大小与抽样值成正函数，大小与抽样值成正比，定义域为全时域（或给定定义域）。比，定义域为全时域（或给定定义域）。时间矩阵：时间矩阵：tt=ones(length(n),1)*t-Ts*n*ones(1,length(t)内插函数矩阵：内插函数矩阵：sinc(fs*tt)函数内插：函数内插：x*sinc(tt)%x为样值函数为样值函数内插函数矩阵内插函数矩阵spline:三次样条内插函数三次样条内插函数SPLINE Cubic spline data interpolation.PP=SPLINE(X,Y)provides the piecewise polynomial form of the cubic spline interpolant to the data values Y at the data sites X,for use with the evaluator PPVAL and the spline utility UNMKPP.X must be a vector.If Y is a vector,then Y(j)is taken as the value to be matched at X(j),hence Y must be of the same length as X -see below for an exception to this.If Y is a matrix or ND array,then Y(:,.,:,j)is taken as the value to be matched at X(j),hence the last dimension of Y must equal length(X)-see below for an exception to this.YY=SPLINE(X,Y,XX)is the same as YY=PVAL(SPLINE(X,Y),XX),thus providing,in YY,the values of the interpolant at XX.For information regarding the size of YY see PPVAL.举例举例clearTs=1;Fs=1/Ts;n=0:10;x=sin(n);t=0:.25:10;x_spline=spline(n,x,t);plot(t,x_spline,b);nTs=0:10;tt=ones(length(n),1)*t-nTs*ones(1,length(t);x_sinc=x*sinc(Fs*tt);hold on;plot(t,x_sinc,r);legend(spline,sinc);hold on；stem(n,x,m);hh=findobj(0,type,line);set(hh,linewidth,2);结果图形结果图形8.8 8.8 连续时间傅里叶变换的符号计算连续时间傅里叶变换的符号计算x1=sym(1/2)*exp(-2*t)*x1=sym(1/2)*exp(-2*t)*heavisideheaviside(t);(t);x2=sym(exp(-4*t)*heaviside(t);x2=sym(exp(-4*t)*heaviside(t);heaviside：单位阶跃函数：单位阶跃函数help heaviside HEAVISIDE Unit Step function f=Heaviside(t)returns a vector f the same size as the input vector,where each element of f is 1 if the corresponding element of t is greater than zero.举例：举例：syms t;y=cos(t)*(heaviside(t+0.5*pi)-heaviside(t-0.5*pi);ezplot(y);解：解：解：解：主要代码主要代码%exe4_7_a.mclear;x1=sym(1/2)*exp(-2*t)*heaviside(t);x2=sym(exp(-4*t)*heaviside(t);subplot(2,1,1);ezplot(x1,0,2);legend(x1);axis(0 2 0 1);subplot(2,1,2);ezplot(x2,0,2);legend(x2);axis(0 2 0 1);fourier函数函数FOURIER Fourier integral transform.F=FOURIER(f)is the Fourier transform of the sym scalar f with default independent variable x.F(w)=int(f(x)*exp(-i*w*x),x,-inf,inf)See also sym/ifourier,sym/laplace,sym/ztrans.主要代码主要代码%exe4_7_a.mx1=sym(1/2)*exp(-2*t)*heaviside(t);x2=sym(exp(-4*t)*heaviside(t);X1=fourier(x1);X2=fourier(x2);subplot(2,1,1);ezplot(abs(X1),-20,20);legend(|X1|)axis(-20 20 0 0.3);subplot(2,1,2);ezplot(abs(X2),-20,20);legend(|X2|)axis(-20 20 0 0.3);练习练习1close close allall;clear;clear allall;symssyms t t a a u u%exersice 1:Fourier transform of exp(-abs(t)%exersice 1:Fourier transform of exp(-abs(t)x1=exp(-abs(t);x1=exp(-abs(t);X1=X1=fourierfourier(x1)(x1)ezplotezplot(X1,-10,10);axis(-10,10 0 2.1);(X1,-10,10);axis(-10,10 0 2.1);x11=x11=ifourierifourier(X1,(X1,ww)figure;figure;x111=x111=simplesimple(x11)(x11)ezplot(x111,-10,10);axis(-10,10 0 1.1);ezplot(x111,-10,10);axis(-10,10 0 1.1);运行结果运行结果1X1=2/(w2+1)x11=(2*pi*exp(-w)*heaviside(w)+2*pi*heaviside(-w)*exp(w)/(2*pi)x111=exp(-w)*heaviside(w)+heaviside(-w)*exp(w)练习练习2%exersice 2:Fourier transform of%exersice 2:Fourier transform of%exp(-a*t)*heaviside(t)%exp(-a*t)*heaviside(t)x2=exp(-a*t)*x2=exp(-a*t)*heavisideheaviside(t);(t);X2=X2=fourierfourier(x2);(x2);a=2;a=2;X22=X22=subssubs(X2)(X2)x22=x22=ifourierifourier(X22)(X22)figure;figure;subplot(2,1,1);ezplot(abs(X22),-10 10);subplot(2,1,1);ezplot(abs(X22),-10 10);subplot(2,1,2);ezplot(angle(X22),-10 10);subplot(2,1,2);ezplot(angle(X22),-10 10);运行结果运行结果2X22=X22=1/(w*i+2)1/(w*i+2)x22=x22=exp(-2*x)*heaviside(x)exp(-2*x)*heaviside(x)H31=H31=2/5-i/5 2/5-i/5练习练习3%exersice 3:sinusoidal signals pass the RC low%exersice 3:sinusoidal signals pass the RC low pass filterpass filterx3=sin(t)+sin(3*t);x3=sin(t)+sin(3*t);H31=subs(X22,H31=subs(X22,ww,1),1)H33=subs(X22,H33=subs(X22,ww,3),3)x33=abs(H31)*sin(t+angle(H31)+abs(H33)*sin(3*tx33=abs(H31)*sin(t+angle(H31)+abs(H33)*sin(3*t+angle(H33);+angle(H33);figure;figure;subplot(2,1,1);ezplot(x3,-10 10);subplot(2,1,1);ezplot(x3,-10 10);subplot(2,1,2);ezplot(x33,-10 10);subplot(2,1,2);ezplot(x33,-10 10);运行结果运行结果3H31=H31=2/5-i/52/5-i/5 H33=H33=2/13-2/13-(3*i)/13(3*i)/13练习练习4%exersice 4:Ideal low pass filter%exersice 4:Ideal low pass filterx4=x3;x4=x3;H4=heaviside(u+2)-heaviside(u-2);H4=heaviside(u+2)-heaviside(u-2);H41=subs(H4,H41=subs(H4,uu,1);,1);H43=subs(H4,H43=subs(H4,uu,3);,3);x44=abs(H41)*sin(t+angle(H41)+abs(H43)*sin(3*tx44=abs(H41)*sin(t+angle(H41)+abs(H43)*sin(3*t+angle(H43);+angle(H43);figurefiguresubplot(2,1,1);ezplot(x4,-10 10);subplot(2,1,1);ezplot(x4,-10 10);subplot(2,1,2);ezplot(x44,-10 10);subplot(2,1,2);ezplot(x44,-10 10);运行结果运行结果4练习练习5%exersice 5:Hilbert transform%exersice 5:Hilbert transformx5=sin(t);x5=sin(t);H5=-i*sign(u);H5=-i*sign(u);H51=subs(H5,H51=subs(H5,uu,1);,1);x55=abs(H51)*sin(t+angle(H51);x55=abs(H51)*sin(t+angle(H51);figurefiguresubplot(2,1,1);subplot(2,1,1);ezplot(x5,-10 10);ezplot(x5,-10 10);subplot(2,1,2);subplot(2,1,2);ezplot(x55,-10 10);ezplot(x55,-10 10);运行结果运行结果5函数总结函数总结soundloadsincsplinefft/ifftfftshiftfreqssemilogyconjabsanglefourier/ifourierfourier/ifourierdiracdiracheavisideheaviside