课程设计（论文）-数字信号的处理与滤波设计.doc

广西工学院课程设计用纸 装 订 线 数字信号的处理与滤波 摘要： 数字信号处理是利用计算机或专用处理设备，以数值计算的方法，对信号进行采集、滤波、增强、压缩、估值和识别等加工处理，借以达到提取信息和便于应用的目的，其应用范围涉及几乎所有的工程技术领域。语音信号的处理与滤波是数字信号处理中一个重要的组成部分。本课题主要研究的就是语音信号的处理与滤波，然后设计一个用户界面实现语音的处理与滤波的功能。使用matlab进行滤波器的设计，使用matlab的guide工具设计用户界面。关键词：数字滤波器；双线性变换；窗函数；凯泽窗；guide一 设计题目所选题目为设计五：语音信号的处理与滤波 内容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。基本要求：1语音信号的采集利用windows下的录音机，录制一段自己的话音，时间在s内。然后在matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，要求理解掌握采样频率、采样位数等概念。2语音信号的频谱分析要求首先画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性，从而加深对频谱特性的理解。其程序如下：3设计数字滤波器和画出其频率响应给出各滤波器的性能指标：(1)低通滤波器性能指标fb1 000 hz，fc1 200 hz，as100 db，ap1 db。(2)高通滤波器性能指标fc4 800 hz，fb5 000 hz as100 db，ap1 db。(3)带通滤波器性能指标fb11 200 hz，fb23 000 hz，fc11 000 hz，fc23 200 hz，as100 db，ap1 db。要求用窗函数法和双线性变换法设计上面要求的3种滤波器。在matlab中，可以利用函数fir1设计fir滤波器，可以利用函数butte，cheby1和ellip设计iir滤波器；利用matlab中的函数freqz画出各滤波器的频率响应。4用滤波器对信号进行滤波要求用自己设计的各种滤波器分别对采集的信号进行滤波，在matlab中，fir滤波器利用函数fftfilt对信号进行滤波，iir滤波器利用函数filter对信号进行滤波。 5比较滤波前后语音信号的波形及频谱要求在一个窗口同时画出滤波前后的波形及频谱。 6回放语音信号在matlab中，函数sound可以对声音进行回放。其调用格式：sound(x，fs，bits)；可以感觉滤波前后的声音有变化。 7设计系统界面（用matlab中guide工具，要求自学）为了使编制的程序操作方便，要求有能力的学生，设计处理系统的用户界面。在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数，显示滤波器的频率响应，选择信号等。二 设计所涉及原理1 对信号进行采样抽样原理：要保证从信号抽样后的离散时间信号无失真地恢复原始时间连续信号（即抽样不会导致任何信息丢失），必须满足：信号是频带受限的(信号频率区间有限)；采样率 至少是信号最高频率的两倍。2 iir滤波器 语音信号的处理与滤波一 设计题目所选题目为设计五：语音信号的处理与滤波 内容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。基本要求：1语音信号的采集利用windows下的录音机，录制一段自己的话音，时间在s内。然后在matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，要求理解掌握采样频率、采样位数等概念。2语音信号的频谱分析要求首先画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性，从而加深对频谱特性的理解。其程序如下：3设计数字滤波器和画出其频率响应给出各滤波器的性能指标：(1)低通滤波器性能指标fb1 000 hz，fc1 200 hz，as100 db，ap1 db。(2)高通滤波器性能指标fc4 800 hz，fb5 000 hz as100 db，ap1 db。(3)带通滤波器性能指标fb11 200 hz，fb23 000 hz，fc11 000 hz，fc23 200 hz，as100 db，ap1 db。要求用窗函数法和双线性变换法设计上面要求的3种滤波器。在matlab中，可以利用函数fir1设计fir滤波器，可以利用函数butte，cheby1和ellip设计iir滤波器；利用matlab中的函数freqz画出各滤波器的频率响应。4用滤波器对信号进行滤波要求用自己设计的各种滤波器分别对采集的信号进行滤波，在matlab中，fir滤波器利用函数fftfilt对信号进行滤波，iir滤波器利用函数filter对信号进行滤波。 5比较滤波前后语音信号的波形及频谱要求在一个窗口同时画出滤波前后的波形及频谱。 6回放语音信号在matlab中，函数sound可以对声音进行回放。其调用格式：sound(x，fs，bits)；可以感觉滤波前后的声音有变化。 7设计系统界面（用matlab中guide工具，要求自学）为了使编制的程序操作方便，要求有能力的学生，设计处理系统的用户界面。在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数，显示滤波器的频率响应，选择信号等。二 设计所涉及原理1 对信号进行采样抽样原理：要保证从信号抽样后的离散时间信号无失真地恢复原始时间连续信号（即抽样不会导致任何信息丢失），必须满足：信号是频带受限的(信号频率区间有限)；采样率 至少是信号最高频率的两倍。2 iir滤波器(1) 脉冲响应不变法：用数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t),让h(n)正好等于ha(t)的采样值，即 h(n)=ha(nt) 其中t为采样间隔，如果以ha(s)及h(z)分别表示ha(t)的拉式变换及h(n)的z变换，则 (2)双线性变换法：s平面与z平面之间满足以下映射关系：s平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在混叠问题。 双线性变换时一种非线性变换 ，这种非线性引起的幅频特性畸变可通过预畸而得到校正。 iir低通、高通、带通数字滤波器设计采用双线性原型变换公式：变换类型 变换关系式 备 注 低通 高通 带通 带通的上下边带临界频率 以低通数字滤波器为例，将设计步骤归纳如下：1.确定数字滤波器的性能指标：通带临界频率fp、阻带临界频率fr；通带内的最大衰减ap；阻带内的最小衰减ar；采样周期t； （3）.确定相应的数字角频率，p=2fpt；r=2frt；（4）.计算经过预畸的相应模拟低通原型的频率， ； （5）.根据p和r计算模拟低通原型滤波器的阶数n，并求得低通原型的传递函数ha(s)； （6）.用上面的双线性变换公式代入ha(s)，求出所设计的传递函数h(z)； （7）.分析滤波器特性，检查其指标是否满足要求。3 fir滤波器 线性相位实系数fir滤波器按其n值奇偶和h(n)的奇偶对称性分为四种： （1）、h(n)为偶对称，n为奇数h(ej)的幅值关于=0，2成偶对称。 （2）、h(n)为偶对称，n为偶数h(ej)的幅值关于=成奇对称，不适合作高通。 （3）、h(n)为奇对称，n为奇数h(ej)的幅值关于=0，2成奇对称，不适合作高通和低通。 （4）、h(n)为奇对称，n为偶数h(ej) =0、20，不适合作低通。(5) 窗口法 窗函数法设计线性相位fir滤波器步骤 确定数字滤波器的性能要求：临界频率k，滤波器单位脉冲响应长度n； 根据性能要求，合理选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应hd(ej)的幅频特性和相频特性； 求理想单位脉冲响应hd(n)，在实际计算中，可对hd(ej)按m(m远大于n)点等距离采样，并对其求idft得hm(n)，用hm(n)代替hd(n)； 选择适当的窗函数w(n)，根据h(n)= hd(n)w(n)求所需设计的fir滤波器单位脉冲响应； 求h(ej)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度n，重复上述设计过程，以得到满意的结果。 窗函数的傅式变换w(ej)的主瓣决定了h(ej)过渡带宽。w(ej)的旁瓣大小和多少决定了h(ej)在通带和阻带范围内波动幅度，常用的几种窗函数有： 矩形窗 w(n)=rn(n)； hanning窗 ； hamming窗 ； blackmen窗 ； kaiser窗 。 式中io(x)为零阶贝塞尔函数。三 设计工具matlab 7.0.1国际标准版四 设计进程说明设计主要分为两大块：（一）主程序 这部分主要实现了设计要求的第16个要求（见附录一） 1 程序中首先用函数wavread对所录制的语音信号（大小为5秒左右）进行采样，采样频率为fs=22050.2 画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性。3 利用函数fir1设计fir滤波器，并且用设计的滤波器对语音信号进行滤波（注：主程序中所使用滤波器为fir滤波器hamming窗）4 比较滤波前后语音信号的波形及频谱 5 用sound(x，fs，bits)函数对所滤声音进行回放，并与原始声音进行比较（二）用户界面（包含所界面程序 见附录二） 这一部分主要实现了设计要求的第7步，界面中使用了所给定的三种指标分别对录制语音进行滤波。界面中可以选择滤波器的类型，可以显示滤波前的声音频域图形，滤波后的频域波形。还可以对原始波形和滤波后的波形进行比较等功能。 附录一 主程序 fs=22050;x1=wavread(f:io4.wav);t=0:1/22050:(size(x1)-1)/22050;au=0.03;d=au*cos(2*pi*5000*t);x2=x1+d;%加入噪声%用滤波器进行滤波（汉宁窗）wp=0.25*pi;ws=0.3*pi;wdelta=ws-wp;n=ceil(6.6*pi/wdelta); %取整wn=(0.2+0.3)*pi/2; %数字截止频率b=fir1(n,wn/pi,hamming(n+1); %选择窗函数，并归一化截止频率figure(1)freqz(b,1,512);f2=fftfilt(b,x2);%显示波形图figure(2)subplot(2,1,1);plot(t,x2);title(滤波前的时域波形);subplot(2,1,2);plot(t,f2);title(滤波后的时域波形);%显示频谱图f0=fft(f2,1024);f=fs*(0:511)/1024;figure(3)y2=fft(x2,1024);subplot(2,1,1);plot(f,abs(y2(1:512);title(滤波前的频谱);xlabel(hz);ylabel(fuzhi); %画出滤波前的频谱图subplot(2,1,2);plot(f,abs(f0(1:512);title(滤波后的频谱);xlabel(hz);ylabel(fuzhi); %画出滤波后的频谱图%声音比较chooise=input(1、原始语音2、加噪后语音3、滤波后语音);if(chooise=1)sound(x1,22050); %播放原始语音endif(chooise=2) sound(x2,22050); %播放加噪后的语音endif(chooise=3) sound(f2,22050); %播放滤波后的信号end 附录二 界面与程序 程序摘要如下： function pushbutton1_callback(hobject, eventdata, handles)%滤波后的声音if handles.but1=1sound(handles.sound1, 22050)endif handles.but2=1sound(handles.sound2, 22050)endif handles.but3=1sound(handles.sound3*5, 22050)end if handles.but4=1sound(handles.sound4, 22050)end if handles.but5=1sound(handles.sound5, 22050)end if handles.but6=1sound(handles.sound6, 22050)end %滤波前的语音function pushbutton2_callback(hobject, eventdata, handles)fs=8000;y, fs=wavread(io4.wav);sound(y, fs); %滤波前的信息function pushbutton3_callback(hobject, eventdata, handles)fs=8000;y, fs=wavread(io4.wav);time=(1:length(y)/fs; _set(gcf, currentaxes,handles.axes1);plot(time, y); title(时域);y=fft(y,8000);amplitude=abs(y);set(gcf, currentaxes,handles.axes2);plot(amplitude);title(频域);%滤波后的信息 function pushbutton4_callback(hobject, eventdata, handles) if handles.but1=1fs=22050;wp=1000;ws=1200;rp=1;rs=100;wc=2*fs*tan(2*pi*wp/(2*fs) %先转化角频率再转换为模拟域频率wt=2*fs*tan(2*pi*ws/(2*fs) n,wn=cheb1ord(wc,wt,1,100,s); %在模拟域计算阶数和固有频率b,a=cheby1(n,1,wn,low,s); %利用上面得到的结果设计滤波器bz,az=bilinear(b,a,8000); %用双线性变换法得到数字滤波器的原型y, fs=wavread(io4.wav);x=filter(bz,az,y);x=fft(x,4096);set(gcf, currentaxes,handles.axes3)plot(x);title(滤波后语音信号波形);set(gcf, currentaxes,handles.axes4)plot(abs(x);title(滤波后语音信号的谱)handles.sound1=x;guidata(hobject, handles);end if handles.but2=1fs=22050;wp=1200 3000/4000;ws=1000 3200/4000;rp=1;rs=100;n,wn=buttord(wp,ws,rp,rs);b,a=butter(n,wn);y, fs=wavread(io4.wav);x=filter(b,a,y);x=fft(x,4096);set(gcf, currentaxes,handles.axes3)plot(x);title(滤波后语音信号波形);set(gcf, currentaxes,handles.axes4)plot(abs(x);title(滤波后语音信号的谱)handles.sound2=x;guidata(hobject, handles);endif handles.but3=1fs=22050;wp=5000;ws=4800;rp=1;rs=100;wc=2*fs*tan(2*pi*wp/(2*fs) %先转化角频率再转换为模拟域频率wt=2*fs*tan(2*pi*ws/(2*fs) n,wn=cheb1ord(wc,wt,1,100,s); %在模拟域计算阶数和固有频率b,a=cheby1(n,1,wn,high,s); %利用上面得到的结果设计滤波器bz,az=bilinear(b,a,22050); %用双线性变换法得到数字滤波器的原型 y, fs=wavread(io4.wav);x=filter(bz,az,y);x=fft(x);set(gcf, currentaxes,handles.axes3)plot(x);title(滤波后语音信号波形);set(gcf, currentaxes,handles.axes4)plot(abs(x);title(滤波后语音信号的谱)handles.sound3=x;guidata(hobject, handles);endif handles.but4=1 y=wavread(io4.wav); fs=1000;fp=1200;fs=22050;ws=fs/fs*2*pi; % wp=fp/fs*2*pi;w1=(ws+wp)/2; width=wp-ws; % n=ceil(12.8*pi/width)+1; window=kaiser(n,10.056);%w1=w1/pi;b=fir1(n-1,w1,window);y1=fftfilt(b,y);% y1=fft(y1,2048);set(gcf,currentaxes,handles.axes3);plot(y1); % title(滤波后的时域波形);set(gcf,currentaxes,handles.axes4);plot(abs(y1); title(滤波后的频域波形); handles.sound4=y1;guidata(hobject, handles);endif handles.but6=1fs=22050;y,fs=wavread(io4.wav); fs=5000;fp=4800;fs=22050;ws=fs/fs*2*pi; wp=fp/fs*2*pi;w1=(ws+wp)/2; %理想高通滤波器中心频率width=ws-wp;n=ceil(12.8*pi/width)+1; window=kaiser(n,10.056);%window=kaiser(n,handles.beta);w1=w1/pi;b=fir1(n-1,w1,high,window);y1=fftfilt(b,y);y1=fft(y1,2048);set(gcf,currentaxes,handles.axes3);plot(y1); % title(滤波后的时域波形);set(gcf,currentaxes,handles.axes4);plot(abs(y1); title(滤波后的频域波形); handles.sound6=y1;guidata(hobject, handles);endif handles.but5=1 %通过fir带通滤波器%输入带通滤波器的通带截止频率fp1，fp2阻带截止频率fs1,fs2和beta值fs=22050;fs=1000;fp=1200;fp2=3000;fs2=3200;y,fs=wavread(io4.wav); ws1=fs/fs*2*pi; wp1=fp/fs*2*pi;wp2=fp2/fs*2*pi; ws2=fs2/fs*2*pi;w1=(ws1+wp1)/2; w2=(wp2+ws2)/2; %理想带通滤波器上下截止频率width= min(wp1-ws1),(ws2-wp2);%过渡带宽n=ceil(12.8*pi/width)+1; window=kaiser(n,10.056);w1=w1/pi;w2=w2/pi;b=fir1(n-1,w1 w2,bandpass,window);y1=fftfilt(b,y);y1=fft(y1,2048);set(gcf,currentaxes,handles.axes3);plot(y1);%滤波后波形title(滤波后的时域波形);set(gcf,currentaxes,handles.axes4);plot(abs(y1);%滤波后频谱title(滤波后的频域波形);handles.sound5=y1;guidata(hobject, handles);end function radiobutton5_callback(hobject, eventdata, handles)%fir低通handles.but4=get(hobject,value)handles.but1=0;handles.but2=0;handles.but3=0;handles.but5=0;handles.but6=0;guidata(hobject, handles); y=wavread(io4.wav); fs=1000;fp=1200;fs=22050;ws=fs/fs*2*pi; % wp=fp/fs*2*pi;w1=(ws+wp)/2; width=wp-ws; % n=ceil(12.8*pi/width)+1; window=kaiser(n,10.056);%w1=w1/pi;b=fir1(n-1,w1,window);y1=fftfilt(b,y);% y1=fft(y1,2048);set(gcf,currentaxes,handles.axes1);plot(y1);axis(0,8000,-300,10); grid; xlabel(频率/hz) ylabel(幅度/db) % - executes on button press in radiobutton5. function radiobutton2_callback(hobject, eventdata, handles)handles.but1=get(hobject,value)handles.but2=0;handles.but3=0;handles.but4=0;handles.but5=0;handles.but6=0;guidata(hobject, handles); fs=8000;wp=1000;ws=1200;rp=1;rs=100;wc=2*fs*tan(2*pi*wp/(2*fs) %先转化角频率再转换为模拟域频率wt=2*fs*tan(2*pi*ws/(2*fs) n,wn=cheb1ord(wc,wt,1,100,s); %在模拟域计算阶数和固有频率b,a=cheby1(n,1,wn,low,s); %利用上面得到的结果设计滤波器num,den=bilinear(b,a,8000); %用双线性变换法得到数字滤波器的原型h,w=freqz(num,den); %求幅频响应set(gcf, currentaxes,handles.axes1);f=w/pi*4000; %将角频率转化为hz，500是采样频率的一半plot(f,20*log10(abs(h); %将幅频响应表示成分贝并画出图形axis(0,5000,-300,10); grid; xlabel(频率/hz) ylabel(幅度/db)% hobject handle to radiobutton2 (see gcbo)% eventdata reserved - to be defined in a future version of matlab% handles structure with handles and user data (see guidata) % hint: get(hobject,value) returns toggle state of radiobutton2 % - executes on button press in radiobutton3.function radiobutton3_callback(hobject, eventdata, handles)handles.but2=get(hobject,value)handles.but1=0;handles.but3=0;handles.but4=0;handles.but5=0;handles.but6=0; guidata(hobject, handles); wp=1200 3000/4000;ws=1000 3200/4000;rp=1;rs=100;n,wn=buttord(wp,ws,rp,rs);b,a=butter(n,wn);h w=freqz(b,a); set(gcf, currentaxes,handles.axes1);f=w/pi*4000; plot(f,20*log10(abs(h); %将幅频响应表示成分贝并画出图形grid;axis(0,5000,-500,10); function radiobutton4_callback(hobject, eventdata, handles)handles.but3=get(hobject,value)handles.but1=0;handles.but2=0;handles.but4=0;handles.but5=0;handles.but6=0;guidata(hobject, handles);fs=22050;wp=5000;ws=4800;rp=1;rs=100;wc=2*fs*tan(2*pi*wp/(2*fs) %先转化角频率再转换为模拟域频率wt=2*fs*tan(2*pi*ws/(2*fs) n,wn=cheb1ord(wc,wt,1,100,s); %在模拟域计算阶数和固有频率b,a=cheby1(n,1,wn,high,s); %利用上面得到的结果设计滤波器num,den=bilinear(b,a,22050); %用双线性变换法得到数字滤波器的原型h,w=freqz(num,den); %求幅频响应set(gcf, currentaxes,handles.axes1);f=w/pi*11025; %将角频率转化为hz，500是采样频率的一半plot(f,20*log10(abs(h); %将幅频响应表示成分贝并画出图形axis(0,8000,-300,10); grid; xlabel(频率/hz) ylabel(幅度/db) function radiobutton6_callback(hobject, eventdata, handles)%通过fir带通滤波器%输入带通滤波器的通带截止频率fp1，fp2阻带截止频率fs1,fs2和beta值handles.but5=get(hobject,value)handles.but1=0;handles.but2=0;handles.but4=0;handles.but6=0;handles.but3=0;guidata(hobject, handles);fs=22050;fs=1000;fp=1200;fp2=3000;fs2=3200;y,fs=wavread(io4.wav); ws1=fs/fs*2*pi; wp1=fp/fs*2*pi;wp2=fp2/fs*2*pi; ws2=fs2/fs*2*pi;w1=(ws1+wp1)/2; w2=(wp2+ws2)/2; %理想带通滤波器上下截止频率width= min(wp1-ws1),(ws2-wp2);%过渡带宽n=ceil(12.8*pi/width)+1; window=kaiser(n,10.