数字信号处理与滤波设计课程设计.docx

数字信号处理与滤波设计课程设计 装 订 线 数字信号的处理与滤波 摘要： 数字信号处理是利用计算机或专用处理设备，以数值计算的方法，对信号进行采集、滤波、增加、压缩、估值和识别等加工处理，借以达到提取信息和便于应用的目的，其应用范围涉及几乎全部的工程技术领域。语音信号的处理与滤波是数字信号处理中一个重要的组成部分。本课题主要探讨的就是语音信号的处理与滤波，然后设计一个用户界面实现语音的处理与滤波的功能。运用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); %播放原始语音 end if(chooise=2) sound(x2,22050); %播放加噪后的语音 end if(chooise=3) sound(f2,22050); %播放滤波后的信号 end 附录二 界面与程序 程序摘要如下： function pushbutton1_Callback(hObject, eventdata, handles)%滤波后的声音 if handles.but1=1 sound(handles.sound1, 22050) end if handles.but2=1 sound(handles.sound2, 22050) end if handles.but3=1 sound(handles.sound3*5, 22050) end if handles.but4=1 sound(handles.sound4, 22050) end if handles.but5=1 sound(handles.sound5, 22050) end if handles.but6=1 sound(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=1 fs=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=1 fs=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); end if handles.but3=1 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); %利用上面得到的结果设计滤波器 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); end if 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); end if handles.but6=1 Fs=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); end if 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.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.axes1); plot(y1); axis(0,8000,-300,10); grid; xlabel(频率/Hz) ylabel(幅度/dB) % 高通FIR滤波器 function radiobutton7_Callback(hObject, eventdata, handles) handles.but6=get(hObject,value) handles.but1=0; handles.but2=0; handles.but4=0; handles.but3=0; handles.but5=0; guidata(hObject, handles); Fs=20000; 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.axes1); plot(y1); axis(0,8000,-300,10); grid; xlabel(频率/Hz) ylabel(幅度/dB) 参考文献： 赵红怡 张常年编著 数字信号处理及其MATLAB实现 北京：化学工业出版社 2002年1月 程佩青编著 数字信号处理教程（其次版） 北京:清华高校出版社 2005年6月 钱同惠主编 数字信号处理 北京： 机械工业出版社 2005年2月 王宏编著 MATLAB6.5及其在信号处理中的应用 北京：清华高校出版社 2005年3月 闻新 周露 张鸿编著 MATLAB科学图形构建基础与应用 北京：科学出版社 2002年6月 苏晓生编著 北京：科学出版社 2002年 陈亚勇等编著 MATLAB信号处理详解 北京：人民邮电出版社 2001年9月 楼顺天 李博菡编著 基于MATLAB的系统分析与设计信号处理 西安：西安电子科技高校出版社 2001年9月 张威编著 MATLAB基础与编程入门 西安：西安电子科技高校出版社 2005年1月 闻新等编著 MATLAB模糊逻辑工具箱的分析与应用 北京：科学出版社 2001 孙洪等译 数字信号处理试验指导书（MATLAB版）北京：电子工业出版社 2005年1月 卢光跃等编著 数字信号处理 北京：北京邮电高校出版社 2003 广xi工学院