利用DFT分析信号频谱【信号与系统二】

利用DFT分析信号频谱

• 1. 利用FFT分析信号频谱
• 2. 产生一个淹没在噪声中的信号x(t)，分析信号的频谱
• 3. 利用有限项分析三角波信号
• 4. 利用频域分析方法，重新对第一次实验时的含噪信号进行滤波处理
• 5. 通过频谱分析，识别电话拨号音的电话号码
• 6. 试着将所给的音阶时域波形和频谱图显示出来
• 7. 将女生语音变成男生语音，进行频域分析

1. 利用FFT分析信号频谱

(1) 确定DFT计算的各参数（抽样间隔，截短长度，频谱分辨率等）；
(2) 比较理论值与计算值，分析误差原因，提出改善误差的措施。

f_sam=50;t=1/f_sam;tp=4;n=500;
T=0:t:tp;
w=(-n/2:n/2-1)*(2*pi/n)*f_sam;
x=exp(-2*T);
X=t*fftshift(fft(x,n));
y=1./(1i*w+2);
subplot(2,1,1);
plot(T,x);
xlabel('t');title('时域波形');
subplot(2,1,2);
plot(w,abs(X),w,abs(y),"r-.");
axis([-10,10,0,1]);
title('幅度谱');xlabel('w');
legend("计算值","理论值");

2. 产生一个淹没在噪声中的信号x(t)，分析信号的频谱

f1=50;f2=120;f3=360;
T=0:0.0001:0.1;
x=sin(2*pi*f1*T)+sin(2*pi*f2*T);
noise=0.2*sin(2*pi*f3*T);
y=x+noise;
subplot(2,2,1);
plot(T,x);
subplot(2,2,2);
plot(T,y);

fsam=240;t=1/fsam;tp=0.1;N=length(x);
T=0:t:tp;
w=(-N/2:N/2-1)*(2*pi/N)*fsam;
X=fftshift(fft(y,N));
subplot(2,2,3);
stem(w,abs(X));
axis([-80,80,0,600]);

50Hz和120Hz的正弦成分对应的谱峰位置如红圈所示。根据采样定理对原信号进行采样变成离散信号，再用DFT变换得到检测信号的频谱图。

3. 利用有限项分析三角波信号

t= -2:0.001:2;
n=input('请输入N： ');
cn=zeros(1,length(t));  
for i=1:1:n                        % 偶次谐波为零
    cn=cn+8/(pi*pi)*(-1).^(i-1).*sin((2*i-1)*pi*t)/((2*i-1).^2);
end
plot(t,cn);

4. 利用频域分析方法，重新对第一次实验时的含噪信号进行滤波处理

[x,fs]=audioread('motherLand.wav');
time=(0:length(x)-1)/fs;
subplot(3,2,1);
plot(time,x);
title('原音频')
xlabel('t/s');

fsam=fs;t=1/fsam;tp=10;n=length(x);
T=0:t:tp;
w=(-n/2:n/2-1)*(2*pi/n)*fsam/(2*pi)/1000;
X=fftshift(fft(x));
subplot(3,2,2);
plot(w,abs(X));
title('原音频频谱')
xlabel('f/kHz');
wide=10;
axis([0,wide,0,1800]);

f1=80000;
f2=95000;
noise=(0.15*sin(2*pi*f1*time))'+(0.1*sin(2*pi*f2*time))';
y=x+noise;
audiowrite('before.wav',y,fs);
subplot(3,2,3);
plot(time,y);
title('加入噪声后音频')
xlabel('t/s');

Y=fftshift(fft(y));
subplot(3,2,4);
plot(w,abs(Y));
title('加入噪声后音频频谱')
xlabel('f/kHz');
axis([0,wide,0,1800]);

b=[1.53116389e+03,-1.29990890e-09,7.32176217e+12,-2.03715033e+00,7.71381999e+21];  
a=[1,3.47913978e+04,1.87590501e+09,4.03313474e+13,7.97671668e+17,7.71381999e+21];  
W =linspace(0,14000*2*pi,1000);
H = freqs(b,a,W); 
subplot(3,2,5);
plot(W/(2*pi)/1000,abs(H));
title('滤波器');
xlabel('f/kHz');

sys=tf(b,a);
output=lsim(sys,x,time);
audiowrite('after.wav',output,fs);
OUTPUT=fftshift(fft(output));
subplot(3,2,6);
plot(w,abs(OUTPUT));
title('处理后的音频频谱');
xlabel('f/kHz');
axis([0,wide,0,1800]);

5. 通过频谱分析，识别电话拨号音的电话号码

[x,fs]=audioread('telephone_number.wav');
n=length(x);
time=(0:n-1)/fs;
subplot(6,4,1);
plot(time,x);
title('电话号码音频');xlabel('t/s');
axis([0,6,-0.6,0.6]);
w=(-n/2:n/2-1)*fs/n;
X=fftshift(fft(x));
subplot(6,4,5);
plot(w,abs(X));
axis([0,2500,0,2500]);
title('电话号码频谱');xlabel('f/Hz');

num_ts=0;num_te=1.62;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,2);
plot(num1_t,num1_x);
title('第一个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,6);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第一个数字频谱');xlabel('f/Hz');

num_ts=1.62;num_te=2;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,3);
plot(num1_t,num1_x);
title('第二个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,7);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第二个数字频谱');xlabel('f/Hz');

num_ts=2;num_te=2.6;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,4);
plot(num1_t,num1_x);
title('第三个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,8);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第三个数字频谱');xlabel('f/Hz');

num_ts=2.6;num_te=2.89;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,9);
plot(num1_t,num1_x);
title('第四个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,13);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第四个数字频谱');xlabel('f/Hz');

num_ts=2.89;num_te=3.22;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,10);
plot(num1_t,num1_x);
title('第五个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,14);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第五个数字频谱');xlabel('f/Hz');

num_ts=3.22;num_te=3.73;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,11);
plot(num1_t,num1_x);
title('第六个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,15);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第六个数字频谱');xlabel('f/Hz');

num_ts=3.73;num_te=4.13;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,12);
plot(num1_t,num1_x);
title('第七个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,16);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第七个数字频谱');xlabel('f/Hz');

num_ts=4.13;num_te=4.46;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,17);
plot(num1_t,num1_x);
title('第八个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,21);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第八个数字频谱');xlabel('f/Hz');

num_ts=4.46;num_te=4.86;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,18);
plot(num1_t,num1_x);
title('第九个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,22);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第九个数字频谱');xlabel('f/Hz');

num_ts=4.86;num_te=5.5;
num1_x=x(fs*num_ts+1:fs*num_te);n1=length(num1_x);
num1_t=(0:n1-1)/fs;
subplot(6,4,19);
plot(num1_t,num1_x);
title('第十个数字时域');xlabel('t/s');
w1=(-n1/2:n1/2-1)*fs/n1;
X1=fftshift(fft(num1_x));
subplot(6,4,23);
plot(w1,abs(X1));
axis([0,2500,0,600]);
title('第十个数字频谱');xlabel('f/Hz');

对照标准频率，我们可以通过频谱分析得出电话号码为：5205201314

6. 试着将所给的音阶时域波形和频谱图显示出来

[x,fs]=audioread('yinjie_Adiao.wav');
n=length(x);
time=(0:n-1)/fs;
subplot(2,1,1);
plot(time,x);
title('时域波形');xlabel('t/s');

w=(-n/2:n/2-1)*fs/n;
X=fftshift(fft(x));
subplot(2,1,2);
plot(w,abs(X));
title('频谱图');xlabel('f/Hz');
axis([0,600,0,1200]);

7. 将女生语音变成男生语音，进行频域分析

[x1,fs1]=audioread('为什么要进行频域分析_男声.wav');
x1=x1(:,1);
n1=length(x1);time1=(0:n1-1)/fs1;
[x2,fs2]=audioread('为什么要进行频域分析_女声.wav');
x2=x2(:,1);
n2=length(x2);time2=(0:n2-1)/fs2;
w1=(-n1/2:n1/2-1)*fs1/n1;
X1=fftshift(fft(x1));
w2=(-n2/2:n2/2-1)*fs2/n2;
X2=fftshift(fft(x2));

subplot(2,2,1);
plot(w1,abs(X1));
axis([0,15000,0,800]);
title('男声频谱图');xlabel('f/Hz');
subplot(2,2,2);
plot(w2,abs(X2));
axis([0,15000,0,800]);
title('女声频谱图');xlabel('f/Hz');

fs3=fs2*0.85;
w3=(-n2/2:n2/2-1)*fs3/n2;
subplot(2,2,3);
plot(w3,abs(X2));
axis([0,15000,0,800]);
title('女声变男声频谱图');xlabel('f/Hz');

audiowrite('女声变男声.wav',x1,fs3);