实验五实验六指导书
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验五 MATLAB 实现DFT
MATLAB 为计算数据的离散快速傅时叶变换,提供了一系列丰富的数学函数,主要有fft 、ifft 、fft2、ifft2和czt 等。当所处理的数据的长度为2的幂次时,采用基-2算法进行计算,计算速度会显著增加。所以,要尽可能使所要处理的数据长度为2幂次或者用添零的方式来添补数据使之成为2的幂次。 1.fft 和ifft 函数 调用格式是: (1)()X fft Y =
如果X 是向量,则采用傅时叶变换来求解X 的离散傅里叶变换;如果X 是矩阵,则计算该矩阵每一列的离散傅里叶变换;如果X 是()D N *维数组,则是对第一个非单元素的维进行离散傅里叶变换。
(2)()N X fft Y ,=
N 是进行离散傅里叶变换的X 的数据长度,可以通过对X 进行补零或截取来实现。 (3)[]()dim ,,X fft Y =或()dim ,,N X fft Y =
在参数dim 指定的维上进行离散傅里叶变换;当X 为矩阵时,dim 用来指定变换的实施方向:dim=1,表明变换按列进行;dim=2,表明变换按行进行。
函数ifft 的参数应用与函数fft 完全相同。 2.fft2和ifft2函数
调用格式是: (1)()X fft Y 2=
如果X 是向量,则此傅里叶变换即变成一维傅里叶变换fft ;如果X 是矩阵,则是计算该矩阵的二维快速傅里叶变换;数据二维傅里叶变换fft 2(X )相当于()()''X fft fft ,即先对X 的列做一维傅里叶变换,然后再对变换结果的行做一维傅里叶变换。
(2)()N M X fft Y ,,2=
通过对X 进行补零或截断,使得X 成为()N M *的矩阵。 函数ifft2的参数应用与函数fft2完全相同
fftn 、ifftn 是对数据进行多维快速傅立变换,其应用与fft2、ifft2类似;在此,不再叙述。 3.czt 函数
调用格式是:
()A W M X czt X ,,,=
式中X 是待变换的时域信号()n x ,其长度设为N ,M 是变换的长度,W 确定变换的步
长,A 确定变换的起点。若M=N ,A=1,则CZT 变成DFT 。
4.fftsfift Y=fftshift(X)
用来重新排列X=fft(x)的输出,当X 为向量时,把X 的左右两半进行交换,从而将零频分量移至频谱的中心;如果X 为二维傅立叶变换的结果,它同时将X 的左右和上下部分进行交换。
5.fftfilt y=fftfilt(b,x)
采用重叠相加法FFT 对信号向量x 快速滤波,得到输出序列向量y ,向量b 为FIR 滤波器的单位脉冲响应,h(n)=b(n+1),n=0,1,…,length(b)-1。
y=fftfilt(b,x,N)
自动选取FFT 长度NF=2^nextpow2(N),输入数据x 分段长度M=NF-length(b)+1,其中nextpow2(N)函数求得一个整数,满足:
2^(nextpow2(N)-1)≤N ≤2^nextpow2(N)
N 缺省时,fftfilt 自动选择合适的FFT 长度NF 和对x 的分段长度M 。
一、利用MATLAB 实现的快速傅里的变换
例1 已知有限长序列()n x 的长度4=N ,且()⎪⎪⎩⎪
⎪⎨
⎧==-===3
3
211
20
1n n n n n x ,用FFT 求()k X ,再用IFFT 求()n x 。
解:利用快速傅里叶变换函数求解的MA TLAB 实现程序清单如下
clear xn=[1,2,-1,3]; X=fft(xn) x=ifft(X)
程序运行结果如下
X =
5.0000 2.0000 + 1.0000i -5.0000 2.0000 - 1.0000i x =
1 2 -1 3
例2 设()n x 是由两个正弦信号及白噪声的叠加,试用FFT 文件对其作频谱分析。 解:程序清单如下
% 产生两个正弦加白噪声;
N=256;
f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs;
x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N); % 应用FFT 求频谱; subplot(2,2,1); plot(x(1:N/4)); title('原始信号'); f=-0.5:1/N:0.5-1/N; X=fft(x); y=ifft(X); subplot(2,2,2);
plot(f,abs(X)); %也可求出将零频分量移至频谱的中心的频谱幅值plot(f,fftshift(abs(X))) title('频域信号'); subplot(2,2,3);
plot(real(x(1:N/4))); %y=x title('时域信号');
运行结果如图1所示,该程序同时完成了傅里叶变换与傅里叶反变换。
图1 例2运行结果
二、 MATLAB 实现序列的移位和周期延拓运算
例3:已知()()n R n x n 88.0=,利用MATLAB 生成并图示(),n x ()m n x -,()()()n R n x N 8,其中N=24,N m <<0,()()8n x 表示()n x 以8为周期的延拓。
图2 序列的移位和周期延拓
解:程序清单如下
N=24;M=8;m=3; %设移位值为3 n=0:N-1;
x1=0.8.^n;x2=[(n>=0)&(n xn=x1.*x2; %产生x(n) xm=xn; nm=n+m %产生x(n-m) xc=xn(mod(n,8)+1); %产生x(n)的周期延拓,求余后加1是因为MA TLAB 向量下标从开始 xcm=xn(mod(n-m,8)+1); %产生x(n)移位后的周期延拓 subplot(2,2,1);stem(n,xn,'.');axis([0,length(n),0,1]);title('x(n)') subplot(2,2,2);stem(nm,xm,'.');axis([0,length(nm),0,1]);title('x(n-m)') subplot(2,2,3);stem(n,xc,'.');axis([0,length(n),0,1]);title('x((n)的周期延拓') subplot(2,2,4);stem(n,xcm,'.');axis([0,length(n),0,1]);title('x(n)的循环移位') 程序运行结果如图2所示。 二、 MATLAB 验证N 点DFT 的物理意义 例4:已知()()n R n x 4=,ω ω ω j j j e e n x FT e X ----==11)]([)(4,绘制相应的幅频和相频曲线, 并计算图示N=8和N=16时的DFT 。 解:程序清单如下 N1=8;N2=16; % 两种FFT 的变换长度 n=0:N1-1;k1=0:N1-1; k2=0:N2-1; w=2*pi*(0:2047)/2048; Xw=(1-exp(-j*4*w))./( 1-exp(-j*w)); %对x(n)的频谱函数采样2048个点可以近似的看作是连续的频谱 xn=[(n>=0)&(n<4)]; %产生x(n) X1k=fft(xn,N1); %计算N1=8点的X 1(k) X2k=fft(xn,N2); %计算N2=16点的X 2(k) subplot(3,2,1);plot(w/pi,abs(Xw));xlabel(‘w/π’) subplot(3,2,2);plot(w/pi,angle(Xw));axis([0,2,-pi,pi]);line([0,2],[0,0]); xlabel(‘w/π’) subplot(3,2,3);stem(k1,abs(X1k),’.’); xlabel(‘k(ω=2πk/N1)’);ylabel(‘|X1(k)|’);hold on plot(N1/2*w/pi,abs(Xw)) %图形上叠加连续频谱的幅度曲线 subplot(3,2,4);stem(k1,angle(X1k)); axis([0,N1,-pi,pi]);line([0,N1],[0,0]);