MATLAB语音信号采集及处理.docx
《MATLAB语音信号采集及处理.docx》由会员分享，可在线阅读，更多相关《MATLAB语音信号采集及处理.docx（20页珍藏版）》请在冰点文库上搜索。
MATLAB语音信号采集及处理
MATLAB课程设计报告
课题：
语音信号采集与处理
一、实践目的
本次课程设计的课题为《基于MATLAB的语音信号采集与处理》，学会运用MATLAB的信号处理功能，采集语音信号，并对语音信号进行滤波及变换处理，观察其时域和频域特性，加深对信号处理理论的理解，并为今后熟练使用MATLAB进行系统的分析仿真和设计奠定基础。
此次实习课程主要是为了进一步熟悉对matlab软件的使用，以及学会利用matlab对声音信号这种实际问题进行处理，将理论应用于实际，加深对它的理解。
二、实践原理：
利用MATLAB对语音信号进行分析和处理，采集语音信号后，利用MATLAB软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。
语音信号的“短时谱”对于非平稳信号,它是非周期的,频谱随时间连续变化,因此由傅里叶变换得到的频谱无法获知其在各个时刻的频谱特性。
如果利用加窗的方法从语音流中取出其中一个短断,再进行傅里叶变换,就可以得到该语音的短时谱。
三、课题要求：
利用windows自带的录音机或者其它录音软件，录制几段语音信号（要有几种不同的声音，要有男声、女声）。
对录制的语音信号进行频谱分析，确定该段语音的主要频率范围，由此频率范围判断该段语音信号的特点（低沉or尖锐）。
利用采样定理，对该段语音信号进行采样，观察不同采样频率（过采样、欠采样、临界采样）对信号的影响。
对采集到的语音信号进行调制与解调，观测调制与解调前后信号的变化。
实现语音信号的快放、慢放、倒放、回声、男女变声。
对语音信号加噪，然后进行滤波，分析不同的滤波方式对信号的影响。
利用MATLABGUI制作语音信号采集与分析演示系统。
四、MATLAB仿真
1、频谱分析：
用WINDOWS下的录音机，用单声道录制一段音乐或声音，时间在5S内。
然后MATLAB软件平台下，利用函数WAVREAD对语音信号进行采样，记住采样频率和采样点数。
对语音信号进行快速傅立叶变换，在一个窗口同时画出信号的时域波形图和频谱图，分析语音信号的频谱特点
程序：
fs=22050;
Nbits=16;
[x,fs,Nbits]=wavread（'D:
\matlab\22hexian.wav'）;%读声音文件
n=length（x）;
t=0:
1/fs:
（length（x）-1）/fs;%求出语音信号的长度
y1=fft（x,n）;%傅里叶变换
y2=fftshift（y1）;%对频谱图进行平移
f=0:
fs/n:
fs*（n-1）/n;%得出频点
subplot（2,1,1）;
plot（t/2,x）%做原始语音信号的时域图形
title（'原始信号时域波形图'）;
subplot（2,1,2）;
plot（f,abs（y2））;
title（'原始信号频谱图'）
仿真波形：
门铃：
和弦：
男女声：
2、调制与解调：
首先画出语音信号的时域波形，然后对语音信号进行频谱分析。
在Matlab 中可以利用函数fft 对信号行快速傅里叶变换,得到信号的频谱特性，从而加深对频谱特性的理解。
程序：
clear;
dt=1/44100;
fs=44100;
[f1,fs,nbits]=wavread（'D:
\1huan.wav'）;
figure
（1）;
subplot（1,1,1）;
N=length（f1）;
t=0:
1/fs:
（N-1）/fs;
plot（t,f1）;
title（'信息信号的时域波形'）;
fy1=fft（f1）;
w1=0:
fs/（N-1）:
fs;
figure
（2）;
subplot（1,1,1）;
plot（w1,abs（fy1））;
title（'信息信号的频谱'）;
f2=cos（22000*pi*t）;
figure（3）;
subplot（1,1,1）;
fy2=fft（f2）;
N2=length（f2）;
w2=fs/N*[0:
N-1];
plot（w2,abs（abs（fy2）））;
title（'载波信号的频谱'）;
f1=f1（:
1）;
f3=f1'.*f2;
figure（4）;
subplot（1,1,1）;
fy3=fft（f3）;
plot（w1,abs（abs（fy3）））;
title（'已调信号的频谱'）;
sound（f3,fs,nbits）;
f4=f3.*f2;
figure（5）;
subplot（1,1,1）;
fy4=fft（f4）;
plot（w1,abs（abs（fy4）））;
title（'解调信号的频谱'）;
sound（f4,fs,nbits）;
fp1=0;
fs1=5000;
As1=100;
wp1=2*pi*fp1/fs;
ws1=2*pi*fs1/fs;
BF1=ws1-wp1;
wc1=（wp1+ws1）/2;
M1=ceil（（As1-7.95）/（2.286*BF1））+1;
N1=M1+1;
beta1=0.1102*（As1-8.7）;
Window=（kaiser（N1,beta1））;
b1=fir1（M1,wc1/pi,Window）;
figure（6）;
subplot（1,1,1）;
freqz（b1,1,512）;
title（'FIR低通滤波器的频率响应'）;
f4_low=filter（b1,1,f4）;
plot（t,f4_low）;
title（'滤波后的解调信号时域波形'）;
sound（f4_low,fs,nbits）;
f5=fft（f4_low）;
figure（7）;
subplot（1,1,1）;
plot（w1,abs（f5））;
title（'滤波后的解调信号频谱'）;
仿真波形：
3、信号变化：
快放：
[x,fs,nbits]=wavread（'D:
\1huan.wav'）
w=2；
M=w*fs;
wavplay（x,M）;
慢放：
[x,fs,nbits]=wavread（'D:
\1huan.wav'）
w=0.8
M=w*fs;
wavplay（x,M）;
倒放：
[x,fs,Nbits]=wavread（'D:
\1huan.wav'）;
y0=flipud（x）;
sound（y0,fs）;
回声：
程序：
[x,fs,bits]=wavread（'D:
\3yang.wav',[140000]）;%读取语音信号
n1=0:
2000;
b=x（:
1）;%产生单声道信号
N=3;
yy2=filter（1,[1,zeros（1,80000/（N+1））,0.7],[b',zeros（1,40000）]）;
figure（3）
subplot（2,1,1）;plot（yy2）;%三次回声滤波器时域波形
title（'三次回声滤波器时域波形'）;
YY2=fft（yy2）;%对三次回声信号做FFT变换
subplot（2,1,2）;plot（n1（1:
1000）,YY2（1:
1000））;%三次回声滤波器频谱图
title（'三次回声滤波器频谱图'）;
figure（4）
subplot（2,1,1）;plot（abs（YY2））;%经傅里叶变换之后的信号的幅值
title（'幅值'）;
subplot（2,1,2）;plot（angle（YY2））;%经傅里叶变换之后的信号的相位
title（'相位'）;
sound（2*yy2,fs,bits）;%经三次回声滤波器后的语音信号，乘以2是为了加强信号
仿真波形：
男女变声：
Voice调用函数：
functionY=voice（x,f）%更改采样率使基频改变f>1降低;f<1升高
f=round（f*1000）;
d=resample（x,f,1000）;%时长整合使语音文件恢复原来时长
W=400;
Wov=W/2;
Kmax=W*2;
Wsim=Wov;
xdecim=8;
kdecim=2;X=d';
F=f/1000;
Ss=W-Wov;
xpts=size（X,2）;
ypts=round（xpts/F）;
Y=zeros（1,ypts）;
xfwin=（1:
Wov）/（Wov+1）;
ovix=（1-Wov）:
0;newix=1:
（W-Wov）;
simix=（1:
xdecim:
Wsim）-Wsim;
padX=[zeros（1,Wsim）,X,zeros（1,Kmax+W-Wov）];
Y（1:
Wsim）=X（1:
Wsim）;lastxpos=0;km=0;
forypos=Wsim:
Ss:
（ypts-W）
xpos=round（F*ypos）;
kmpred=km+（xpos-lastxpos）;
lastxpos=xpos;
if（kmpred<=Kmax）
km=kmpred;
else
ysim=Y（ypos+simix）;
rxy=zeros（1,Kmax+1）;
rxx=zeros（1,Kmax+1）;
Kmin=0;
fork=Kmin:
kdecim:
Kmax
xsim=padX（Wsim+xpos+k+simix）;
rxx（k+1）=norm（xsim）;
rxy（k+1）=（ysim*xsim'）;
end
Rxy=（rxx~=0）.*rxy./（rxx+（rxx==0））;
km=min（find（Rxy==max（Rxy））-1）;
end
xabs=xpos+km;
Y（ypos+ovix）=（（1-xfwin）.*Y（ypos+ovix））+（xfwin.*padX（Wsim+xabs+ovix））;
Y（ypos+newix）=padX（Wsim+xabs+newix）;
end
end
变声程序：
[y,fs,nbits]=wavread（'3yang.wav'）;%读取声音文件
x=y（:
1）;%读入的y矩阵有两列,取第1列
sound（voice（x,1.4）,fs,nbits）;%调整voice（）第2个参数转换音调,>1降调,<1升调
4、信号加噪
程序：
fs=22050; %语音信号采样频率为22050
x1=wavread（'3yang.wav'）;%读取语音信号数据赋值给x1
f=fs*（0:
511）/1024; %将0到511，步长为1的序列的值与fs相乘并除以1024的值赋给f
t=0:
1/fs:
（length（x1）-1）/fs;%将0到x1的长度减1后的值除以fs的值，且步长为1/fs的值，的序列的值，赋予t
Au=0.05;%噪声幅值
d=[Au*cos（2*pi*2000*t）]'; %干扰信号构建命令函数，构建了一个余弦函数
x2=x1（:
1）+d;%取原始语音信号的单声部信号然后和噪声信号相加
wavwrite（x2,22050,'jiazaoyang'）;%生成wav文件
sound（x2,22050）; %播放语音信号
y1=fft（x1,1024）; %对信号做1024点的FFT变换
y2=fft（x2,1024）;
figure
（1）; %创建图形窗
subplot（2,1,1）;
plot（t,x1）; %做原始语音信号的时域波形
title（'加噪前的信号'）;
xlabel（'time n'）; %x轴的名字是time
ylabel（'fuzhi n'）; %y轴的名字是fuzhi
subplot（2,1,2）;%创建两行一列绘图区间的第1个绘图区间
plot（t,x2）
title（'加噪后的信号'）;
xlabel（'time n'）;
ylabel（'fuzhi n'）;
figure
（2）
subplot（2,1,1）;
plot（f,abs（y1（1:
512）））;
title（'原始语音信号频谱'）;
xlabel（'Hz'）;
ylabel（'fuzhi'）;
subplot（2,1,2）;
plot（f,abs（y2（1:
512）））;
title（'加噪后的信号频谱'）;
xlabel（'Hz'）;
ylabel（'fuzhi'）;
仿真波形：
5、用窗函数法设计FIR滤波器
根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为Δw，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。
根据待求滤波器的理想频率响应求出理想单位脉冲响应hd（n），如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：
在一般情况下，hd（n）是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2π采样N点，采用离散傅里叶反变换（IDFT）即可求出。
用窗函数wd（n）将hd（n）截断，并进行加权处理，得到
如果要求线性相位特性，则h（n）还必须满足：
根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。
要根据所设计的滤波特性正确选择其中一类。
例如，要设计线性相位低通特性可选择h（n）=h（N-1-n）一类，而不能选h（n）=-h（N-1-n）一类。
验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。
FIR低通滤波器：
程序：
[x1,Fs,bits]=wavread（'jiazapyang'）;
derta_Fs = Fs/length（x1）;%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致
fs=Fs;
fp1=1000;
fs1=1200;
As1=100;
wp1=2*pi*fp1/fs; %
ws1=2*pi*fs1/fs; %
BF1=ws1-wp1;
wc1=（wp1+ws1）/2;
M1=ceil（（As1-7.95）/（2.286*BF1））+1;%按凯泽窗计算滤波器阶数
N1=M1+1;
beta1=0.1102*（As1-8.7）;
Window=（kaiser（N1,beta1））; %求凯泽窗窗函数
b1=fir1（M1,wc1/pi,Window）;%wc1/pi为归一化,窗函数法设计函数
figure
（2）;
freqz（b1,1,512）; %[H,w]=freqz（B,A,N）,
（1）中B和A分别为离散系统的系统函数分子、分母多项式的系数向量，返回量H则包含了离散系统频响在 0~pi范围内N个频率等分点的值（其中N为正整数），w则包含了范围内N个频率等分点。
调用默认的N时，其值是512。
title（'FIR低通滤波器的频率响应'）;
x1_low = filter（b1,1, x1）;%对信号进行低通滤波 ,Y = filter（B,A,X） ，输入X为滤波前序列，Y为滤波结果序列，B/A 提供滤波器系数，B为分子， A为分母
sound（x1_low,Fs,bits）;
figure（3）;
subplot（2,1,1）;
plot（x1_low）;
title（'信号经过FIR低通滤波器（时域）'）;
subplot（2,1,2）;
plot（[-Fs/2:
derta_Fs:
Fs/2-derta_Fs],abs（fftshift（fft（x1_low））））;
title（'信号经过FIR低通滤波器（频域）'）;
仿真波形：
FIR高通滤波器：
程序：
[x1,Fs,bits]=wavread（'jiazaoyang.wav'）;
derta_Fs = Fs/length（x1）;%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致
fs=Fs;
As2=100;
fp2=4000;
fs2=2000;
wp2=2*pi*fp2/fs;
ws2=2*pi*fs2/fs;
BF2=wp2-ws2;
wc2=（wp2+ws2）/2;
M2=ceil（（As2-7.95）/（2.286*BF2））+1;%按凯泽窗计算滤波器阶数
N2=M2+1;
beta2=0.1102*（As2-8.7）;
Window=（kaiser（N2,beta2））; %求凯泽窗窗函数
b2=fir1（M2,wc2/pi,'high',Window）;
figure（4）;
freqz（b2,1,512）;%数字滤波器频率响应
title（'FIR高通滤波器的频率响应'）;
x1_high = filter（b2,1,x1）;%对信号进行高通滤波
sound（x1_high,Fs,bits）;
figure（5）;
subplot（211）;
plot（x1_high）;
title（'信号经过FIR高通滤波器（时域）'）;
subplot（212）;
plot（[-Fs/2:
derta_Fs:
Fs/2-derta_Fs],abs（fftshift（fft（x1_high））））;
title（'信号经过FIR高通滤波器（频域）'）;
仿真结果：
FIR带通滤波：
程序：
[x1,Fs,bits]=wavread（'jiazaoyang.wav'）;
derta_Fs = Fs/length（x1）;%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致
fs=Fs;
As3=100;
fp3=[1200,3000];fs3=[1000,3200];
wp3=2*pi*fp3/fs;
ws3=2*pi*fs3/fs;
BF3=wp3
（1）-ws3
（1）;
wc3=wp3+BF3/2;
M3=ceil（（As3-7.95）/（2.286*BF3））+1;%按凯泽窗计算滤波器阶数
N3=M3+1;
beta3=0.1102*（As3-8.7）;
Window=（kaiser（N3,beta3））; %求凯泽窗窗函数
b3=fir1（M3,wc3/pi,'bandpass',Window）;%带通滤波器
figure（6）;
freqz（b3,1,512）;%数字滤波器频率响应
title（'FIR带通滤波器的频率响应'）;
x1_daitong = filter（b3,1,x1）;%对信号进行带通滤波
sound（x1_daitong,Fs,bits）;
figure（7）;
subplot（211）;
plot（x1_daitong）;
title（'信号经过FIR带通滤波器（时域）'）;
subplot（212）;
plot（[-Fs/2:
derta_Fs:
Fs/2-derta_Fs],abs（fftshift（fft（x1_daitong））））;
title（'信号经过FIR带通滤波器（频域）'）;
仿真结果：
升级会员 