MATLAB语音信号采集与处理DOCWord文件下载.docx
《MATLAB语音信号采集与处理DOCWord文件下载.docx》由会员分享，可在线阅读，更多相关《MATLAB语音信号采集与处理DOCWord文件下载.docx（19页珍藏版）》请在冰豆网上搜索。
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））;
原始信号频谱图'
）
仿真波形：
门铃：
和弦：
男女声：
2、调制与解调：
首先画出语音信号的时域波形，然后对语音信号进行频谱分析。
在Matlab
中可以利用函数fft
对信号行快速傅里叶变换,得到信号的频谱特性，从而加深对频谱特性的理解。
clear;
dt=1/44100;
fs=44100;
[f1,fs,nbits]=wavread（'
\1huan.wav'
figure
（1）;
subplot（1,1,1）;
N=length（f1）;
（N-1）/fs;
plot（t,f1）;
信息信号的时域波形'
fy1=fft（f1）;
w1=0:
fs/（N-1）:
fs;
figure
（2）;
plot（w1,abs（fy1））;
信息信号的频谱'
f2=cos（22000*pi*t）;
figure（3）;
fy2=fft（f2）;
N2=length（f2）;
w2=fs/N*[0:
N-1];
plot（w2,abs（abs（fy2）））;
载波信号的频谱'
f1=f1（:
1）;
f3=f1'
.*f2;
figure（4）;
fy3=fft（f3）;
plot（w1,abs（abs（fy3）））;
已调信号的频谱'
sound（f3,fs,nbits）;
f4=f3.*f2;
figure（5）;
fy4=fft（f4）;
plot（w1,abs（abs（fy4）））;
解调信号的频谱'
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）;
freqz（b1,1,512）;
FIR低通滤波器的频率响应'
f4_low=filter（b1,1,f4）;
plot（t,f4_low）;
滤波后的解调信号时域波形'
sound（f4_low,fs,nbits）;
f5=fft（f4_low）;
figure（7）;
plot（w1,abs（f5））;
滤波后的解调信号频谱'
3、信号变化：
快放：
[x,fs,nbits]=wavread（'
w=2；
M=w*fs;
wavplay（x,M）;
慢放：
w=0.8
倒放：
[x,fs,Nbits]=wavread（'
y0=flipud（x）;
sound（y0,fs）;
回声：
[x,fs,bits]=wavread（'
\3yang.wav'
[140000]）;
%读取语音信号
n1=0:
2000;
b=x（:
%产生单声道信号
N=3;
yy2=filter（1,[1,zeros（1,80000/（N+1））,0.7],[b'
zeros（1,40000）]）;
figure（3）
plot（yy2）;
%三次回声滤波器时域波形
三次回声滤波器时域波形'
YY2=fft（yy2）;
%对三次回声信号做FFT变换
plot（n1（1:
1000）,YY2（1:
1000））;
%三次回声滤波器频谱图
三次回声滤波器频谱图'
figure（4）
plot（abs（YY2））;
%经傅里叶变换之后的信号的幅值
幅值'
plot（angle（YY2））;
%经傅里叶变换之后的信号的相位
相位'
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）;
xabs=xpos+km;
Y（ypos+ovix）=（（1-xfwin）.*Y（ypos+ovix））+（xfwin.*padX（Wsim+xabs+ovix））;
Y（ypos+newix）=padX（Wsim+xabs+newix）;
end
变声程序：
[y,fs,nbits]=wavread（'
3yang.wav'
%读取声音文件
x=y（:
%读入的y矩阵有两列,取第1列
sound（voice（x,1.4）,fs,nbits）;
%调整voice（）第2个参数转换音调,>
1降调,<
1升调
4、信号加噪
fs=22050;
%语音信号采样频率为22050
x1=wavread（'
%读取语音信号数据赋值给x1
f=fs*（0:
511）/1024;
%将0到511，步长为1的序列的值与fs相乘并除以1024的值赋给f
（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）;
%创建图形窗
plot（t,x1）;
%做原始语音信号的时域波形
加噪前的信号'
xlabel（'
time
n'
%x轴的名字是time
ylabel（'
fuzhi
%y轴的名字是fuzhi
%创建两行一列绘图区间的第1个绘图区间
plot（t,x2）
加噪后的信号'
figure
（2）
plot（f,abs（y1（1:
512）））;
原始语音信号频谱'
Hz'
fuzhi'
plot（f,abs（y2（1:
加噪后的信号频谱'
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;
%
%按凯泽窗计算滤波器阶数
%求凯泽窗窗函数
%wc1/pi为归一化,窗函数法设计函数
%[H,w]=freqz（B,A,N）,
（1）中B和A分别为离散系统的系统函数分子、分母多项式的系数向量，返回量H则包含了离散系统频响在
0~pi范围内N个频率等分点的值（其中N为正整数），w则包含了范围内N个频率等分点。
调用默认的N时，其值是512。
x1_low
filter（b1,1,
x1）;
%对信号进行低通滤波
Y
filter（B,A,X）
，输入X为滤波前序列，Y为滤波结果序列，B/A
提供滤波器系数，B为分子，
A为分母
sound（x1_low,Fs,bits）;
plot（x1_low）;
信号经过FIR低通滤波器（时域）'
plot（[-Fs/2:
derta_Fs:
Fs/2-derta_Fs],abs（fftshift（fft（x1_low））））;
信号经过FIR低通滤波器（频域）'
FIR高通滤波器：
jiazaoyang.wav'
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）;
freqz（b2,1,512）;
%数字滤波器频率响应
FIR高通滤波器的频率响应'
x1_high
filter（b2,1,x1）;
%对信号进行高通滤波
sound（x1_high,Fs,bits）;
subplot（211）;
plot（x1_high）;
信号经过FIR高通滤波器（时域）'
subplot（212）;
Fs/2-derta_Fs],abs（fftshift（fft（x1_high））））;
信号经过FIR高通滤波器（频域）'
仿真结果：
FIR带通滤波：
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'
%带通滤波器
freqz（b3,1,512）;
FIR带通滤波器的频率响应'
x1_daitong
filter（b3,1,x1）;
%对信号进行带通滤波
sound（x1_daitong,Fs,bits）;
plot（x1_daitong）;
信号经过FIR带通滤波器（时域）'
Fs/2-derta_Fs],abs（fftshift（fft（x1_daitong））））;
信号经过FIR带通滤波器（频域）'
升级会员 