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

1、通过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);

2、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幅度谱');gr

3、id;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; 以上是基础部分紧接双线性变换iirFt=8000;Fp=1000;Fs=1200;wp=2*

4、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(&#

5、39;幅度');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;subp

6、lot(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;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;%计算理想低通滤波

7、器通带截止频率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;subpl

8、ot(2,2,2);plot(yn);title('语音信号加噪声');grid;subplot(2,2,3);plot(yy);title('凯赛处理后的语音信号波形图');grid;subplot(2,2,4);plot(yyy);title('凯赛处理后的FFt频谱图');grid;到这里是个完整的程序，如果需要进行部分改变 请从下面选取响应部分进行替换冲击响应不变法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,'

9、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(2,2,1);plot(abs(Y);title('冲击响应不变法IIR滤波前信号频谱');grid;subplot(2,2,2);plot(a

10、bs(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('海明窗设计低通滤波器参数:');M=ceil(6.6*pi/tr_width)+1;di

11、sp('滤波器的长度为',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('最小阻带衰减为&#

12、39;,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');

13、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('语音

14、信号加噪声');grid;subplot(2,2,3);plot(yy);title('海明窗处理后的语音信号波形图');grid;subplot(2,2,4);plot(yyy);title('海明窗处理后的FFt频谱图');grid;布莱克曼窗firFt=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)

15、;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;subpl

16、ot(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)');sub

17、plot(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

18、,2,3);plot(yy);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=(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函数文件idea