matlab实现DFT

DFT基于Matlab的实现一、实验目的1．掌握DFT函数的用法。2.利用DFT进行信号检测及谱分析。3．了解信号截取长度对谱分析的影响。二、实验内容1．利用DFT计算信号功率谱。实验程序：t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t)+randn(1,length(t));Y=dft(x,512);P=Y.*conj(Y)/512;f=1000*(0:255)/512;plot(f,P(1:256))2.进行信号检测。分析信号频谱所对应频率轴的数字频率和频率之间的关系。模拟信号)8cos(5)4sin(*2)(tttx，以nt01.010Nn进行取样，求N点DFT的幅值谱。实验程序：subplot(2,2,1)N=45;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title('DFTN=45')subplot(2,2,2)N=50;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title('DFTN=50')subplot(2,2,3)N=55;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title('DFTN=55')subplot(2,2,4)N=60;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title('DFTN=60')3.对2，进一步增加截取长度和DFT点数，如N加大到256，观察信号频谱的变化，分析产生这一变化的原因。在截取长度不变的条件下改变采样频率，观察信号频谱的变化，分析产生这一变化的原因。N加大到256时的程序：N=256;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title('DFTN=256')分析原因：在T=0.01s的情况下，第一个序列的周期是100，第二个序列的周期是50，所以当取样点数小于100时，频率分辨率不够，不能够区分出两个信号。当采样点数足够多（256）时，频率分辨率增加，能够区分出两个频率的信号。将采样间隔变为T=0.1s时，N仍为45的程序：N=45;n=0:N-1;t=0.1*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title('DFTN=45')分析原因：在T=0.1s的情况下，第一个序列的周期是10，第二个序列的周期是5，所以当取样点数为45时，能够区分出两个信号。参数同上，N取64，并在信号中加入噪声w(t)。figure(2)subplot(2,1,1)N=64;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title('DFTN=64')subplot(2,1,2)N=64;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t)+0.8*randn(1,N);y=dft(x,N);plot(q,abs(y));title('DFTN=64(withnoise)')由图可以看出这种噪音不影响信号检测。4.对3，加大噪声到2*randn(1,N)和8*randn(1,N)，画出并比较不同噪声下时域波形和频谱。subplot(2,1,1)N=64;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t)+2*randn(1,N);y=dft(x,N);plot(q,abs(y));title('DFTN=64(withnoise2)')subplot(2,1,2)N=64;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t)+8*randn(1,N);y=dft(x,N);plot(q,abs(y));title('DFTN=64(withnoise8)')subplot(3,2,1)N=64;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t)+0.8*randn(1,N);plot(x);title('噪声为0.8*w的信号')y=dft(x,N);subplot(3,2,2)plot(q,abs(y));title('噪声为0.8*w时的频谱')subplot(3,2,3)N=64;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t)+2*randn(1,N);plot(x);title('噪声为2*w时的信号')y=dft(x,N);subplot(3,2,4)plot(q,abs(y));title('噪声为2*w时的频谱')subplot(3,2,5)N=64;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t)+8*randn(1,N);plot(x);title('噪声为8*w时的信号')y=dft(x,N);subplot(3,2,6)plot(q,abs(y));title('噪声为8*w时的频谱')实验分析：当噪声较小时，不影响信号的检测，但当噪声较大时，就看不出原信号的频率成分了，可以继续加大噪声，可看到其频谱杂乱无章了。5.用一个N点DFT计算两个长度为N的实序列N点离散傅里叶变换，并将结果和直接使用两个N点DFT得到的结果进行比较。x=[123456];y=[654321];[a,b]=dft_2(x,y)a=Columns1through321.0000-3.0000+5.1962i-3.0000+1.7321iColumns4through6-3.0000-3.0000-1.7321i-3.0000-5.1962ib=Columns1through321.00003.0000-5.1962i3.0000-1.7321iColumns4through63.00003.0000+1.7321i3.0000+5.1962i函数文件如下：function[y1,y2]=dft_2(a,b)N=length(a);x=zeros(1,N);x=a+j*b;X=dft(x,N);X0=conj(fliplr(X));X0=[X0(N)X0(1:N-1)];y1=(X+X0)./2;y2=(X-X0)./2./j;直接运行计算：dft(x)ans=Columns1through321.0000-3.0000+5.1962i-3.0000+1.7321iColumns4through6-3.0000-3.0000-1.7321i-3.0000-5.1962idft(y)ans=Columns1through321.00003.0000-5.1962i3.0000-1.7321iColumns4through63.00003.0000+1.7321i3.0000+5.1962i6．比较DFT和DFT的运算时间。(计时函数tic，toc)N分别取256，512，1024，2048，4096，…程序如下：N=256;N=4096;x=randn(1,N);ticy=dft(x,N);tocticz=dft(x);tocN=256：Elapsedtimeis0.172000seconds.Elapsedtimeis0.015000seconds.N=512：Elapsedtimeis0.687000seconds.Elapsedtimeis0.000000seconds.N=1024：Elapsedtimeis3.031000seconds.Elapsedtimeis0.047000seconds.N=2048：Elapsedtimeis13.375000seconds.Elapsedtimeis0.063000seconds.N=4096：Elapsedtimeis59.250000seconds.Elapsedtimeis0.125000seconds.7．对给定语音信号进行谱分析，写出采样频率，画出语音信号的波形及频谱，并分析语音信号的频率分布特点。(1)画时域波形并对整个语音序列做DFT[x,fs]=wavread('C:\ai1.wav');subplot(2,1,1)N=length(x);n=0:N-1;plot(n,x);xlabel('n');ylabel('x');title('时域波形');subplot(2,1,2);N=length(x);n=0:N-1;t=0.01*n;q=n*2*pi/N;y=dft(x,N);plot(q,abs(y));xlabel('n');ylabel('ai1');title('DFT');(2)并分别求出k=300，3500所对应的信号频率（Hz）[x,fs]=wavread('C:\ai1.wav')N=length(x);n=0:N-1;t=n*(1/fs);q=n*2*pi/N;n1=300;q1=n1*2*pi/N;f1=q1*fs/(2*pi);n2=3500;q2=n2*2*pi/N;f2=q2*fs/(2*pi);fs=16000（3）从语音中截取一段语音（256点）做DFT，得出频谱，画出时域波形及频谱。分别求出k=5，60时所对应的信号频率（Hz）程序如下：[x,fs]=wavread('C:\ai1.wav');subplot(2,1,1);N=256;n=0:N-1;x=x(1:256);plot(n,x);xlabel('n');ylabel('x');title('256点时域波形');subplot(2,1,2);N=256;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=x(1:256);y=dft(x,N);plot(q,abs(y));xlabel('n');ylabel('ai1');title('DFT-256');(4)分别求出k=5，60时所对应的信号频率（Hz）的程序。[x,fs]=wavread('C:\ai1.wav');subplot(2,1,1)N=256;n=0:N-1;x=x(1:256);t=0.01*n;q=n*2*pi/N;x=x(1:256);y=dft(x,N);n1=5;q1=n1*2*pi/N;f1=q1*fs/(2*pi);n2=60;q2=n2*2*pi/N;f2=q2*fs/(2*pi);