matlab实现DFT.doc

DFT基于Matlab的实现一、 实验目的1 •掌握DFT函数的用法2. 利用DFT进行信号检测及谱分析3 •了解信号截取长度对谱分析的影响二、 实验内容1 •利用DFT计算信号功率谱实验程序：t=0:0.001:0.6;x=si n( 2*pi*50*t)+si n(2*pi*120*t)+ra ndn (1,le ngth(t)); Y=dft(x,512);P=Y.*conj(Y)/512; f=1000*(0:255)/512;plot(f,P(1:256))2.进行信号检测分析信号频谱所对应频率轴的数字频率和频率之间的关系模拟信号X(t) =2* sin(4p 5cos(8二t),以t =0.01 n 0乞n乞N -1进行取样,求N点DFT勺幅 值谱实验程序：subplot(2,2,1)N=45 ;n=0:N-1;t=0.01* n;q=n *2*pi/N;x=2*si n(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title( 'DFT N=45')subplot(2,2,2)N=50; n=0:N-1;t=0.01* n;q=n*2*pi/N;x=2*si n(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title( 'DFT N=50')subplot(2,2,3)N=55 ;n=0:N-1;t=0.01* n;q=n *2*pi/N;x=2*si n(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title( 'DFT N=55')subplot(2,2,4)N=60; n=0:N-1;t=0.01* n;q=n *2*pi/N;x=2*si n(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title( 'DFT N=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*si n(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title( 'DFT N=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*si n(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title( 'DFT N=45')FF I分析原因：在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*si n(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y));title( 'DFT N=64')subplot(2,1,2)N=64; n=0:N-1;t=0.01* n;q=n *2*pi/N;x=2*si n(4*pi*t)+5*cos(8*pi*t)+0.8*ra ndn (1,N);y=dft(x,N); plot(q,abs(y));title( 'DFT N=64(with noise)' )由图可以看出这种噪音不影响信号检测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*s in (4*pi*t)+5*cos(8*pi*t)+2*ra ndn (1,N);y=dft(x,N);plot(q,abs(y));title( 'DFT N=64(with noise2)' )subplot(2,1,2)N=64; n=0:N-1;t=0.01* n;q=n *2*pi/N;x=2*s in (4*pi*t)+5*cos(8*pi*t)+8*ra ndn (1,N);y=dft(x,N); plot(q,abs(y));title( 'DFT N=64(with noise8)' )FFT hl=M(with noise2)FFT N-64(with noiseS)subplot(3,2,1)N=64; n=O:N-1;t=O.O1* n;q=n *2*pi/N;x=2*s in (4*pi*t)+5*cos(8*pi*t)+0.8*ra ndn (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=O:N-1;t=O.O1* n;q=n *2*pi/N;x=2*s in (4*pi*t)+5*cos(8*pi*t)+2*ra ndn (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*s in (4*pi*t)+5*cos(8*pi*t)+8*ra ndn (1,N); plot(x);title( '噪声为8*w时的信号') y=dft(x,N);subplot(3,2,6)plot(q,abs(y));title( '噪声为 8*w时的频谱')噪声为U旷曲的信号噪声为M频谱□ 2CJ 40 60 BO0 2 4 6 8-10噪声为卉时的信号噪声対于利时的频谱实验分析：当噪声较小时，不影响信号的检测，但当噪声较大时， 就看不出原信号的频率成分了，可以继续加大噪声，可看到其频谱杂乱无章了。

5. 用一个N点DFT计算两个长度为N的实序列N点离散傅里叶变换，并将结果和直 接使用两个N点DFT得到的结果进行比较x=[1 2 3 4 5 6];y=[6 5 4 3 2 1];[a,b]=dft_2(x,y)a =Columns 1 through 321.0000 -3.0000 + 5.1962i -3.0000 + 1.7321iColu mns 4 through 6-3.0000 -3.0000 - 1.7321i -3.0000 - 5.1962ib =Columns 1 through 321.0000 3.0000 - 5.1962i 3.0000 - 1.7321i Columns 4 through 63.0000 3.0000 + 1.7321i 3.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 =Columns 1 through 321.0000 -3.0000 + 5.1962i -3.0000 + 1.7321i Columns 4 through 6-3.0000 -3.0000 - 1.7321i -3.0000 - 5.1962i dft(y)ans =Columns 1 through 321.0000 3.0000 - 5.1962i 3.0000 - 1.7321iColumns 4 through 63.0000 3.0000 + 1.7321i 3.0000 + 5.1962i6. 比较DFT和DFT勺运算时间。

计时函数tic , toe)N分别取 256, 512, 1024, 2048, 4096,…程序如下：N=256;N=4096;x=randn(1,N);ticy=dft(x,N);tocticz=dft(x);tocN=256：Elapsed time is 0.172000 seconds.Elapsed time is 0.015000 seconds.N=512：Elapsed time is 0.687000 seconds.Elapsed time is 0.000000 seconds. N=1024： Elapsed time is 3.031000 seconds.Elapsed time is 0.047000 seconds.N=2048： Elapsed time is 13.375000 seconds.Elapsed time is 0.063000 seconds.N=4096： Elapsed time is 59.250000 seconds.Elapsed time is 0.125000 seconds.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');吋域液形n⑵并分别求出k=300, 3500所对应的信号频率(Hz)[x,fs]=wavread('C:\ai1.wav')N=le ngth(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=n 2*2*pi/N;f2。