matlab實(shí)現(xiàn)DFT

1、DFT基于Matlab的實(shí)現(xiàn)一、實(shí)驗(yàn)?zāi)康?掌握DFT函數(shù)的用法。2. 利用DFT進(jìn)行信號(hào)檢測(cè)及譜分析。3了解信號(hào)截取長(zhǎng)度對(duì)譜分析的影響。二、實(shí)驗(yàn)內(nèi)容1利用DFT計(jì)算信號(hào)功率譜。實(shí)驗(yàn)程序：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. 進(jìn)行信號(hào)檢測(cè)。分析信號(hào)頻譜所對(duì)應(yīng)頻率軸的數(shù)字頻率和頻率之間的關(guān)系。模擬信號(hào)，以 進(jìn)行取樣，求N點(diǎn)DFT的幅值譜。實(shí)驗(yàn)程序：subplot(2,2,1)

2、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(DFT N=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(DFT N=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=

3、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*sin(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y);title(DFT N=60)3. 對(duì)2，進(jìn)一步增加截取長(zhǎng)度和DFT點(diǎn)數(shù)，如N加大到256，觀察信號(hào)頻譜的變化，分析產(chǎn)生這一變化的原因。在截取長(zhǎng)度不變的條件下改變采樣頻率，觀察信號(hào)頻譜的變化，分析產(chǎn)生這一變化的原因。N加大到256時(shí)的程序：N=256;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*si

4、n(4*pi*t)+5*cos(8*pi*t);y=dft(x,N);plot(q,abs(y);title(DFT N=256)分析原因：在T=0.01s的情況下，第一個(gè)序列的周期是100，第二個(gè)序列的周期是50，所以當(dāng)取樣點(diǎn)數(shù)小于100時(shí)，頻率分辨率不夠，不能夠區(qū)分出兩個(gè)信號(hào)。當(dāng)采樣點(diǎn)數(shù)足夠多（256）時(shí)，頻率分辨率增加，能夠區(qū)分出兩個(gè)頻率的信號(hào)。將采樣間隔變?yōu)門(mén)=0.1s時(shí)，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(DFT N=

5、45)分析原因：在T=0. 1s的情況下，第一個(gè)序列的周期是10，第二個(gè)序列的周期是5，所以當(dāng)取樣點(diǎn)數(shù)為45時(shí)，能夠區(qū)分出兩個(gè)信號(hào)。參數(shù)同上，N取64，并在信號(hào)中加入噪聲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(DFT N=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*ra

6、ndn(1,N); y=dft(x,N);plot(q,abs(y);title(DFT N=64(with noise)由圖可以看出這種噪音不影響信號(hào)檢測(cè)。4. 對(duì)3，加大噪聲到2*randn(1,N)和8*randn(1,N)，畫(huà)出并比較不同噪聲下時(shí)域波形和頻譜。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(DFT N=64(with noise2)subplot(2,1,2)N=64;n=0:

7、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(DFT N=64(with noise8)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的信號(hào))y=dft(x,N);subplot(3,2,2)plot(q,abs(y);title(噪聲為0.8*w時(shí)的頻譜)sub

8、plot(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時(shí)的信號(hào))y=dft(x,N);subplot(3,2,4)plot(q,abs(y);title(噪聲為2*w時(shí)的頻譜)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時(shí)的信號(hào))y=dft(x,N);subp

9、lot(3,2,6)plot(q,abs(y);title(噪聲為8*w時(shí)的頻譜)實(shí)驗(yàn)分析：當(dāng)噪聲較小時(shí)，不影響信號(hào)的檢測(cè)，但當(dāng)噪聲較大時(shí)，就看不出原信號(hào)的頻率成分了，可以繼續(xù)加大噪聲，可看到其頻譜雜亂無(wú)章了。5. 用一個(gè)N點(diǎn)DFT計(jì)算兩個(gè)長(zhǎng)度為N的實(shí)序列N點(diǎn)離散傅里葉變換，并將結(jié)果和直接使用兩個(gè)N點(diǎn)DFT得到的結(jié)果進(jìn)行比較。x=1 2 3 4 5 6; y=6 5 4 3 2 1; a,b=dft_2(x,y) a = Columns 1 through 3 21.0000 -3.0000 + 5.1962i -3.0000 + 1.7321i Columns 4 through 6 -3.

10、0000 -3.0000 - 1.7321i -3.0000 - 5.1962i b = Columns 1 through 3 21.0000 3.0000 - 5.1962i 3.0000 - 1.7321i Columns 4 through 6 3.0000 3.0000 + 1.7321i 3.0000 + 5.1962i 函數(shù)文件如下：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

11、; y2=(X-X0)./2./j; 直接運(yùn)行計(jì)算 ：dft(x) ans = Columns 1 through 3 21.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 3 21.0000 3.0000 - 5.1962i 3.0000 - 1.7321i Columns 4 through 6 3.0000 3.0000 + 1.7321i 3.0000 + 5.19

12、62i 6比較DFT和DFT的運(yùn)算時(shí)間。(計(jì)時(shí)函數(shù) tic， toc)N分別取256，512，1024，2048，4096， 程序如下： N=256;N=4096; x=randn(1,N); tic y=dft(x,N); toc tic z=dft(x); toc N=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 tim

13、e 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對(duì)給定語(yǔ)音信號(hào)進(jìn)行譜分析，寫(xiě)出采樣頻率，畫(huà)出語(yǔ)音信號(hào)的波形及頻譜，并分析語(yǔ)音信號(hào)的頻率分布特點(diǎn)。(1)畫(huà)時(shí)域波形并對(duì)整個(gè)語(yǔ)音序列做DFT x,fs=wavread(C:ai1.wav

14、); subplot(2,1,1) N=length(x); n=0:N-1; plot(n,x); xlabel(n); ylabel(x); title(時(shí)域波形);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所對(duì)應(yīng)的信號(hào)頻率（Hz） x,fs=wavread(C:ai1.wav)N=length(x);n=0:N-1; t=n*(1/fs);q=n*2*pi

15、/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）從語(yǔ)音中截取一段語(yǔ)音（256點(diǎn)）做DFT，得出頻譜，畫(huà)出時(shí)域波形及頻譜。分別求出k=5，60時(shí)所對(duì)應(yīng)的信號(hào)頻率（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點(diǎn)時(shí)域波形); 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);