一、实验目的
（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)的频率依次呈递增趋势。
