x1n=exp(j*pi*n/8); %产生x1(n)
X1k=fft(x1n,N); %计算N点DFT[x1(n)] X2k=fft(x1n,N1); %计算N1点DFT[x1(n)] x2n=cos(pi*n/8); %产生x2(n)
X3k=fft(x2n,N); %计算N点DFT[x2(n)] X4k=fft(x2n,N1); %计算N1点DFT[x2(n)]
Subplot(2,2,1);stem(n,abs(X1k),'.');axis([0,20,0,20]);ylabel('|X1(k)|') title('16点的DFT[x1(n)]')
Subplot(2,2,2);stem(n,abs(X3k),'.');axis([0,20,0,20]);ylabel('|X2(k)|') title('16点的DFT[x2(n)]')
Subplot(2,2,3);stem(k,abs(X2k),'.');axis([0,20,0,20]);ylabel('|X1(k)|') title('8点的DFT[x1(n)]')
Subplot(2,2,4);stem(k,abs(X4k),'.');axis([0,20,0,20]);ylabel('|X2(k)|') title('8点的DFT[x2(n)]') 程序运行结果如图5–1所示。
N点离散傅里叶变换的一种物理解释就是X?k?是x?n?以N为周期的周期延拓序列的
~~(k)?X(k)R(k)离散傅里叶级数系数X(k)的主值区间序列，即X。当N?16时， x1(n)N和x2(n)正好分别是ej?n8、cos(?n)的一个周期，所以x1(n)和x2(n)的周期延拓序列就是8j这两个单一频率的正弦序列，其离散傅里叶级数的系数分别如图5–1中的16点的DFT[x1(n)]和16点的DFT[x2(n)]所示。而当N?8时， x1(n)和x2(n)正好分别是e?n8、cos(?n)的8半个周期，所以x1(n)和x2(n)的周期延拓序列就不再是单一频率的正弦序列，而是含有丰富的谐波成分，其离散傅里叶级数的系数与N?16时的差别很大，因此对信号进行谱分析的时候，一定要截取整个周期，否则得到错误的频谱。
16点的DFT[x1(n)]2015201516点的DFT[x2(n)]|X1(k)|1050|X2(k)|051015201050051015208点的DFT[x1(n)]201520158点的DFT[x2(n)]|X1(k)|105005101520|X2(k)|105005101520
图5–1 基本序列的离散傅里叶变换
3．MATLAB验证N点DFT的物理意义
?j4?1?e????(e)?FT[x(n)]??例5–2： 已知x，Xn?Rn4?，绘制相应的幅频和相频曲
1?ejj?线，并计算图示N?8和N?16时的DFT。 解：程序清单如下：
clear all;close all;
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点的X1(k) X2k=fft(xn,N2); %计算N2=16点的X2(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]); xlabel('k(ω=2πk/N1)') ;ylabel('Arg[X1(k)]');hold on
plot(N1/2*w/pi,angle(Xw)) %图形上叠加连续频谱的相位曲线 subplot(3,2,5);stem(k2,abs(X2k),'.');xlabel('k(ω=2πk/N2)');ylabel('|X2(k)|');hold on plot(N2/2*w/pi,abs(Xw))
subplot(3,2,6);stem(k2,angle(X2k),'.');axis([0,N2,-pi,pi]);line([0,N2],[0,0]); xlabel('k(ω=2πk/N2)') ;ylabel('Arg[X2(k)]');hold on plot(N2/2*w/pi,angle(Xw))
程序运行结果如图5–2所示。
42020-200.51w/π1.5200.51w/π1.52420Arg[X1(k)]0246k(ω=2πk/N1)8|X1(k)|20-20246k(ω=2πk/N1)8420Arg[X2(k)]051015k(ω=2πk/N2)2020-20510k(ω=2πk/N2)15|X2(k)|
图5–2 离散傅里叶变换和傅里叶变换的采样关系
序列x?n?的N点DFT的物理意义是对X(ej?)在[0，2π]上进行N点的等间隔采样。图5–2可以直观的看到X?k?与X(e)之间的采样关系。
j?注意：用DFT（或FFT）对模拟信号分析频谱时，最好将X?k?的自变量k换算成对应的模拟频率fk，作为横坐标绘图，便于观察频谱。这样，不管变换区间N取信号周期的几倍，画出的频谱图中有效离散谐波谱线所在的频率值不变。变换公式如下式：
fk?三、实验内容
Fs11k?k?k, k?0,1,2,?,N?1 NNTTp1．设x?n??R6(n)，用MATLAB程序计算N取不同长度(8,32,64)时x?n?的离散傅里叶变换，画图并分析三种情况下异同，并分析原因。
?n?1,0?n?3?2．对以下序列进行谱分析x1(n)??8?n,4?n?7?0,其它n??4?n,0?n?3?x2(n)??n?3,4?n?7
?0,其它n? 选择FFT的变换区间N为8和16 两种情况进行频谱分析。分别画出幅频特性曲线。
并进行对比、分析和讨论。
n，x5(n)?cos(?n/4)?cos(?n/8) ，选择4FFT的变换区间N为8和16 两种情况分别对以上序列进行频谱分析。分别画出其幅频特性曲线。并进行对比、分析和讨论。
3．对以下周期序列进行谱分析x4(n)?cos4．设x(n)?cos(0.48?n)?cos(0.52?n) （1）取0?n?9，求X1?k?；
（2）将（1）中的x?n?补零加长到0?n?99，求X2?k?； （3）增加采样值的个数，取0?n?99，求X3?k?；
（4）要求画出（1）（2）（3）中的时域和频域序列图，分析频域序列之间区别，并分析原因。
?四、实验要求
1．在实验报告中写出完整的自编程序，并给出实验结果。要求自编实验程序不能只采用两
种方法中的一种。
2．实验结果用图像打印显示，结果分析要详细。
3．在绘制图形时要求图形的完整性，每个图形都要有图标标示。 4．在绘制图形时，坐标最好要一致，以方便对比。
五、思考题
1．对于周期序列，如果周期不知道，如何用FFT进行谱分析？ 2．如何选择FFT的变换区间？（包括非周期信号和周期信号）
实验6 序列的离散傅里叶变换（DFT）的应用
一、实验目的
1．加深对序列离散傅里叶变换的理解
2．掌握使用MATLAB实现卷积的两种方法
3．掌握使用MATLAB对连续时间信号进行频谱分析
二．实验方法和示例
1．MATLAB实现快速卷积
例6–1：用快速卷积法计算x(n)?sin(0.4n)R15(n)，h(n)?0.9R20(n)两个序列的卷积。 快速卷积计算原理框图如图6–1所示。
x(n) L点的FFT L点 IDFT h(n) L点的FFT y(n) n
图6–1 快速卷积计算原理框图
解：程序清单如下：
M=15;N=20;nx=1:15;nh=1:20; xn=sin(0.4*nx);hn=0.9.^nh;
L=pow2(nextpow2(M+N-1));%L变为2的整数次幂 Xk=fft(xn,L); Hk= fft(hn,L); Yk=Xk.*Hk;
yn=ifft(Yk,L);ny=1:L
subplot(3,1,1);stem(nx,xn,'.');title('x(n)'); subplot(3,1,2);stem(nh,hn,'.');title('h(n)');
subplot(3,1,3);stem(ny,real(yn),'.');title('y(n)'); 程序运行结果如图6–2所示。
x(n)10-110.5050-505h(n)10150246810y(n)1214161820010203040506070
图6–2 x(n),h(n)及其线性卷积波形
若N和M分别是x(n)和h(n)的长度，FFT的变换长度L必须满足L?N?M?1，输出y(n)才等于x(n)和h(n)的线性卷积。计算两个序列的卷积时，也可以直接调用函数conv来计算，因为MATLAB中的计时比较粗糙，所以只有当N和M较大的时候，才能比较两种方法的执行时间快慢。
2．MATLAB实现用DFT对连续信号作谱分析
例6–2：设xa(t)?cos(200?t)?sin(100?t)?cos(50?t)，用DFT分析xa(t)的频谱结构，选择不同的截取长度TP，观察存在的截断效应。
T?1fs；选取的参数：（1）采样频率fs?400Hz，（2） 采样信号序列x?n??xa?nT?w?n?，
w?n?是窗函数，选取矩形窗为窗函数：函数w?n??RN?n?；（3）对x(n)作2048点DFT作
为xa(t)的近似连续频谱Xa?jf?。其中N为采样点数，N?fs?Tp，Tp为截取时间长度，
取三种长度：0.04s、4?0.04s、8?0.04s。 解：程序清单如下：
fs=400;T=1/fs; %采样频率和采样间隔 Tp=0.04;N=Tp*fs; %采样点数N
N1=[N,4*N,8*N]; %设定三种截取长度 for m=1:3
n=1:N1(m);
xn=cos(200*pi*n*T)+ sin(100*pi*n*T)+ cos(50*pi*n*T); Xk=fft(xn,4096); fk=[0:4095]/4096/T;
subplot(3,1,m);plot(fk,abs(Xk)/max(abs(Xk)));
if m==1 title('时域截取长度为0.04s的信号的幅频特性');end if m==2 title('时域截取长度为4*0.04s的信号的幅频特性');end if m==3 title('时域截取长度为8*0.04s的信号的幅频特性');end end
程序运行结果如图6–3所示。
时域截取长度为0.04s的信号的幅频特性10.5010.5010.50050100150200250300350400时域截取长度为4*0.04s的信号的幅频特性050100150200250300350400时域截取长度为8*0.04s的信号的幅频特性050100150200250300350400