数字信号处理上机实验dft

本科实验报告实验名称：数字信号处理上机实验作业1：用DFT分析周期序列的频谱任务：设周期序列()cos(0.48)cos(0.52)xnnn截取N点长得到()()()NxnxnRn（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)|的图形，并进行比较、分析和讨论。程序：clearall;n=0:1000;xn=cos(pi*0.48*n)+cos(pi*0.52*n);Xk1=fft(xn(1:10),10);X1=abs(Xk1);subplot(3,1,1);stem(X1,'.');xlabel('k');ylabel('|X1(k)|');title('N=10,10点DFT');Xk2=fft(xn(1:10),100);X2=abs(Xk2);subplot(3,1,2);stem(X2,'.');xlabel('k');ylabel('|X2(k)|');title('N=10,100点补零DFT');Xk3=fft(xn(1:100),100);X3=abs(Xk3);subplot(3,1,3);stem(X3,'.');xlabel('k');ylabel('|X3(k)|');title('N=100,100点DFT');运行结果：1234567891002468k|X1(k)|N=10,10点DFT01020304050607080901000246810k|X2(k)|N=10,100点补零DFT01020304050607080901000102030405060k|X3(k)|N=100,100点DFT分析：从幅度谱中我们可以明显看出，X1(k)的相邻谱线间隔大，栅栏效应明显，频率分辨率低。X2（k）是在采样数据10点后面补零再作DFT，降低了X1(k)中存在的栅栏效应，相邻谱线间隔变小，但并没有提高频率分辨率，没有改变频谱函数的主瓣宽度，只是将X1(k)的已存在的频谱函数进行内插，显示出更多细节。X3(k)与X1(k)相比，信号采样数据增加，幅度谱栅栏效应减小，频率分辨率提高。X2（k）与X3(k)对照，X3(k)相当于时域上乘以一个较长的矩形窗，幅度谱的主瓣与X2(k)相比减小(主瓣宽度为N/4)，频谱泄露明显减轻，频率分辨率提高,但计算复杂度增加。作业2：用DFT对模拟信号进行谱分析任务：设()0.3cos(200)0.5sin(100)cos(50)axtttt用FFT分析其频谱结构，选择不同的截取长度Tp，观察截断效应，试用加窗的方法降低旁瓣。程序：fs=400;n=0:1000;xn=0.3*cos(200*pi*n/fs)+0.5*sin(100*pi*n/fs)+cos(50*pi*n/fs);Vk1=fft(xn(1:16),2048);Vk1=abs(Vk1)./max(abs(Vk1));subplot(3,1,1);plot((0:length(Vk1)-1)*fs/length(Vk1),abs(Vk1));xlabel('HZ');ylabel('|V1(k)|');title('截取时间长度Tp=0.04s，N=16,2048点DFT');Vk2=fft(xn(1:64),2048);Vk2=abs(Vk2)./max(abs(Vk2));subplot(3,1,2);plot((0:length(Vk2)-1)*fs/length(Vk2),abs(Vk2));xlabel('HZ');ylabel('|V2(k)|');title('截取时间长度Tp=4*0.04s，N=64,2048点DFT');Vk3=fft(xn(1:256),2048);Vk3=abs(Vk3)./max(abs(Vk3));subplot(3,1,3);plot((0:length(Vk3)-1)*fs/length(Vk3),abs(Vk3));xlabel('HZ');ylabel('|V3(k)|');title('截取时间长度Tp=16*0.04s，N=256,2048点DFT');运行结果：05010015020025030035040000.20.40.60.81HZ|V1(k)|选取矩形窗函数截取时间长度Tp=0.04s，N=16,2048点DFT05010015020025030035040000.20.40.60.81HZ|V2(k)|截取时间长度Tp=4*0.04s，N=64,2048点DFT05010015020025030035040000.20.40.60.81HZ|V3(k)|截取时间长度Tp=16*0.04s，N=256,2048点DFT程序：fs=400;n=1:1000;xn=0.3*cos(200*pi*n/fs)+0.5*sin(100*pi*n/fs)+cos(50*pi*n/fs);Vk1=fft(xn(1:16).*hamming(16)',2048);subplot(3,1,1);plot((0:length(Vk1)-1)*fs/length(Vk1),abs(Vk1));xlabel('HZ');ylabel('|V1(k)|');title('N=16,2048点DFT');Vk2=fft(xn(1:64).*hamming(64)',2048);subplot(3,1,2);plot((0:length(Vk2)-1)*fs/length(Vk2),abs(Vk2));xlabel('HZ');ylabel('|V2(k)|');title('N=64,2048点DFT');Vk3=fft(xn(1:256).*hamming(256)',2048);subplot(3,1,3);plot((0:length(Vk3)-1)*fs/length(Vk3),abs(Vk3));xlabel('HZ');ylabel('|V3(k)|');title('N=256,2048点DFT');运行结果：05010015020025030035040000.20.40.60.81HZ|V1(k)|选取海明窗函数截取时间长度Tp=0.04s，N=16,2048点DFT05010015020025030035040000.20.40.60.81HZ|V2(k)|截取时间长度Tp=4*0.04s，N=64,2048点DFT05010015020025030035040000.20.40.60.81HZ|V3(k)|截取时间长度Tp=16*0.04s，N=256,2048点DFT分析：从以上两图可以看出，无论是采取矩形窗还是海明窗，随着截取时间长度增大(窗谱压缩变窄)，旁瓣高度降低，能量向主瓣集中，泄露误差减少，同时主瓣变窄。当N增大到256时，能量主要集中在25Hz、50Hz和100Hz，幅度比值为1:0.5:0.3，与题目给出的信号函数相符。300Hz、350Hz和375Hz处为频率延拓过来的，无实际意义。再对比经过矩形窗截取和海明窗截取处理后的图形。程序：fs=400;n=1:1000;xn=0.3*cos(200*pi*n/fs)+0.5*sin(100*pi*n/fs)+cos(50*pi*n/fs);Vk1=fft(xn(1:16).*hamming(16)',2048);Vk1=abs(Vk1)./max(abs(Vk1));subplot(3,1,1);plot((0:length(Vk1)-1)*fs/length(Vk1),abs(Vk1),'y');holdon;Vk1=fft(xn(1:16),2048);Vk1=abs(Vk1)./max(abs(Vk1));subplot(3,1,1);plot((0:length(Vk1)-1)*fs/length(Vk1),abs(Vk1),'g');xlabel('HZ');ylabel('|V1(k)|');title('截取时间长度Tp=0.04s，N=16,2048点DFT');Vk2=fft(xn(1:64).*hamming(64)',2048);Vk2=abs(Vk2)./max(abs(Vk2));subplot(3,1,2);plot((0:length(Vk2)-1)*fs/length(Vk2),abs(Vk2),'y');holdon;Vk2=fft(xn(1:64),2048);Vk2=abs(Vk2)./max(abs(Vk2));subplot(3,1,2);plot((0:length(Vk2)-1)*fs/length(Vk2),abs(Vk2),'g');xlabel('HZ');ylabel('|V2(k)|');title('截取时间长度Tp=4*0.04s，N=64,2048点DFT');Vk3=fft(xn(1:256).*hamming(256)',2048);Vk3=abs(Vk3)./max(abs(Vk3));subplot(3,1,3);plot((0:length(Vk3)-1)*fs/length(Vk3),abs(Vk3),'y');holdon;Vk3=fft(xn(1:256),2048);Vk3=abs(Vk3)./max(abs(Vk3));subplot(3,1,3);plot((0:length(Vk3)-1)*fs/length(Vk3),abs(Vk3),'g');xlabel('HZ');ylabel('|V3(k)|');title('截取时间长度Tp=16*0.04s，N=256,2048点DFT');运行结果：05010015020025030035040000.20.40.60.81HZ|V1(k)|截取时间长度Tp=0.04s，N=16,2048点DFT05010015020025030035040000.20.40.60.81HZ|V2(k)|截取时间长度Tp=4*0.04s，N=64,2048点DFT05010015020025030035040000.20.40.60.81HZ|V3(k)|截取时间长度Tp=16*0.04s，N=256,2048点DFT分析：图中绿线是矩形窗处理后的频谱，黄线是海明窗处理后的频谱。可以看出，用矩形窗截取后的信号主瓣比较集中，但是旁瓣较高，变换中带进了高频干扰和泄漏，频率识别精度较高，幅值识别精度较低。海明窗截取后的信号，其主瓣宽约等于矩形窗的两倍，频率识别精度较低，但旁瓣降低，能量泄露少。可见两种窗函数各有利弊，应用于不同场合。作业3：用DFT分析音乐信号的频谱任务：读取一个音乐文件并播放，取其中一段数据显示其时域波形，利用DFT进行频谱分析并显示其幅度谱。程序：clearall；[data,fs,bits]=wavread('传奇.WAV');sound(data,fs);n=500000;%截取下限N=1000000;%截取上限N0=1000000;%DFT点数subplot(2,1,1);plot(data);axis([nN-22]);title('声音信号波形图');xn=data(n:N);Xk=fft(xn,N0);subplot(2,1,2);plot((0:length(Xk)-1)*fs/length(Xk),abs(Xk));title(['N=',num2str(N-n),'',num2str(N0),'点DFT']);xlabel('hz');ylabel('幅值');grid;运行结果：55.566.577.588.599.510x105-2-1012声音信号波形图00.511.522.533.544.5x10401234x104N=5000001000000点DFThz幅值由于运行的wav文件较大，故只截取了从500