通过matlab实现语音信号输入抽样量化转换.docx

1、通过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(录音

2、中.); 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;subpl

3、ot(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; 以上是基础部分紧接双线性变换iirFt=8000;Fp=1000;Fs=1200;wp=2*pi*Fp/Ft;ws=2*pi*Fs/Ft;fp=2*Ft*tan(wp/

4、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)

5、; 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; 凯赛窗firFt=5000

6、;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);%计算凯塞窗所需阶数Mwc=(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)|)

7、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;到这里是个完整的程序，如果需要进行部分改变 请从下面选取响应

8、部分进行替换冲击响应不变法iirFt=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和截止频率wcB,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

9、(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;海明窗firFt=5000;Fp=1000;Fs=1200;wp=2*Fp/Ft;ws=2*Fs/Ft;tr_width=ws-wp;disp(海明窗设计低通滤波器参数

10、:);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(最小阻带衰减为

11、,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

12、,-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;布莱克曼窗fi

13、rFt=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);%求出实际通带波

14、动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

15、,-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

16、);title(布莱克曼处理后的语音信号波形图);grid;subplot(2,2,4);plot(yyy);title(布莱克曼处理后的FFt频谱图);grid; 以下2个需要单独保存functiondb,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

17、=(w(1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);gfd=grpdelay(b,a,w);function hd=ideal_lp(wc,M)%理想低通滤波器单位抽样响应hd(n)的m函数文件ideal.m% Ideal LowPass filter computation% -% hd = ideal_lp(wc,M);% hd = ideal impulse response between 0 to M-1% wc = cutoff frequency in radians% M = length of the ideal filter%alpha = (M-1)/2;n = (0:1:(M-1);m = n - alpha +eps; % add smallest number to avoi divided by zerohd = sin(wc*m)./(pi*m);