计算机处理方法

计算机处理方法第一页，共一百五十五页，编辑于2023年，星期二6.1双门限法第二页，共一百五十五页，编辑于2023年，星期二图6-1-1语音波形语音短时平均能量语音短时平均过零率第三页，共一百五十五页，编辑于2023年，星期二进行判决的具体步骤：第一级判决：第四页，共一百五十五页，编辑于2023年，星期二第2级判决第五页，共一百五十五页，编辑于2023年，星期二计算短时平均能量例：读入数据文件bluesky1.wav（内容为男声“蓝天，白云，碧绿的大海”），求其短时平均能量第六页，共一百五十五页，编辑于2023年，星期二filedir=[];%设置路径filename='bluesky1.wav';%设置文件名fle=[filedirfilename];%构成完整的路径和文件名[x,Fs]=wavread(fle);%读入数据文件wlen=200;inc=80;%给出帧长和帧移win=hanning(wlen);%给出海宁窗N=length(x);%信号长度X=enframe(x,win,inc)';%分帧fn=size(X,2);%求出帧数time=(0:N-1)/Fs;%计算出信号的时间刻度fori=1:fnu=X(:,i);%取出一帧

u2=u.*u;%求出能量

En(i)=sum(u2);%对一帧累加求和endsubplot211;plot(time,x,'k');%画出时间波形title('语音波形');ylabel('幅值');xlabel(['时间/s'10'(a)']);frameTime=frame2time(fn,wlen,inc,Fs);%求出每帧对应的时间subplot212;plot(frameTime,En,'k')%画出短时能量图title('短时能量');ylabel('幅值');xlabel(['时间/s'10'(b)']);第七页，共一百五十五页，编辑于2023年，星期二第八页，共一百五十五页，编辑于2023年，星期二例：读入数据文件bluesky1.wav（内容为男声“蓝天，白云，碧绿的大海”），求其短时平均过零率filedir=[];%设置路径filename='bluesky1.wav';%设置文件名fle=[filedirfilename];%构成完整的路径和文件名[xx,Fs]=wavread(fle);%读入数据文件x=xx-mean(xx);%消除直流分量wlen=200;inc=80;%设置帧长、帧移win=hanning(wlen);%窗函数N=length(x);%求数据长度X=enframe(x,win,inc)';%分帧fn=size(X,2);%获取帧数zcr1=zeros(1,fn);%初始化fori=1:fnz=X(:,i);%取得一帧数据

forj=1:(wlen-1);%在一帧内寻找过零点

ifz(j)*z(j+1)<0%判断是否为过零点

zcr1(i)=zcr1(i)+1;%是过零点，记录1次

endendendtime=(0:N-1)/Fs;%计算时间坐标frameTime=frame2time(fn,wlen,inc,Fs);%求出每帧对应的时间%作图subplot211;plot(time,x,'k');grid;title('语音波形');ylabel('幅值');xlabel(['时间/s'10'(a)']);subplot212;plot(frameTime,zcr1,'k');grid;title('短时平均过零率');ylabel('幅值');xlabel(['时间/s'10'(b)']);第九页，共一百五十五页，编辑于2023年，星期二第十页，共一百五十五页，编辑于2023年，星期二用短时平均能量和短时过零率提取语音端点位置第十一页，共一百五十五页，编辑于2023年，星期二第十二页，共一百五十五页，编辑于2023年，星期二第十三页，共一百五十五页，编辑于2023年，星期二第十四页，共一百五十五页，编辑于2023年，星期二%pr6_1_1clearall;clc;closeall;filedir=[];%指定文件路径filename='bluesky1.wav';%指定文件名fle=[filedirfilename]%构成路径和文件名的字符串[x,fs]=wavread(fle);%读入数据文件x=x/max(abs(x));%幅度归一化N=length(x);%取信号长度time=(0:N-1)/fs;%计算时间pos=get(gcf,'Position');%作图set(gcf,'Position',[pos(1),pos(2)-100,pos(3),(pos(4)-200)]);plot(time,x,'k');title('男声“蓝天，白云，碧绿的大海”的端点检测');ylabel('幅值');axis([0max(time)-11]);grid;xlabel('时间/s');设置无话段长度IS：0.1s第十五页，共一百五十五页，编辑于2023年，星期二wlen=200;inc=80;%分帧参数IS=0.1;overlap=wlen-inc;%设置ISNIS=fix((IS*fs-wlen)/inc+1);%计算NISfn=fix((N-wlen)/inc)+1;%求帧数frameTime=frame2time(fn,wlen,inc,fs);%计算每帧对应的时间[voiceseg,vsl,SF,NF]=vad_ezm1(x,wlen,inc,NIS);%端点检测fork=1:vsl%画出起止点位置

nx1=voiceseg(k).begin;nx2=voiceseg(k).end;nxl=voiceseg(k).duration;fprintf('%4d%4d%4d%4d\n',k,nx1,nx2,nxl);line([frameTime(nx1)frameTime(nx1)],[-1.51.5],'color','k','LineStyle','-');line([frameTime(nx2)frameTime(nx2)],[-1.51.5],'color','k','LineStyle','--');end第十六页，共一百五十五页，编辑于2023年，星期二voiceseg=1x3structarraywithfields:beginendduration>>vslvsl=3第十七页，共一百五十五页，编辑于2023年，星期二>>SF(1:10)ans=0000000000>>SF(25:34)ans=0000111111有话帧标记：1：有话帧，0：无话帧第十八页，共一百五十五页，编辑于2023年，星期二>>NF(1:10)ans=1111111111>>NF(25:34)ans=1111000000无话帧标记：1：无话帧，0：有话帧第十九页，共一百五十五页，编辑于2023年，星期二fle=bluesky1.wav129123952131222923245377133有话段起始点（第29帧）有话段终止点（第123帧）共3段有话帧有话帧长度第二十页，共一百五十五页，编辑于2023年，星期二第二十一页，共一百五十五页，编辑于2023年，星期二课堂练习设置无话段长度IS：0.2s,0.3s,0.4s,用vad_ezm1函数对bluesky1.wav（男声：“蓝天，白云，碧绿的大海”）数据进行端点检测，并比较结果第二十二页，共一百五十五页，编辑于2023年，星期二6.2双门限法的改进和推广第二十三页，共一百五十五页，编辑于2023年，星期二第二十四页，共一百五十五页，编辑于2023年，星期二第二十五页，共一百五十五页，编辑于2023年，星期二第二十六页，共一百五十五页，编辑于2023年，星期二%pr6_2_1clearall;clc;closeall;filedir=[];%指定文件路径filename='bluesky1.wav';%指定文件名fle=[filedirfilename]%构成路径和文件名的字符串[xx,fs]=wavread(fle);%读入数据文件xx=xx/max(abs(xx));%幅度归一化N=length(xx);%取信号长度time=(0:N-1)/fs;%计算时间刻度x=Gnoisegen(xx,20);%把白噪声叠加到信号上wlen=200;inc=80;%设置帧长和帧移IS=0.25;overlap=wlen-inc;%设置前导无话段长度NIS=fix((IS*fs-wlen)/inc+1);%计算前导无话段帧数fn=fix((N-wlen)/inc)+1;%求出总帧数frameTime=frame2time(fn,wlen,inc,fs);%计算每帧对应的时间[voiceseg,vsl,SF,NF]=vad_ezr(x,wlen,inc,NIS);%端点检测第二十七页，共一百五十五页，编辑于2023年，星期二%作图subplot211;plot(time,xx,'k');holdontitle('纯语音男声“蓝天，白云，碧绿的大海”波形');ylabel('幅值');axis([0max(time)-11]);xlabel('(a)');fork=1:vslnx1=voiceseg(k).begin;nx2=voiceseg(k).end;fprintf('%4d%4d%4d\n',k,nx1,nx2);line([frameTime(nx1)frameTime(nx1)],[-1.51.5],'color','k','LineStyle','-');line([frameTime(nx2)frameTime(nx2)],[-1.51.5],'color','k','LineStyle','--');endsubplot212;plot(time,x,'k');title('加噪语音波形(信噪比20dB)');ylabel('幅值');axis([0max(time)-11]);xlabel(['时间/s'10'(b)']);第二十八页，共一百五十五页，编辑于2023年，星期二13013021322303245370第二十九页，共一百五十五页，编辑于2023年，星期二第三十页，共一百五十五页，编辑于2023年，星期二第三十一页，共一百五十五页，编辑于2023年，星期二课堂练习：（1）产生400点高斯分布随机信号并求它的幅频特性（2）产生400点均匀分布随机信号并求它的幅频特性第三十二页，共一百五十五页，编辑于2023年，星期二第三十三页，共一百五十五页，编辑于2023年，星期二第三十四页，共一百五十五页，编辑于2023年，星期二%%pr5_3_1clearall;clc;closeall;filedir=[];%指定文件路径filename='bluesky3.wav';%指定文件名fle=[filedirfilename];[s,fs]=wavread(fle);%读入数据文件s=s-mean(s);%消除直流分量s=s/max(abs(s));%幅值归一化N=length(s);%求出数据长度time=(0:N-1)/fs;%求出时间刻度subplot411;plot(time,s,'k');%画出纯语音信号的波形图title('纯语音信号');ylabel('幅值')第三十五页，共一百五十五页，编辑于2023年，星期二SNR=[1550];%信噪比的取值区间fork=1:3snr=SNR(k);%设定信噪比

[x,noise]=Gnoisegen(s,snr);%求出相应信噪比的高斯白噪声，构成带噪语音

subplot(4,1,k+1);plot(time,x,'k');ylabel('幅值');%作图

snr1=SNR_singlech(s,x);%计算出带噪语音中的信噪比

fprintf('k=%4dsnr=%5.1fsnr1=%5.4f\n',k,snr,round(snr1*1e4)/1e4);title(['带噪语音信号设定信噪比='num2str(snr)'dB计算出信噪比='...num2str(round(snr1*1e4)/1e4)'dB']);endxlabel('时间/s')第三十六页，共一百五十五页，编辑于2023年，星期二k=1snr=15.0snr1=15.0000k=2snr=5.0snr1=5.0000k=3snr=0.0snr1=0.0000第三十七页，共一百五十五页，编辑于2023年，星期二加15，5，0dB的高斯白噪声的波形图第三十八页，共一百五十五页，编辑于2023年，星期二课堂练习读入数据文件redriver2.wav

，调用Gnoisegen函数分别对语音信号加15，5，0dB的高斯白噪声，并用SNR_singlech函数计算出信噪比第三十九页，共一百五十五页，编辑于2023年，星期二例5-3-11(pr5_3_11)读入数据文件bluesky3.wav（内容为男声“蓝天，白云”），调用Nnoisegen函数分别对语音信号加15，5，0dB的均匀白噪声，并用SNR_singlech函数计算出信噪比%%pr5_3_11clearall;clc;closeall;filedir=[];%指定文件路径filename='bluesky3.wav';%指定文件名fle=[filedirfilename];[s,fs]=wavread(fle);%读入数据文件s=s-mean(s);%消除直流分量s=s/max(abs(s));%幅值归一化N=length(s);%求出数据长度time=(0:N-1)/fs;%求出时间刻度subplot411;plot(time,s,'k');%画出纯语音信号的波形图title('纯语音信号');ylabel('幅值')第四十页，共一百五十五页，编辑于2023年，星期二SNR=[1550];%信噪比的取值区间fork=1:3snr=SNR(k);%设定信噪比

[x,noise]=Nnoisegen(s,snr);%求出相应信噪比的高斯白噪声，构成带噪语音

subplot(4,1,k+1);plot(time,x,'k');ylabel('幅值');%作图

snr1=SNR_singlech(s,x);%计算出带噪语音中的信噪比

fprintf('k=%4dsnr=%5.1fsnr1=%5.4f\n',k,snr,round(snr1*1e4)/1e4);title(['带噪语音信号设定信噪比='num2str(snr)'dB计算出信噪比='...num2str(round(snr1*1e4)/1e4)'dB']);endxlabel('时间/s')第四十一页，共一百五十五页，编辑于2023年，星期二k=1snr=15.0snr1=15.0000k=2snr=5.0snr1=5.0000k=3snr=0.0snr1=0.0000第四十二页，共一百五十五页，编辑于2023年，星期二加15，5，0dB的均匀白噪声的波形图第四十三页，共一百五十五页，编辑于2023年，星期二课堂练习读入数据文件redriver2.wav

，调用Nnoisegen函数分别对语音信号加15，5，0dB的均匀白噪声，并用SNR_singlech函数计算出信噪比第四十四页，共一百五十五页，编辑于2023年，星期二(3)生成叠加任意噪声数据的带噪语音第四十五页，共一百五十五页，编辑于2023年，星期二Add_noisedata函数清单（1）function[signal,noise]=add_noisedata(s,data,fs,fs1,snr)s=s(:);%把信号转换成列数据s=s-mean(s);%消除直流分量sL=length(s);%求出信号的长度iffs~=fs1%若纯语音信号的采样频率与噪声的采样频率不相等

x=resample(data,fs,fs1);%对噪声重采样，使噪声采样频率与纯语音信号的采样频率相同elsex=data;Endx=x(:);%把噪声数据转换成列数据x=x-mean(x);%消除直流分量xL=length(x);%求噪声数据长度ifxL>=sL%如果噪声数据长度与信号数据长度不等，把噪声数据截断或补足

x=x(1:sL);elsedisp('Warning:噪声数据短于信号数据，以补0来补足！')x=[x;zeros(sL-xL,1)];end第四十六页，共一百五十五页，编辑于2023年，星期二Add_noisedata函数清单（2）Sr=snr;Es=sum(x.*x);%求出信号的能量Ev=sum(s.*s);%求出噪声的能量a=sqrt(Ev/Es/(10^(Sr/10)));%计算出噪声的比例因子noise=a*x;%调整噪声的幅值signal=s+noise;%构成带噪语音第四十七页，共一百五十五页，编辑于2023年，星期二%pr5_3_2clearall;clc;closeall;filedir=[];%指定文件路径filename='bluesky3.wav';%指定文件名fle=[filedirfilename];[s,fs]=wavread(fle);%读入数据文件s=s/max(abs(s));%幅值归一化N=length(s);%求出数据长度time=(0:N-1)/fs;%求出时间刻度subplot411;plot(time,s,'k');%画出纯语音信号的波形图title('纯语音信号');ylabel('幅值')采样频率8000Hz100第四十八页，共一百五十五页，编辑于2023年，星期二SNR=[50-5];%信噪比的取值区间fork=1:3snr=SNR(k);%设定信噪比

data=sin(2*pi*100*time);%产生一个正弦信号

[x,noise]=add_noisedata(s,data,fs,fs,snr);%按信噪比构成正弦信号叠加到语音上

subplot(4,1,k+1);plot(time,x,'k');ylabel('幅值');%作图

ylim([-22]);snr1=SNR_singlech(s,x);%计算出带噪语音中的信噪比

fprintf('k=%4dsnr=%5.1fsnr1=%5.4f\n',k,snr,snr1);title(['带噪语音信号设定信噪比='num2str(snr)'dB计算出信噪比='...num2str(round(snr1*1e4)/1e4)'dB']);endxlabel('时间/s')第四十九页，共一百五十五页，编辑于2023年，星期二k=1snr=5.0snr1=5.0000k=2snr=0.0snr1=-0.0000k=3snr=-5.0snr1=-5.0000第五十页，共一百五十五页，编辑于2023年，星期二第五十一页，共一百五十五页，编辑于2023年，星期二课堂练习读入数据文件redriver2.wav，叠加一个以不同信噪比5，0，-5dB的正弦波，正弦波频率为10Hz，采样频率为8000Hz。第五十二页，共一百五十五页，编辑于2023年，星期二（4）把噪声数据文件叠加生成任意信噪比的带噪语音第五十三页，共一百五十五页，编辑于2023年，星期二add_noisefile函数清单（1）function[y,noise]=add_noisefile(s,filepath_name,SNR,fs)s=s(:);%把信号转换成列数据s=s-mean(s);%消除直流分量[wavin,fs1,nbits]=wavread(filepath_name);%读入噪声文件的数据wavin=wavin(:);%把噪声数据转换成列数据iffs1~=fs%纯语音信号的采样频率与噪声的采样频率不相等

wavin1=resample(wavin,fs,fs1);%对噪声重采样，使噪声采样频率与纯语音信号的采样频率相同elsewavin1=wavin;endwavin1=wavin1-mean(wavin1);%消除直流分量第五十四页，共一百五十五页，编辑于2023年，星期二add_noisefile函数清单（2）ns=length(s);%求出s的长度noise=wavin1(1:ns);%把噪声长度截断为与s等长noise=noise-mean(noise);%噪声去除直流分量signal_power=1/ns*sum(s.*s);%求出信号的能量noise_power=1/ns*sum(noise.*noise);%求出噪声的能量noise_variance=signal_power/(10^(SNR/10));%求出噪声设定的方差值noise=sqrt(noise_variance/noise_power)*noise;%调整噪声幅值y=s+noise;%构成带噪语音第五十五页，共一百五十五页，编辑于2023年，星期二%pr5_3_3clearall;clc;closeall;filedir=[];%指定文件路径filename='bluesky3.wav';%指定文件名fle=[filedirfilename];[s,fs]=wavread(fle);%读入数据文件s=s/max(abs(s));%幅值归一化N=length(s);%求出数据长度time=(0:N-1)/fs;%求出时间刻度subplot411;plot(time,s,'k');%画出纯语音信号的波形图title('纯语音信号');ylabel('幅值')filepath_name='factory1.wav';第五十六页，共一百五十五页，编辑于2023年，星期二SNR=[50-5];%信噪比的取值区间fork=1:3snr=SNR(k);%设定信噪比

[x,noise]=add_noisefile(s,filepath_name,snr,fs);%按信噪比构成噪声叠加到语音上

subplot(4,1,k+1);plot(time,x,'k');ylabel('幅值');%作图

ylim([-22]);snr1=SNR_singlech(s,x);%计算出带噪语音中的信噪比

fprintf('k=%4dsnr=%5.1fsnr1=%5.4f\n',k,snr,snr1);title(['带噪语音信号设定信噪比='num2str(snr)'dB计算出信噪比='...num2str(round(snr1*1e4)/1e4)'dB']);endxlabel('时间/s')

第五十七页，共一百五十五页，编辑于2023年，星期二例pr5_3_3在MATLAB命令栏中输出结果k=1snr=5.0snr1=5.0000k=2snr=0.0snr1=-0.0000k=3snr=-5.0snr1=-5.0000第五十八页，共一百五十五页，编辑于2023年，星期二例pr5_3_3波形第五十九页，共一百五十五页，编辑于2023年，星期二factory1.wav波形及频谱%pr5_3_3clearall;clc;closeall;filedir=[];%指定文件路径filename='factory1.wav';%指定文件名fle=[filedirfilename];[s,fs]=wavread(fle);%读入数据文件subplot(211)plot(s);subplot(212)S=fft(s,1024);plot(1:512,abs(S(1:512)));第六十页，共一百五十五页，编辑于2023年，星期二factory1.wav波形及频谱第六十一页，共一百五十五页，编辑于2023年，星期二课堂练习：读入数据文件redriver2.wav

，叠加噪声文件名factory1.wav噪声。设不同信噪比10，0，-10dB第六十二页，共一百五十五页，编辑于2023年，星期二5.4语音信号的预处理一—消除趋势项和直流分量第六十三页，共一百五十五页，编辑于2023年，星期二5.4.1最小二乘法拟合趋势项的原理第六十四页，共一百五十五页，编辑于2023年，星期二第六十五页，共一百五十五页，编辑于2023年，星期二消除常数趋势项：第六十六页，共一百五十五页，编辑于2023年，星期二消除线性趋势项：第六十七页，共一百五十五页，编辑于2023年，星期二消除线性趋势项：第六十八页，共一百五十五页，编辑于2023年，星期二5.4.2最小二乘法拟合消除趋势项的函数第六十九页，共一百五十五页，编辑于2023年，星期二程序清单function[y,xtrend]=polydetrend(x,fs,m)x=x(:);%把语音信号x转换为列数据N=length(x);%求出x的长度t=(0:N-1)'/fs;%按x的长度和采样频率设置时间序列a=polyfit(t,x,m);%用最小二乘法拟合语音信号x的多项式系数axtrend=polyval(a,t);%用系数a和时间序列t构成趋势项y=x-xtrend;%从语音信号x中清除趋势项第七十页，共一百五十五页，编辑于2023年，星期二%%pr5_4_1clearall;clc;closeall;[x,fs,nbit]=wavread('bluesky31.wav');%读入bluesky31.wav文件t=(0:length(x)-1)/fs;%设置时间y=detrend(x);%消除线性趋势项y=y/max(abs(y));%幅值归一化subplot211;plot(t,x,'k');%画出带有趋势项的语音信号xtitle('带趋势项的语音信号');xlabel('时间/s');ylabel('幅值');subplot212;plot(t,y,'k');%画出消除趋势项的语音信号yxlabel('时间/s');ylabel('幅值');title('消除趋势项的语音信号');第七十一页，共一百五十五页，编辑于2023年，星期二第七十二页，共一百五十五页，编辑于2023年，星期二%pr5_4_2clearall;clc;closeall;[x,fs,nbit]=wavread('bluesky32.wav');%读入bluesky32.wav文件[y,xtrend]=polydetrend(x,fs,2);%调用polydetrend消除趋势项t=(0:length(x)-1)/fs;%设置时间subplot211;plot(t,x,'k');%画出带有趋势项的语音信号xline(t,xtrend,'color',[.6.6.6],'linewidth',3);%画出趋势项曲线ylim([-1.51]);title('带趋势项的语音信号');legend('带趋势项的语音信号','趋势项信号',4)xlabel('时间/s');ylabel('幅值');subplot212;plot(t,y,'k');%画出消除趋势项的语音信号yxlabel('时间/s');ylabel('幅值');title('消除趋势项的语音信号');第七十三页，共一百五十五页，编辑于2023年，星期二第七十四页，共一百五十五页，编辑于2023年，星期二例:读入数据文件bluesky3.wav（内容为男声“蓝天，白云”），加入线性趋势项x=0.00001+0.00002*n,并消除该趋势项clearall;clc;closeall;filedir=[];%指定文件路径filename='bluesky3.wav';%指定文件名fle=[filedirfilename];[s,fs]=wavread(fle);%读入数据文件t=(0:length(s)-1)/fs;%设置时间N=length(s);n=0:N-1;x=0.00001+0.00002*n;y=x'+s;第七十五页，共一百五十五页，编辑于2023年，星期二subplot(411)plot(s);subplot(412)plot(y);y=detrend(y);%消除线性趋势项y=y/max(abs(y));%幅值归一化subplot413;plot(t,x,'k');%画出带有趋势项的语音信号xtitle('趋势项');xlabel('时间/s');ylabel('幅值');subplot414;plot(t,y,'k');%画出消除趋势项的语音信号yxlabel('时间/s');ylabel('幅值');title('消除趋势项的语音信号');第七十六页，共一百五十五页，编辑于2023年，星期二第七十七页，共一百五十五页，编辑于2023年，星期二读入数据文件redriver2.wav

，加入多项式趋势项x=0.00002+0.00004*n;并消除该趋势项课堂练习：第七十八页，共一百五十五页，编辑于2023年，星期二例:读入数据文件bluesky3.wav（内容为男声“蓝天，白云”），加入多项式趋势项x=0.000001+0.000005*n+0.000000004*n.^2;并消除该趋势项filedir=[];%指定文件路径filename='bluesky3.wav';%指定文件名fle=[filedirfilename];[s,fs]=wavread(fle);%读入数据文件N=length(s);n=0:N-1;x=0.000001+0.000005*n+0.000000004*n.^2;y=x'+s;subplot(411)plot(s);subplot(412)plot(y);第七十九页，共一百五十五页，编辑于2023年，星期二[y,xtrend]=polydetrend(y,fs,2);%调用polydetrend消除趋势项t=(0:length(s)-1)/fs;%设置时间subplot413;plot(t,x,'k');%画出趋势项line(t,xtrend,'color',[.6.6.6],'linewidth',3);%画出趋势项曲线ylim([-1.51]);title('带趋势项的语音信号');legend('带趋势项的语音信号','趋势项信号',4)xlabel('时间/s');ylabel('幅值');subplot414;plot(t,y,'k');%画出消除趋势项的语音信号yxlabel('时间/s');ylabel('幅值');title('消除趋势项的语音信号');第八十页，共一百五十五页，编辑于2023年，星期二第八十一页，共一百五十五页，编辑于2023年，星期二课堂练习2读入数据文件redriver2.wav,加入多项式趋势项x=0.000001+0.000002*n+0.000000003*n.^2;,并消除该趋势项第八十二页，共一百五十五页，编辑于2023年，星期二5.5语音信号的预处理二—数字滤波器第八十三页，共一百五十五页，编辑于2023年，星期二5.5.1IIR低通、高通、带通和带阻滤波器设计第八十四页，共一百五十五页，编辑于2023年，星期二1.设计IIR数字滤波器的函数[N,Wn]=buttord(Wp,Ws,Rp,Rs);%设计参数为通带截止频率Wp，阻带截止Ws，通带最大衰减Rp，阻带最小衰减Rs的巴特沃斯滤波器。Wp,Ws为对归一化。例如当Wp=0.1 ，在程序中Wp=0.1。第八十五页，共一百五十五页，编辑于2023年，星期二第八十六页，共一百五十五页，编辑于2023年，星期二2.IIR数字滤波器的阶次选择函数第八十七页，共一百五十五页，编辑于2023年，星期二第八十八页，共一百五十五页，编辑于2023年，星期二例：设计一个巴特沃斯低通滤波器，它的Wp为500Hz，Ws为750Hz，采样频率fs=8000Hz，Rp=3dB,Rs=50dB,并画出幅频特性。clearall;clc;closeall;fp=500;fs=750;%设置滤波器的通带和阻带频率Fs=8000;Fs2=Fs/2;%采样频率FsWp=fp/Fs2;Ws=fs/Fs2;%把通带和阻带频率归一化Rp=3;Rs=50;%通带波纹和阻带衰减[n,Wn]=buttord(Wp,Ws,Rp,Rs);%求取滤波器阶数[b,a]=butter(n,Wn);%设计巴特沃斯型低通滤波器系数[db,mag,pha,grd,w]=freqz_m(b,a);%求滤波器的频率响应曲线plot(w*Fs2/pi,db);gridon;第八十九页，共一百五十五页，编辑于2023年，星期二巴特沃斯低通滤波器幅频特性第九十页，共一百五十五页，编辑于2023年，星期二课堂练习设计一个巴特沃斯低通滤波器，它的Wp为800Hz，Ws为1000Hz，采样频率fs=8000Hz，Rp=3dB,Rs=40dB,并画出幅频特性。第九十一页，共一百五十五页，编辑于2023年，星期二例：设计一个巴特沃斯带通滤波器，它的Wp为[500Hz,1000Hz]，Ws=[250Hz,1250Hz]，为采样频率fs=8000Hz，Rp=3dB,Rs=50dB,并画出幅频特性clearall;clc;closeall;fp=[500,1000];fs=[250,1250];%设置滤波器的通带和阻带频率Fs=8000;Fs2=Fs/2;%采样频率Wp=fp/Fs2;Ws=fs/Fs2;%把通带和阻带频率归一化Rp=3;Rs=50;%通带波纹和阻带衰减[n,Wn]=buttord(Wp,Ws,Rp,Rs);%求取滤波器阶数[b,a]=butter(n,Wn);%设计巴特沃斯型低通滤波器系数[db,mag,pha,grd,w]=freqz_m(b,a);%求滤波器的频率响应曲线plot(w*Fs2/pi,db);gridon;第九十二页，共一百五十五页，编辑于2023年，星期二巴特沃斯带通滤波器幅频特性第九十三页，共一百五十五页，编辑于2023年，星期二课堂练习：设计一个巴特沃斯带通滤波器，它的Wp为[600Hz,900Hz]，Ws=[30Hz,1300Hz]，为采样频率fs=8000Hz，Rp=3dB,Rs=40dB,并画出幅频特性第九十四页，共一百五十五页，编辑于2023年，星期二例：设计一个巴特沃斯高通滤波器，它的Wp为1000Hz，Ws=500Hz，为采样频率fs=8000Hz，Rp=3dB,Rs=50dB,并画出幅频特性clearall;clc;closeall;fp=1000;fs=500;%设置滤波器的通带和阻带频率Fs=8000;Fs2=Fs/2;%采样频率Wp=fp/Fs2;Ws=fs/Fs2;%把通带和阻带频率归一化Rp=3;Rs=50;%通带波纹和阻带衰减[n,Wn]=buttord(Wp,Ws,Rp,Rs);%求取滤波器阶数[b,a]=butter(n,Wn,'high');%设计巴特沃斯型高通滤波器系数[db,mag,pha,grd,w]=freqz_m(b,a);%求滤波器的频率响应曲线plot(w*Fs2/pi,db);gridon;第九十五页，共一百五十五页，编辑于2023年，星期二巴特沃斯高通滤波器幅频特性第九十六页，共一百五十五页，编辑于2023年，星期二课堂练习：设计一个巴特沃斯高通滤波器，它的Wp为800Hz，Ws=400Hz，为采样频率fs=8000Hz，Rp=3dB,Rs=30dB,并画出幅频特性第九十七页，共一百五十五页，编辑于2023年，星期二例：设计一个巴特沃斯带阻滤波器，它的Ws为[500Hz,1000Hz]，Wp=[250Hz,1250Hz]，为采样频率fs=8000Hz，Rp=3dB,Rs=50dB,并画出幅频特性clearall;clc;closeall;fp=[250,1250];fs=[500,1000];%设置滤波器的通带和阻带频率Fs=8000;Fs2=Fs/2;%采样频率Wp=fp/Fs2;Ws=fs/Fs2;%把通带和阻带频率归一化Rp=3;Rs=50;%通带波纹和阻带衰减[n,Wn]=buttord(Wp,Ws,Rp,Rs);%求取滤波器阶数[b,a]=butter(n,Wn,'stop');%设计巴特沃斯型带阻滤波器系数[db,mag,pha,grd,w]=freqz_m(b,a);%求滤波器的频率响应曲线plot(w*Fs2/pi,db);gridon;第九十八页，共一百五十五页，编辑于2023年，星期二巴特沃斯带阻滤波器幅频特性第九十九页，共一百五十五页，编辑于2023年，星期二课堂练习：设计一个巴特沃斯带阻滤波器，它的Ws为[600Hz,9000Hz]，Wp=[400Hz,1300Hz]，为采样频率fs=8000Hz，Rp=3dB,Rs=20dB,并画出幅频特性第一百页，共一百五十五页，编辑于2023年，星期二例：设计一个切比雪夫1型低通滤波器，它的Wp为500Hz，Ws为750Hz，采样频率fs=8000Hz，Rp=3dB,Rs=50dB,并画出幅频特性。clearall;clc;closeall;fp=500;fs=750;%设置滤波器的通带和阻带频率Fs=8000;Fs2=Fs/2;%采样频率Wp=fp/Fs2;Ws=fs/Fs2;%把通带和阻带频率归一化Rp=3;Rs=50;%通带波纹和阻带衰减[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs);%求取滤波器阶数[b,a]=cheby1(n,Rp,Wn);%设计契比雪夫II型低通滤波器系数[db,mag,pha,grd,w]=freqz_m(b,a);%求滤波器的频率响应曲线plot(w*Fs2/pi,db);gridon;第一百零一页，共一百五十五页，编辑于2023年，星期二切比雪夫1型低通滤波器幅频特性第一百零二页，共一百五十五页，编辑于2023年，星期二课堂练习：设计一个切比雪夫1型低通滤波器，它的Wp为400Hz，Ws为800Hz，采样频率fs=8000Hz，Rp=3dB,Rs=45dB,并画出幅频特性。第一百零三页，共一百五十五页，编辑于2023年，星期二例：设计一个切比雪夫2型低通滤波器，它的Wp为500Hz，Ws为750Hz，采样频率fs=8000Hz，Rp=3dB,Rs=50dB,并画出幅频特性。clearall;clc;closeall;fp=500;fs=750;%设置滤波器的通带和阻带频率Fs=8000;Fs2=Fs/2;%采样频率Wp=fp/Fs2;Ws=fs/Fs2;%把通带和阻带频率归一化Rp=3;Rs=50;%通带波纹和阻带衰减[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs);%求取滤波器阶数[b,a]=cheby2(n,Rs,Wn);%设计契比雪夫II型低通滤波器系数[db,mag,pha,grd,w]=freqz_m(b,a);%求滤波器的频率响应曲线plot(w*Fs2/pi,db);gridon;第一百零四页，共一百五十五页，编辑于2023年，星期二切比雪夫2型低通滤波器幅频特性第一百零五页，共一百五十五页，编辑于2023年，星期二课堂练习：设计一个切比雪夫2型低通滤波器，它的Wp为400Hz，Ws为800Hz，采样频率fs=8000Hz，Rp=3dB,Rs=45dB,并画出幅频特性。第一百零六页，共一百五十五页，编辑于2023年，星期二例：设计一个椭圆型低通滤波器，它的Wp为500Hz，Ws为750Hz，采样频率fs=8000Hz，Rp=3dB,Rs=50dB,并画出幅频特性。clearall;clc;closeall;fp=500;fs=750;%设置滤波器的通带和阻带频率Fs=8000;Fs2=Fs/2;%采样频率Wp=fp/Fs2;Ws=fs/Fs2;%把通带和阻带频率归一化Rp=3;Rs=50;%通带波纹和阻带衰减[n,Wn]=ellipord(Wp,Ws,Rp,Rs);%求取滤波器阶数[b,a]=ellip(n,Rp,Rs,Wn);%设计巴特沃斯型低通滤波器系数[db,mag,pha,grd,w]=freqz_m(b,a);%求滤波器的频率响应曲线plot(w*Fs2/pi,db);gridon;第一百零七页，共一百五十五页，编辑于2023年，星期二椭圆型低通滤波器幅频特性第一百零八页，共一百五十五页，编辑于2023年，星期二课堂练习：设计一个椭圆型低通滤波器，它的Wp为400Hz，Ws为800Hz，采样频率fs=8000Hz，Rp=3dB,Rs=45dB,并画出幅频特性。第一百零九页，共一百五十五页，编辑于2023年，星期二%pr5_5_1clearall;clc;closeall;fp=500;fs=750;%设置滤波器的通带和阻带频率Fs=8000;Fs2=Fs/2;%采样频率Wp=fp/Fs2;Ws=fs/Fs2;%把通带和阻带频率归一化Rp=3;Rs=50;%通带波纹和阻带衰减[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs);%求取滤波器阶数[b,a]=cheby2(n,Rs,Wn);%设计契比雪夫II型低通滤波器系数[db,mag,pha,grd,w]=freqz_m(b,a);%求滤波器的频率响应曲线第一百一十页，共一百五十五页，编辑于2023年，星期二filedir=[];%指定文件路径filename='bluesky3.wav';%指定文件名fle=[filedirfilename]%构成路径和文件名的字符串[s,fs]=wavread(fle);%读入数据文件s=s-mean(s);%消除直流分量s=s/max(abs(s));%幅值归一化N=length(s);%求出信号长度t=(0:N-1)/fs;%设置时间第一百一十一页，共一百五十五页，编辑于2023年，星期二y=filter(b,a,s);%把语音信号通过滤波器wlen=200;inc=80;nfft=512;%设置帧长,帧移和nfft长win=hann(wlen);%设置窗函数d=stftms(s,win,nfft,inc);%原始信号的STFT变换fn=size(d,2);%获取帧数frameTime=(((1:fn)-1)*inc+nfft/2)/Fs;%计算每帧对应的时间--时间轴刻度W2=1+nfft/2;%计算频率轴刻度n2=1:W2;freq=(n2-1)*Fs/nfft;d1=stftms(y,win,nfft,inc);%滤波后信号的STFT变换第一百一十二页，共一百五十五页，编辑于2023年，星期二%作图figure(1)plot(w/pi*Fs2,db,'k','linewidth',2)grid;axis([04000-1005]);title('低通滤波器的幅值响应曲线')xlabel('频率/Hz');ylabel('幅值/dB');figure(2)subplot211;plot(t,s,'k');title('纯语音信号：男声“蓝天，白云”')xlabel(['时间/s'10'(a)']);ylabel('幅值')subplot212;imagesc(frameTime,freq,abs(d(n2,:)));axisxytitle('纯语音信号的语谱图')xlabel(['时间/s'10'(b)']);ylabel('频率/Hz')第一百一十三页，共一百五十五页，编辑于2023年，星期二m=256;LightYellow=[0.60.60.6];MidRed=[000];Black=[0.50.71];Colors=[LightYellow;MidRed;Black];colormap(SpecColorMap(m,Colors));figure(3)subplot211;plot(t,y,'k');title('滤波后的语音信号')xlabel(['时间/s'10'(a)']);ylabel('幅值')subplot212;imagesc(frameTime,freq,abs(d1(n2,:)));axisxytitle('滤波后语音信号的语谱图')xlabel(['时间/s'10'(b)']);ylabel('频率/Hz')m=256;LightYellow=[0.60.60.6];MidRed=[000];Black=[0.50.71];Colors=[LightYellow;MidRed;Black];colormap(SpecColorMap(m,Colors));ylim([01000]);第一百一十四页，共一百五十五页，编辑于2023年，星期二第一百一十五页，共一百五十五页，编辑于2023年，星期二第一百一十六页，共一百五十五页，编辑于2023年，星期二第一百一十七页，共一百五十五页，编辑于2023年，星期二课堂练习：1、设计一个巴特沃斯低通滤波器，它的Wp为750Hz,Ws为1000Hz，采样频率fs为8000Hz，Rp和Rs分别为3dB和50dB，并对bluesky3。Wav数据进行滤波。2、设计一个椭圆型低通滤波器，它的Wp为750Hz,Ws为1000Hz，采样频率fs为8000Hz，Rp和Rs分别为3dB和50dB，并对bluesky3。Wav数据进行滤波。第一百一十八页，共一百五十五页，编辑于2023年，星期二5.5.2FIR低通、高通、带通和带阻滤波器的设计第一百一十九页，共一百五十五页，编辑于2023年，星期二第一百二十页，共一百五十五页，编辑于2023年，星期二2.用窗函数设计FIR滤波器步骤第一百二十一页，共一百五十五页，编辑于2023年，星期二第一百二十二页，共一百五十五页，编辑于2023年，星期二第一百二十三页，共一百五十五页，编辑于2023年，星期二3.用窗函数设计FIR滤波器的方法第一百二十四页，共一百五十五页，编辑于2023年，星期二第一百二十五页，共一百五十五页，编辑于2023年，星期二第一百二十六页，共一百五十五页，编辑于2023年，星期二第一百二十七页，共一百五十五页，编辑于2023年，星期二第一百二十八页，共一百五十五页，编辑于2023年，星期二例：用窗函数法设计数字FIR低通滤波器,窗长N取偶数：它的Wp为500Hz，Ws为1000Hz，采样频率fs=8000Hz，Rp=3dB,Rs=40dB,并画出幅频特性。通过查表，选汉宁窗：As=40;Fs=8000;Fs2=Fs/2;%阻带最小衰减和采样频率fp=500;fs=1000;%通带阻带频率wp=fp/Fs2*pi;ws=fs/Fs2*pi;%转为圆频率deltaw=ws-wp;%求取过渡带wc=(wp+ws)/2;%求取截止频率N0=ceil(6.2*pi/deltaw); %按汉宁窗计算所需滤波器长度N=N0+mod(N0,2);%采用第