DFT频谱分析作业

DFT频谱分析作业通信工程2013040150作业一设周期序列𝑥̃(𝑛)=cos(0.48𝜋𝑛)+cos(0.52𝜋𝑛)截取N点长得到𝑥(𝑛)=𝑥̃(𝑛)𝑅𝑁(n)（1）N=10，做10点DFT，得到X1(k);（2）N=10，做100点补零DFT,得到X2(k);（3）N=100，做100点DFT,得到X3(k)。要求：针对以上三种情况，分别输出|X1(k)|、|X2(k)|、|X3(k)|的图形，并进行比较、分析和讨论。解：代码如下%双频周期信号的DFT分析n=0:1000;x=cos(0.48*pi*n)+cos(0.52*pi*n);%时域的离散时间信号%-----------N=1010点DFT-------------%N=10;%截取的序列长度N0=10;%DFT点数N0Nxn=x(1:N);%取一个周期的长度Xk=fft(xn,N0);subplot(3,1,1)stem(0:1:N0-1,abs(Xk),'fill','MarkerSize',3.75)%取模画图title('N=1010点DFT')xlabel('k');ylabel('|X1(k)|')grid%----------N=10100点补零DFT----------%N=10;%截取的序列长度N0=100;%DFT点数N0Nxn=x(1:N);%取一个周期的长度Xk=fft(xn,N0);subplot(3,1,2)stem(0:1:N0-1,abs(Xk),'fill','MarkerSize',3.75)%取模画图gridtitle('N=10100点补零DFT')xlabel('k');ylabel('|X2(k)|')%----------N=100100点DFT------------%N=100;%截取的序列长度N0=100;%DFT点数N0Nxn=x(1:N);%取一个周期的长度Xk=fft(xn,N0);subplot(3,1,3)stem(0:1:N0-1,abs(Xk),'fill','MarkerSize',3.75)%取模画图gridtitle('N=100100点DFT')xlabel('k');ylabel('|X3(k)|')运行,得结果如图1所示：图1DFT频谱观察对比：1)10点DFT的|X1(k)|出现明显的栅栏效应，相邻谱线间隔很大，而100点DFT的间隔较小，说明通过增大DFT点数N0即频域抽样点数，可以减小栅栏效应；2)|X2(k)|与|X3(k)|相比出现了明显的频谱泄漏，即截断效应。由于DFT需要对时域序列进行截断，对x(n)乘一个N点长的矩形窗函数RN(n)(主瓣宽度为4π/N)，N越大主瓣越小引起的泄漏就越小。作业二设𝑥𝑎(𝑡)=cos(200𝜋𝑡)+sin(100𝜋𝑡)+cos(50𝜋𝑡)。用FFT分析其频谱结构，选择不同的截取长度Tp，观察截断效应，试用加窗的方法降低旁瓣。解：加矩形窗的FFT代码如下%单频周期信号的DFT分析cos(pi*t)0.5Hzfs=400;%采样频率HZf1=100;f2=50;f3=25;%模拟频率HZn=0:1000;x=cos(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+cos(2*pi*f3*n/fs);%时域的离散时间01234567890510N=1010点DFTk|X1(k)|01020304050607080901000510N=10100点补零DFTk|X2(k)|0102030405060708090100050100N=100100点DFTk|X3(k)|信号%----------Tp=0.04Tp=0.04;N=fs*Tp;%截取的序列长度N0=2048;%DFT点数N0Nxn=x(1:N);%取一个周期的长度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk));subplot(3,1,1)stem(0:1:N0-1,Xk,'fill','MarkerSize',3.75)%取模画图gridtitle('Tp=0.04s');xlabel('k');ylabel('|X(k)|')axis([0220001.2])%----------Tp=4*0.04Tp=4*0.04;N=fs*Tp;%截取的序列长度N0=2048;%DFT点数N0Nxn=x(1:N);%取一个周期的长度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk));subplot(3,1,2)stem(0:1:N0-1,Xk,'fill','MarkerSize',3.75)%取模画图gridtitle('Tp=4*0.04s');xlabel('k');ylabel('|X(k)|')axis([0220001.2])%----------Tp=16*0.04Tp=16*0.04;N=fs*Tp;%截取的序列长度N0=2048;%DFT点数N0Nxn=x(1:N);%取一个周期的长度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk));subplot(3,1,3)stem(0:1:N0-1,Xk,'fill','MarkerSize',3.75)%取模画图gridtitle('Tp=16*0.04s');xlabel('k');ylabel('|X(k)|')axis([0220001.2])运行,得0200400600800100012001400160018002000220000.51Tp=0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=4*0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=16*0.04sk|X(k)|图2DFT频谱图3DFT频谱（模拟频率）由模拟频率图可以看出，在做到Tp=16×0.04S也即N取256时，出现了三个频点，就是原信号的三个频点100π、50π、25π。当Tp=0.04S、4×0.04S也即N取16、64时频率分辨率较低。加Hamming窗进行分析：%单频周期信号的DFT分析cos(pi*t)0.5Hzfs=400;%采样频率HZf1=100;f2=50;f3=25;%模拟频率HZn=0:1000;x=cos(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+cos(2*pi*f3*n/fs);%时域的离散时间信号%----------Tp=0.04Tp=0.04;N=fs*Tp;%截取的序列长度N0=2048;%DFT点数N0Nxn=x(1:N).*hamming(N)';%取一个周期的长度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk));subplot(3,1,1)stem(0:1:N0-1,Xk,'fill','MarkerSize',3.75)%取模画图gridtitle('Tp=0.04s');xlabel('k');ylabel('|X(k)|')axis([0220001.2])0200400600800100012001400160018002000220000.51Tp=0.04s模拟频率/Hz|X(k)|0200400600800100012001400160018002000220000.51Tp=4*0.04s模拟频率/Hz|X(k)|0200400600800100012001400160018002000220000.51Tp=16*0.04s模拟频率/Hz|X(k)|%----------Tp=4*0.04Tp=4*0.04;N=fs*Tp;%截取的序列长度N0=2048;%DFT点数N0Nxn=x(1:N).*hamming(N)';%取一个周期的长度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk));subplot(3,1,2)stem(0:1:N0-1,Xk,'fill','MarkerSize',3.75)%取模画图gridtitle('Tp=4*0.04s');xlabel('k');ylabel('|X(k)|')axis([0220001.2])%----------Tp=16*0.04Tp=16*0.04;N=fs*Tp;%截取的序列长度N0=2048;%DFT点数N0Nxn=x(1:N).*hamming(N)';%取一个周期的长度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk));subplot(3,1,3)stem(0:1:N0-1,Xk,'fill','MarkerSize',3.75)%取模画图gridtitle('Tp=16*0.04s');xlabel('k');ylabel('|X(k)|')axis([0220001.2])0200400600800100012001400160018002000220000.51Tp=0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=4*0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=16*0.04sk|X(k)|图4加Hamming窗的DFT频谱对比可以明显看出，加窗后得到的频谱图上非信号频点处幅度很低。这说明旁瓣确实低了。下图对比更明显：右侧的为加矩形窗，左侧为加海明窗。也可以看出虽然海明窗降低了旁瓣幅度、但是主瓣宽度变宽，主瓣泄漏厉害，本来能够分开的两个频率分量幅值很相近。这显示了加窗的矛盾。作业三读取一个音乐文件并播放，取其中一段数据显示其时域波形，利用DFT进行频谱分析并显示其幅度谱。解：clearall%声音文件的DFT分析%[data,fs,bits]=wavread('warning');4804905005105205305400.780.8Tp=0.04sk|X(k)|4804905005105205305400.981Tp=4*0.04sk|X(k)|4804905005105205305400.981Tp=16*0.04sk|X(k)|4905005105205305405500.840.86Tp=0.04sk|X(k)|5005055105155205255300.980.99Tp=4*0.04sk|X(k)|4905005105205305405500.981Tp=16*0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=4*0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=16*0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=4*0.04sk|X(k)|0200400600800100012001400160018002000220000.51Tp=16*0.04sk|X(k)|[data,fs,bits]=wavread('clock');sound(data,fs)subplot(2,1,1)plot(data);axis([0length(data)-1.21.2])title('声音信号波形图')N=65536;%截取的序列长度N0=65536*2;%DFT点数N0Nxn=data(1:N);%取一个