应用MatLab对语音信号进行频谱分析及滤波

探究性学习之二
--应用MatLab
对语音信号进行频谱分析及滤波
一、课题目的
1．通过课下的自我的学习，加深对书本理论知识的理解，提升自身的实际应用能力；
2．巩固所学的数字信号处理理论知识，使自身对信号的采集、处理、传输、显示和存储等有一个系统的掌握和理解；
3．培养自我学习的能力和对相关课程的兴趣；
二、课题内容
录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。

三、具体实现
1．语音信号的采集
利用Windows下的录音机，录制一段自己的话音，时间在1 s内。

然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。

通过wavread函数的使用，我们很快理解了采样频率、采样位数等概念。

2．语音信号的频谱分析
要求学生首先画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性，从而加深学生对频谱特性的理解。

程序如下：
clc;clear;close all;
fs=22050; %语音信号采样频率为22050
x1=wavread('wo.wav'); %读取语音信号的数据，赋给变量x1
sound(x1,22050); %播放语音信号
y1=fft(x1,1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;figure(1);
plot(x1) %做原始语音信号的时域图形
title('原始语音信号');xlabel('time n');ylabel('fuzhi n');
figure(2);freqz(x1) %绘制原始语音信号的频率响应图
title('频率响应图')
figure(3);subplot(2,1,1);
plot(abs(y1(1:512))) %做原始语音信号的FFT频谱图？？
title('原始语音信号FFT频谱');
subplot(2,1,2);plot(f,abs(y1(1:512)));％同一个图？？？
title('原始语音信号频谱')
xlabel('Hz');ylabel('fuzhi'); 波形如下：
-1-0.8-0.6-0.4-0.200.2
0.40.60.81原始语音信号
time n
f u z h i n
-10
-505
4
Normalized Frequency (⨯π rad/sample)
P h a s e (d e g r e e s
)
-20
0204060
Normalized Frequency (⨯π rad/sample)
M a g n i t u d e (d B )
频率响应图
010203040原始语音信号FFT 频谱
0102030
40原始语音信号频谱
Hz
f u z h i
3． 滤波电路设计
给出低通滤波器性能指标fb ＝1000 Hz ，fc ＝1200 Hz ，Rs ＝15 dB ，Rp ＝1 dB 。

并使语音信号经 过该滤波器。

程序如下：
clc;clear;close all;
fs=22050;x1=wavread('wo.wav'); t=0:1/22050:(size(x1)-1)/22050;
Au=0.03;d=[Au*cos(2*pi*5000*t)]';x2=x1+d; wp=0.25*pi;ws=0.3*pi;Rp=1;Rs=15; Fs=22050;Ts=1/Fs;
wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标 ws1=2/Ts*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s'); %选择滤波器的最小阶数 [Z,P,K]=buttap(N); %创建butterworth 模拟滤波器 [Bap,Aap]=zp2tf(Z,P,K);[b,a]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 [H,W]=freqz(bz,az); %绘制频率响应曲线 figure(1);plot(W*Fs/(2*pi),abs(H));grid
xlabel('频率／Hz');ylabel('频率响应幅度');title('Butterworth') f1=filter(bz,az,x2);
figure(2);subplot(2,1,1);plot(t,x2); %画出滤波前的时域图 title('滤波前的时域波形');
subplot(2,1,2);plot(t,f1); %画出滤波后的时域图 title('滤波后的时域波形'); sound(f1,22050); %播放滤波后的信号 F0=fft(f1,1024);f=fs*(0:511)/1024; figure(3) y2=fft(x2,1024);
subplot(2,1,1);plot(f,abs(y2(1:512))); %画出滤波前的频谱图 title('滤波前的频谱');xlabel('Hz');ylabel('fuzhi');
subplot(2,1,2);F1=plot(f,abs(F0(1:512))); %画出滤波后的频谱图 title('滤波后的频谱');xlabel('Hz');ylabel('fuzhi'); 波形如下：
0.2
0.4
0.6
0.8
1
1.2
1.4
频率／Hz
频率响应幅度
Butterworth
-1-0.500.51滤波前的时域波形
-1-0.500.51滤波后的时域波形
0102030
40滤波前的频谱
Hz
f u z h i
2000
4000
60008000
10000
12000
0102030
40滤波后的频谱
Hz
f u z h i
4．回放语音信号
在Matlab 中，函数sound 可以对声音进行回放。

其调用格式：sound(x ，fs ，bits)；在运行Matlab 程序时，可以听到回放的声音，经过比较，明显感觉滤波前后的声音有变化。

四、结语
应用MATLAB 进行语音信号的处理是与我们所学课程及专业紧密相连的，有着很强的实践性，正因为如此，在做这个课题的过程中我们也遇到了一些困难，但在我们小组同学的共同努力下克服了这些困难，达到了预期的效果。

刚开始，在进行WAV 格式语音信号处理的时候，程序始终现实如下错误提示： ??? Error using ==> wavread Error using ==> wavread
Data compression format (IMA ADPCM) is not supported.
我们在查阅了很多资料，在网上也查阅相关信息，花费了大量时间也没找出结果，最后发现在WAV 格式的语音文件有两种格式，即PCM 格式和IMA ADPCM 格式，而在MATLAB 中用wavread 函数进行语音处理时，并不能直接处理IMA ADPCM 格式的语音信号，经过格式转换之后（选择PCM 格式），我们运行出了正确的结果。

实验过程中，我们很清楚的感觉到初始语音信号和滤波输出后的语音信号在音色上有一定的差别，这说明了信号在处理、传输过程中有损耗。

另外，在设计滤波器的过程中我们采用的是双线性变化法，事实上，采用窗函数法也可以进行，因此，不管对于多么简单的课题，其实也是有很多东西可以发掘的，这需要我们在平时多积累，多思考，只有这样，才能取得更大的进步，才能学有所用，学有所长。

加矩形窗的短时能量函数： a=wavread('F:\WO.wav');
N=32;
for i=2:6
h=linspace(1,1, (i-1)*N);
%形成一个矩形窗，长度为N
En=conv(h,a.*a);
%求卷积得其短时能量函数En
subplot(6,1,i),plot(En);
if(i==2) legend('N=32');
elseif(i==3) legend('N=64');
elseif(i==4) legend('N=128');
elseif(i==5) legend('N=256');
elseif(i==6) legend('N=512');
end
end
加hamming窗的短时能量函数：
把h=linspace(1,1, (i-1)*N);
改为h1=hamming((i-1)*N);
加矩形窗的短时平均幅度：
a=wavread('F:\WO.wav');
subplot(6,1,1),plot(a);
N=32;
for i=2:6
h=linspace(1,1,(i-1)*N);
%形成一个矩形窗，长度为N
En=conv(h,abs(a));
%求卷积得其短时能量函数En
subplot(6,1,i),plot(En);
if(i==2) legend('N=32');
elseif(i==3) legend('N=64');
elseif(i==4) legend('N=128');
elseif(i==5) legend('N=256');
elseif(i==6) legend('N=512');
end
end
短时过零率：
a=wavread('F:\WO.wav');
n=length(a);
N=320;
subplot(3,1,1),plot(a);
h=linspace(1,1,N);%形成一个矩形窗，长度为N En=conv(h,a.*a);%求卷积得其短时能量函数En
for i=1:n-1
if a(i)>=0
b(i)= 1;
else
b(i) = -1;
end
if a(i+1)>=0
b(i+1)=1;
else
b(i+1)=-1;
end
w(i)=abs(b(i+1)-b(i));
end%求出每相邻两点符号的差值的绝对值
k=1;
j=0;
while (k+N-1)<n
Zm(k)=0;
for i=0:N-1;
Zm(k)=Zm(k)+w(k+i);
end
j=j+1;
k=k+160; %每次移动半个窗
end
for w=1:j
Q(w)=Zm(160*(w-1)+1)/640;%短时平均过零率end
subplot(3,1,3),plot(Q);。