通过matlab实现语音信号输入抽样量化转换.docx
《通过matlab实现语音信号输入抽样量化转换.docx》由会员分享，可在线阅读，更多相关《通过matlab实现语音信号输入抽样量化转换.docx（14页珍藏版）》请在冰豆网上搜索。
通过matlab实现语音信号输入抽样量化转换
公司内部档案编码：
[OPPTR-OPPT28-OPPTL98-OPPNN08]
通过matlab实现语音信号输入抽样量化转换
通过matlab实现
1、语言信号的抽样、量化（A/D）；
2、数字语言信号的存储；
3、一维随机噪声的产生；
4、IIR、FIR低通滤波器设计
IIR数字滤波器的设计方法a冲激响应不变法b双线性变换法
FIR数字滤波器设计Kaiser（凯塞窗）海明窗布莱克曼窗
5附上fr，和ideal_lp，
Fs=8000;
time=input（'输入录音时间（整数秒）：
'）;
fprintf（'按回车开始录音%ds',time）;
pause;
disp（'录音中......'）;
x=audiorecorder（Fs,8,1）;
recordblocking（x,time）;
disp（'录音结束\n'）;
fprintf（'按回车回放录音'）;
pause;
play（x）;
myrecording=getaudiodata（x）;
audiowrite（'lll.wav',myrecording,Fs）;
y=audioread（'lll.wav'）;
figure;
plot（y）;
title（'录音语音信号'）;
time=10;
N=Fs*time;
u=fft（y,N）;
figure;
subplot（2,2,1）;
plot（abs（u））;
title（'N点FFT幅度谱'）;
grid;
subplot（2,2,2）;
plot（y）;
title（'语音信号'）;
grid;
L=length（y）;
noise=0.05*randn（L,1）;
yn=y+noise;
audioplayer（yn,Fs）;
N=Fs*time;
m=fft（yn,N）;
subplot（2,2,3）;
plot（abs（m））;
title（'加噪N点FFT幅度谱'）;
grid;
subplot（2,2,4）;
plot（yn）;
title（'语音信号加噪声'）;
grid;
以上是基础部分
紧接双线性变换iir
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*pi*Fp/Ft;
ws=2*pi*Fs/Ft;
fp=2*Ft*tan（wp/2）;
fs=2*Fs*tan（ws/2）;
[n,Wn]=buttord（wp,ws,1,50,'s'）;
[Z,P,K]=buttap（n）;
[Bap,Aap]=zp2tf（Z,P,K）;
[b,a]=lp2lp（Bap,Aap,Wn）;
[bz,ba]=bilinear（b,a,1）;
[H,W]=freqz（bz,ba,64）;
W1=W/pi;
H1=abs（H）;
H2=20*log10（H1）;
figure;
subplot（2,1,1）;
stem（W1,H1）;
grid;
xlabel（'频率'）;
ylabel（'幅度'）;
subplot（2,1,2）;
stem（W1,H2）;
grid;
xlabel（'频率'）;ylabel（'幅度（db）'）;
n=length（yn）;
Y=fft（yn,n）;
z22=filter（bz,ba,yn）;
sound（z22）;
m22=fft（z22）;
figure;
subplot（2,2,1）;
plot（abs（Y））;
title（'双线性IIR滤波前信号频谱'）;
grid;
subplot（2,2,2）;
plot（abs（m22））;
title（'双线性IIR滤波后信号频谱'）;
grid;
subplot（2,2,3）;
plot（z22）;
title（'双线性IIR滤波前信号波形'）;
grid;
subplot（2,2,4）;
plot（y）;
title（'双线性IIR滤波后号波形'）;
grid;
凯赛窗fir
Ft=5000;
Fp=1000;
Fs=1200;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
DB=ws-wp;%计算过渡带宽度
Rs=49;
beta=0.5842*（Rs-21）^0.4+0.07886*（Rs-21）;%计算凯塞窗的控制参数
M=ceil（（Rs-8）/2.285/DB）;%计算凯塞窗所需阶数M
wc=（wp+ws）/2/pi;%计算理想低通滤波器通带截止频率
hn=fir1（M,wc,kaiser（M+1,beta））;%调用firl函数计算低通FIRDF的h（n）
figure;
plot（hn）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
title（'数字滤波器幅频响应|H（ejOmega）|'）
figure;
freqz（hn,1,512）;
yy=fftfilt（hn,yn）;
yy1=fft（yy,L）;
yyy=abs（yy1）;
figure;
subplot（2,2,1）;
plot（abs（m））;
title（'加噪N点FFT幅度谱'）;
grid;
subplot（2,2,2）;
plot（yn）;
title（'语音信号加噪声'）;
grid;subplot（2,2,3）;
plot（yy）;
title（'凯赛处理后的语音信号波形图'）;
grid;subplot（2,2,4）;
plot（yyy）;
title（'凯赛处理后的FFt频谱图'）;
grid;
到这里是个完整的程序，如果需要进行部分改变请从下面选取响应部分进行替换
冲击响应不变法iir
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*pi*Fp/Ft;
ws=2*pi*Fs/Ft;
T=1;
[n,Wn]=buttord（wp,ws,1,50,'s'）;%计算模拟滤波器阶数N和截止频率wc
[B,A]=butter（n,Wn,'s'）;%计算相应模拟滤波器系统函数
[bz,ba]=impinvar（B,A,T）;%用脉冲响应不变法将模拟转化成数字滤波器
figure;
freqz（bz,ba,512）;
n=length（yn）;
Y=fft（yn,n）;
z22=filter（bz,ba,yn）;
sound（z22）;
m22=fft（z22）;
figure;
subplot（2,2,1）;
plot（abs（Y））;
title（'冲击响应不变法IIR滤波前信号频谱'）;
grid;
subplot（2,2,2）;
plot（abs（m22））;
title（'冲击响应不变法IIR滤波后信号频谱'）;
grid;
subplot（2,2,3）;
plot（z22）;
title（'冲击响应不变法IIR滤波前信号波形'）;
grid;
subplot（2,2,4）;
plot（y）;
title（'冲击响应不变法IIR滤波后号波形'）;
grid;
海明窗fir
Ft=5000;
Fp=1000;
Fs=1200;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
tr_width=ws-wp;
disp（'海明窗设计低通滤波器参数:
'）;
M=ceil（6.6*pi/tr_width）+1;
disp（['滤波器的长度为',num2str（M）]）;
n=0:
M-1;
wc=（ws+wp）/2;%理想LPF的截止频率
hd=ideal_lp（wc,M）;
w_ham=（hamming（M））';
hn=hd.*w_ham;
[db,mag,pha,gfd,w]=fr（hn,1）;
delta_w=2*pi/1000;
Rp=-（min（db（1:
1:
wp/delta_w+1）））;%求出实际通带波动
disp（['实际带通波动为',num2str（Rp）]）;
As=-round（max（db（ws/delta_w+1:
1:
501）））;%求出最小阻带衰减
disp（['最小阻带衰减为',num2str（As）]）;
figure;
subplot（2,2,1）;
plot（n,hd）;
title（'理想冲击响应'）;
axis（[0,M-1,-0.1,0.3]）;
ylabel（'hd（n）'）;
subplot（2,2,2）;
plot（n,w_ham）;
title（'海明窗'）;
axis（[0,M-1,0,1.1]）;
ylabel（'w（n）'）;
subplot（2,2,3）;
plot（n,hn）;
title（'实际冲激响应'）;
axis（[0,M-1,-0.1,0.3]）;
xlabel（'n'）;
ylabel（'h（n）'）;
subplot（2,2,4）;
plot（w/pi,db）;
title（'幅度响应（db）'）;
axis（[0,1,-100,10]）;
grid;
xlabel（'以pi为单位的频率'）;
ylabel（'分贝数'）;
yy=fftfilt（hn,yn）;
yy1=fft（yy,L）;
yyy=abs（yy1）;
figure;
subplot（2,2,1）;
plot（abs（m））;
title（'加噪N点FFT幅度谱'）;
grid;
subplot（2,2,2）;
plot（yn）;
title（'语音信号加噪声'）;
grid;subplot（2,2,3）;
plot（yy）;
title（'海明窗处理后的语音信号波形图'）;
grid;subplot（2,2,4）;
plot（yyy）;
title（'海明窗处理后的FFt频谱图'）;
grid;
布莱克曼窗fir
Ft=5000;
Fp=1000;
Fs=1200;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
r_width=ws-wp;
disp（'布莱克曼窗设计低通滤波器的参数:
'）;
M=ceil（11.0*pi/tr_width）+1;
disp（['滤波器的长度为',num2str（M）]）;
n=0:
M-1;%理想LPF的截止频率?
wc=（ws+wp）/2;
hd=ideal_lp（wc,M）;
w_bla=（blackman（M））';
hn=hd.*w_bla;
[db,mag,pha,gfd,w]=fr（hn,1）;
delta_w=2*pi/1000;
Rp=-（min（db（1:
1:
wp/delta_w+1）））;%求出实际通带波动
disp（['实际带通波动为',num2str（Rp）]）;
As=-round（max（db（ws/delta_w+1:
1:
501）））;%求出最小阻带衰减
disp（['最小阻带衰减db',num2str（As）]）;
figure;
subplot（2,2,1）;
plot（n,hd）;
title（'理想冲击响应'）;
axis（[0,M-1,-0.1,0.3]）;
ylabel（'hd（n）'）;
subplot（2,2,2）;
plot（n,w_bla）;
title（'布莱克曼窗'）;
axis（[0,M-1,0,1.1]）;
ylabel（'w（n）'）;
subplot（2,2,3）
plot（n,hn）;
title（'实际冲激响应'）;
axis（[0,M-1,-0.1,0.3]）;
xlabel（'n'）;
ylabel（'hn（n）'）;
subplot（2,2,4）;
plot（w/pi,db）;
title（'幅度响应（db）'）;
axis（[0,1,-100,10]）;
grid;
xlabel（'以pi为单位的频率'）;
ylabel（'分贝数'）;
yy=fftfilt（hn,yn）;
yy1=fft（yy,L）;
yyy=abs（yy1）;
figure;
subplot（2,2,1）;
plot（abs（m））;
title（'加噪N点FFT幅度谱'）;
grid;
subplot（2,2,2）;
plot（yn）;
title（'语音信号加噪声'）;
grid;subplot（2,2,3）;
plot（yy）;
title（'布莱克曼处理后的语音信号波形图'）;
grid;subplot（2,2,4）;
plot（yyy）;
title（'布莱克曼处理后的FFt频谱图'）;
grid;
以下2个需要单独保存
function[db,mag,pha,gfd,w]=fr（b,a）%N阶差分方程所描述的系统频响函数的m函数文件fr.m
%求解系统响应
%db为相位振幅（db）
%mag为绝对振幅
%pha为相位响应
%grd为群延时
%w为频率样本点矢量
%b为Ha（z）分析多项式系数（对FIR而言，b=h）
%a为Hz（z）分母多项式系数（对FIR而言，a=1）
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
501））';
w=（w（1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
gfd=grpdelay（b,a,w）;
functionhd=ideal_lp（wc,M）%理想低通滤波器单位抽样响应hd（n）的m函数文件ideal.m
%IdealLowPassfiltercomputation
%--------------------------------
%[hd]=ideal_lp（wc,M）;
%hd=idealimpulseresponsebetween0toM-1
%wc=cutofffrequencyinradians
%M=lengthoftheidealfilter
%
alpha=（M-1）/2;
n=（0:
1:
（M-1））;
m=n-alpha+eps;%addsmallestnumbertoavoidividedbyzero
hd=sin（wc*m）./（pi*m）;
升级会员 