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

第八章 连续时间傅里叶变换,连续时间傅里叶变换（CTFT）,将连续时间傅里叶级数(CTFS)推广到既能对周期连续时间信号，又能对非周期连续时间信号进行频谱分析。这是一种重要而强有力的方法，因为有很多信号当从时域来看时呈现出很复杂的结构，但从频域来看却很简单。另外，许多LTI系统的特性行为在频域要比在时域容易理解得多。为了更有效地应用频域方法，重要的是要将信号的时域特性是如何与它的频域特性联系起来的建立直观的认识。,频谱计算中的问题,连续离散（抽样，抽样间隔如何选取？） 无穷积分有限长（截断）,8.1 连续时间傅里叶变换的数值近似,傅里叶变换的近似表示,8.2 连续时间信号的采样,8.3 理想抽样信号的傅里叶变换（利用卷积定理）,冲激抽样信号的频谱,说明,抽样定理,8.4 DTFT的引出（利用时移性质）,DTFT:Discrete-time Fourier transform,为研究离散时间系统的频率响应作准备，从抽样信号的傅里叶变换引出：,离散时间信号的傅里叶变换DTFT就是抽样信号的傅立叶变换。,比较,利用快速傅里叶变换计算频谱,fft函数,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 Ts=0.05; t=-5:Ts:5; x=exp(-2*abs(t); subplot(2,1,1); h=plot(t,x); set(h,linewidth,2); xlabel(t /s); ylabel(exp(-2|t|); N=256; w=-pi/Ts+(0:N-1)/N*(2*pi/Ts); X=4./(w.*w+4); subplot(2,1,2); h=plot(w,X); set(h,linewidth,2); xlabel(omega rad/s);ylabel(X(jomega);,抽样间隔如何选取？,(b),% exe4_2_bcde.m clear; T=10;Ts=0.01; t=(-T/2):Ts:(T/2-Ts); N=length(t); x=exp(-2*abs(t);,X=fft(x,N); X=Ts*fftshift(X); w=-pi/Ts+(0:N-1)/N*(2*pi/Ts);,SEMILOGY Semi-log scale plot. SEMILOGY(.) is the same as PLOT(.), except a logarithmic (base 10) scale is used for the Y-axis.,abs_X=4./(4+w.*w); subplot(2,1,1);h=semilogy(w,abs(X); set(h,linewidth,2); xlabel(omega rad/s); ylabel(log_1_0(|X(jomega)|); hold on semilogy(w,abs_X,r:); legend(fft,real); subplot(2,1,2);h=plot(w,unwrap(angle(X); set(h,linewidth,1); xlabel(omega rad/s);ylabel(phi(omega);,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.0. Values outside that range are clipped. Stereo sounds are played, on platforms that support it, when Y is an N-by-2 matrix. SOUND(Y) plays the sound at the default sample rate of 8192 Hz. SOUND(Y,FS,BITS) plays the sound using BITS bits/sample if possible. Most platforms support BITS=8 or 16.,举例,% 8ex4_3_a.m load splat N=8192; y=y(1:N); fs=8192; sound(y,fs); t=(0:N-1)/fs; plot(t,y); xlabel(t/s); ylabel(y);title(waveform);,The egg landed on my cheek with a splat. 鸡蛋啪嚓一声打在我脸上。,(b),% ex4_3_b.m clear; load splat N=length(y); Y=fft(y); Y2=abs(Y); f=(0:N-1)/N*Fs; plot(f,Y2); xlabel(f (Hz);ylabel(|Y|); y1=ifft(Y); figure;plot(abs(y-y1),幅度频谱图,(c),MATLAB实现,% ex4_3_c.m y,fs=wavread(hello.wav); %load splat N=length(y); t=(0:N-1)/fs; subplot(2,1,1); plot(t,y); %sound(y/max(y),fs);pause,Y=fft(y); Y1=conj(Y); y1=ifft(Y1); y1=real(y1); sound(y1,fs); subplot(2,1,2); plot(t,y1);,此图对应 splat,对应 hello.wav,证明,(c1),(c2),实函数,实函数,(d)(e)(f)编程实现,Y2=abs(Y); Y3=exp(j*angle(Y); w=-pi:(2*pi/N):(pi-pi/N); y2=ifft(Y2); y2=real(y2); y3=ifft(Y3); y3=real(y3); sound(y2/max(y2);pause; sound(y3/max(y3);,问题,思考题：如何截断y2信号，即将信号幅度大于1的部分置为1，小于-1的部分置为-1？,思路： 选出大于1的成分； 将大于1的部分置为1。,实现： 选出大于1的成分：position=find(y21); 将大于1的部分置为1：y2(position)=1;,8.6 幅度调制和连续时间傅里叶变换,本地载波,解调,举例,莫尔斯电报编码,(a),% exe4_6_a.m clear; load ctftmod.mat Z=dash dash dot dot; plot(t,Z,r);,(b),freqs(bf,af);,freqs,FREQS 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 B(s) b(1)s + b(2)s + . + b(nb) H(s) = - = - na-1 na-2 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 由欠采样引起的混叠,基本题,MATLAB实现,% exe7_1_a.m T=1/8192; n=0:8191; t=n*T; f0=1000; x=sin(2*pi*f0*t);,(b),取前50个样本：x(1:50),(c),% exe7_1_c.m,fs=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),提示：通过修改exe7_1_c.m中的数据来实现,深入题,% exe7_1_g.m fs=8192; T=1/fs; n=0:fs*10; t=n*T; %f0=3000/2/pi; %bate=2000; f0=100; bate=5000; x=sin(2*pi*f0*t+bate*t.*t/2);,sound(x); specgram(x,8192);,SPECTROGRAM,SPECTROGRAM 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 由样本重建信号,零阶保持 一阶保持 抽样函数,demo,9sam_inversesam_inverse.m,Sa函数作为内插函数（理想化）,Sa函数作为内插函数（理想化）,Sa函数作为内插函数（理想化）,sinc函数内插,sinc函数内插的MATLAB实现,分析：在各抽样值处插入一个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.,举例,clear Ts=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 连续时间傅里叶变换的符号计算,x1=sym(1/2)*exp(-2*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.m clear; 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.m x1=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);,练习1,close all;clear all; syms t a u % exersice 1: Fourier transform of exp(-abs(t) x1=exp(-abs(t); X1=fourier(x1) ezplot(X1,-10,10);axis(-10,10 0 2.1); x11=ifourier(X1,w) figure; x111=simple(x11) ezplot(x111,-10,10);axis(-10,10 0 1.1);,运行结果1,X1 = 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 % exp(-a*t)*heaviside(t) x2=exp(-a*t)*heaviside(t); X2=fourier(x2); a=2; X22=subs(X2) x22=ifourier(X22) figure; subplot(2,1,1);ezplot(abs(X22),-10 10); subplot(2,1,2);ezplot(angle(X22),-10 10);,运行结果2,X22 = 1/(w*i + 2) x22 = exp(-2*x)*heaviside(x) H31 = 2/5 - i/5,练习3,% exersice 3: sinusoidal signals pass the RC low pass filter x3=sin(t)+sin(3*t); H31=subs(X22,w,1) H33=subs(X22,w,3) x33=abs(H31)*sin(t+angle(H31)+abs(H33)*sin(3*t+angle(H33); figure; subplot(2,1,1);ezplot(x3,-10 10); subplot(2,1,2);ezplot(x33,-10 10);,运行结果3,H31 = 2/5 - i/5 H33 = 2/13 - (3*i)/13,练习4,% exersice 4: Ideal low pass filter x4=x3; H4=heaviside(u+2)-heaviside(u-2); H41=subs(H4,u,1); H43=subs(H4,u,3); x44=abs(H41)*sin(t+angle(H41)+abs(H43)*sin(3*t+angle(H43); figure subplot(2,1,1);ezplot(x4,-10 10); subplot(2,1,2);ezplot(x44,-10 10);,运行结果4,练习5,% exersice 5: Hilbert transform x5=sin(t); H5=-i*sign(u); H51=subs(H5,u,1); x55=abs(H51)*sin(t+angle(H51); figure subplot(2,1,1); ezplot(x5,-10 10); subplot(2,1,2); ezplot(x55,-10 10);,运行结果5,函数总结,sound load sinc spline fft/ifft fftshift freqs semilogy conj abs angle,fourier/ifourier dirac heaviside,