一、實驗目的

（1）掌握利用DFT近似計算不同類型信號頻譜的原理和方法。

（2）理解誤差產生的原因及減小誤差的方法。

（3）培養學生自主學習能力，以及發現問題、分析問題和解決問題的能力。

二、知識點及背景知識

（1）利用DFT分析連續信號的頻譜，DFT參數 ?

（2） 聲音包括語音、樂音、噪音等。樂音是發音物體有規律地振動而產生的具有固定音高的音，如音樂中的1(Do)、2(Re)、3(Mi)。按照音高順次排列的一串樂音就是音階，如大家熟悉的1(Do )2(Re)3(Mi) 4(Fa)5(So)6(La)7(Si)就是音階。樂音由不同頻率的正弦信號構成，其最簡單的數學模型是cos(2pft)，如C大調音階各樂音對應的頻率如下表：

樂音	1	2	3	4	5	6	7
對應頻率	261.63	293.66	329.63	349.23	392	440	493.88

三、研討內容

1．利用DFT分析x(t)=Acos(2pf1t)+Bsin(2pf2t)的頻譜，其中f1=200Hz，f2=220Hz。分析題目，給出合適的DFT參數，并對實驗結果進行分析，討論窗口的長度和窗口的類型對譜分析有何影響。

?(1)A=B=1; (2)A=1,B=0.1。

• 代碼：

A = 1;B =0.1;f1 =200;f2 = 220;

t = [0:0.001:0.5];

y = A*cos(2*pi*f1*t)+B*sin(2*pi*f2*t);

subplot(2,2,1);plot(t,y);title('原始信號');

axis([0,0.5,-2,2])

y2 = fftshift(fft(y));

fs = linspace(-1000/2,1000/2,length(y));

subplot(2,2,2);plot(fs,abs(y2));title('原始頻譜')

bm = blackman(length(y));

win = y.*bm';

subplot(2,2,3);plot(t,win);axis([0,0.5,-2,2]);

title('加blackman窗后信號')

subplot(2,2,4);

win_fft = fftshift(fft(win));

fs = linspace(-1000/2,1000/2,length(win));

plot(fs,abs(win_fft));

title('加blackman窗后頻譜')

• 結果：

A=B=1

? ? ??A=1,B=0.1

? ?

?

• 分析：

實驗中DFT點數為信號長度，從圖中可以看出，blackman窗的頻譜泄露要比矩形窗(原始帶限信號)的小。

• 代碼：

w0 = 12*pi/64;w1 = 13*pi/64;

k = 0:63;L =64;

xk = cos(w0*k)+1*cos(w1*k);plot(xk);

xk_f = fftshift(fft(xk,L));f1 = (0:L-1)/L;

figure(1);plot(f1,abs(xk_f));title('64點DFT')

k = 0:127;L =128;

xk = cos(w0*k)+1*cos(w1*k);

xk_f = fftshift(fft(xk,L));f1 = (0:L-1)/L;

figure(2);plot(f1,abs(xk_f));title('128點DFT')

k = 0:511;L =512;

xk = cos(w0*k)+1*cos(w1*k);

xk_f = fftshift(fft(xk,L));f1 = (0:L-1)/L;

figure(3);plot(f1,abs(xk_f));title('512點DFT')

• 結果：

?

• 分析：

64點DFT時，兩個譜峰混在一起，無法分辨出來；128點DFT時，兩個譜峰依舊混在一起，無法分辨出來；512點DFT時，兩個譜峰分離開了。因為DFT是對離散信號頻譜DTFT的等間隔抽樣，DFT點數越多，譜線間隔越小，頻譜會顯示更多的細節，也就能夠區分出相鄰的譜峰了。

3．(*)利用DFT分析音階信號yueyin1.wav的頻譜。要求讀取該信號的抽樣頻率，獲得時域抽樣點數N，確定信號的持續時間以及合適的DFT點數，并根據譜分析的結果，判斷是什么調的音階。

C大調對應頻率

樂音	1	2	3	4	5	6	7
對應頻率	261.63	293.66	329.63	349.23	392	440	493.88

• 代碼：

[y,fs] = audioread('yueyin1.wav')%y為時域抽樣點數，fs為抽樣頻率8000Hz

N = length(y);L = N

FFT = fftshift(fft(y,L));Wsam = 2*pi*fs

W = (-Wsam/2+(0:L-1)*Wsam/L)/(2*pi)

plot(W,abs(FFT));axis([0,600,0,1000]);title('yueyin1幅度譜')

• 結果：

• 分析:

從圖中可以看出，yueyin1對應的是C大調。實驗中得出，信號抽樣頻率為8KHz，抽樣點數32000點，從而計算出信號持續時間為32000/8000=4S，這與播放器顯示的時長一致。本實驗中采取的DFT點數為信號時域點數。

（1）利用DFT分析和弦信號hexian1.wav頻譜，確定構成該和弦是哪幾個樂音（即什么頻率分量）

• 代碼：

[y,fs] = audioread('hexian1.wav');

N = length(y);L = N

FFT = fftshift(fft(y,L));Wsam = 2*pi*fs;W = (-Wsam/2+(0:L-1)*Wsam/L)/(2*pi)

figure(1);plot(W,abs(FFT));axis([-500,500,0,1500]);title('hexian1幅度譜')

plot(W,abs(FFT));axis([0,600,0,1000]);title('yueyin1幅度譜')

• 結果：

?

• 分析:

與樂音表對比，該和弦應該是由2、5、6樂音組成。

（2）若樂曲全音符的持續時間為0.2s， 16分音符，從理論上分析利用DFT分析其樂音構成會出現什么問題？設計實驗驗證一下你的判斷，并給出解決問題的方案。

• 代碼：

%0.2s的16分音符

N=100;L=1024;f1=100;f2=200;f3=300;

fs=1000;ws=2*pi*fs;

T=1/fs;t=(0:(N-1))*T;

y=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t);

Y=fftshift(fft(y,L));w=(-ws/2+(0:L-1)*ws/L)/(2*pi);

figure(2);subplot(2,1,1);plot(w,abs(Y));

axis([-350,350,0,80]);title('矩形窗截短后幅度譜');

wh = (hann(N))';y = y.*wh;

Y2 = fftshift(fft(y,L))

subplot(2,1,2);plot(w,abs(Y2));axis([-350,350,0,50]);title('hann窗截短后頻譜')

• 結果：

?

• 分析:

• 代碼：

[y,fs] = audioread('yueyin2.wav')

ws = 2*pi*fs;N = length(y);L = N

w = (-ws/2+(0:L-1)*ws/L)/(2*pi);Y = fftshift(fft(y,L))

plot(w,abs(Y));axis([0,2500,0,1100]);title('yueyin2頻譜')

for i =1:8

???k=y((i-1)*4000+1:i*4000)

???yueyin(k,i,ws)

end

%yueyin.m

function f = yueyin(y,i,ws)

????Y =fftshift(fft(y))

????N = length(y)

????L = N

????w = (-ws/2+(0:L-1)*ws/L)/(2*pi)

????subplot(2,4,i)

????plot(w,abs(Y))

????axis([200,2100,0,1000])

????title(i)

end

• 結果：

???????

?

?

• 分析:

不能直接確定各樂音的頻譜組成，因為頻譜不包含時間信息，不能確定各樂音的頻譜構成。解決方法：將樂音信號按時間分為八小段，再對每一小段進行譜分析，得到相應的頻譜如上圖所示。

從上圖可以得到各個樂音諧波分量，1-7(i)的頻率依次呈遞增趨勢。

總結

以上是生活随笔為你收集整理的近代数字信号处理实验-DFT分析信号的频谱的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。