林数字信号处理课程设计

广西工学院

数字信号处理 课程设计

说 明 书

设计题目 语音信号的处理与滤波系 别 计算机工程系 专业班级 通信062班 学生姓名 林芳慧 学 号 200600402059 指导教师 周坚和 日 期 2009年1月2日

2

目录

(一) (二) (三)

(四) (五) (六) (七)

设计目的…………………………………………………………………3

设计内容…………………………………………………………………3

设计原理…………………………………………………………………4

设计步骤与程序………………………………………………………5

设计思考………………………………………………………………19

结论………………………………………………………………………19

参考文献………………………………………………………………20

3

数字信号处理课程设计

一、设计目的

这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。 二．设计内容

语音信号的处理与滤波:

1 熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数。在MATLAB环境中，使用声音相关函数录制2秒左右自己的声音，抽样率是8000Hz/s。并考虑如何解决一个实际问题：录制刚开始时，常会出现实际发出声音落后录制动作半拍的现象，如何排除对这些无效点的采样？ 2．语音信号的频谱分析

分别取8000个和16000个数据进行频谱分析，得到幅度和相位谱，比较二者异同并分析原因。要求首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速傅立叶变换，得到信号的频谱特性；从而加深对频谱特性的理解。 3．设计数字滤波器和画出频率响应

针对电话信道（最高3500Hz），设计一个FIR或IIR滤波器进行滤波，把抽样率转变为7000Hz/s，并进行频谱分析，得到幅度和相位谱。根据语音信号的特点给出有关滤波器的性能指标：（1）低通滤波器性能指标，fp=1000Hz，fc=1200 Hz， As=80dB，Ap=1dB；（2）高通滤波器性能指标，fc=2800 Hz，fp=3000 Hz As=80dB，Ap=1dB；（3）带通滤波器性能指标，fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz，As=80dB，Ap=1dB。要求学生首先用窗函数法设计上面要求的三种滤波器，在MATLAB中，可以利用函数fir1设计FIR滤波器，然后在用双线性变换法设计上面要求的三种滤波器；之后再利用函数butter和cheby1设计上面要求的三种IIR滤波器。最后，利用MATLAB中的函数freqz画出各滤

4

波器的频率响应。

4．用滤波器对信号进行滤波

比较FIR和IIR两种滤波器的性能，然后用性能好的各滤波器分别对采集的信号进行滤波，在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。 5．比较滤波前后语音信号的波形及频谱

要求在一个窗口同时画出滤波前后的波形及频谱。 6．回放语音信号

在MATLAB中，函数sound可以对声音进行回放。其调用格式： sound(x,fs,bits)；可以感觉滤波前后的声音有变化。

三．设计原理

1． IIR滤波器双线性变换法设计原理：双线性变换法中S平面与z平面之间满足以下映射关系：

s平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在混叠问题。双线性变换时一种非线性变换

，这种非线性引起的幅频特性畸变可通过预畸而得到校正。

双线性法设计IIR数字滤波器的步骤：(1）将数字滤波器的频率指标{ k}由Wk＝（2/T）*tan(wk),转换为模拟滤波器的频率指标{k}. (2) 根据通带临界频率和阻带临界频率计算模拟低通原型滤波器的阶数N，并求得低通原型的传递函数H(s)； (3) 用上面的双线性变换公式代入H(s)，求出所设计的传递函数H(z).(4) 分析滤波器特性，检查其指标是否满足要求。 2. 窗函数法设计线性相位FIR滤波器步骤

(1)确定数字滤波器的性能要求：临界频率{ωk}，滤波器单位脉冲响应长度N； 根据性能要求，合理选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应Hd(ejω)的幅频特性和相频特性；

(2)求理想单位脉冲响应hd(n)，在实际计算中，可对Hd(ejω)按M(M远大于N)点等距离采样，并对其求IDFT得hM(n)，用hM(n)代替hd(n)；

5

(3)选择适当的窗函数w(n)，根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为Δw，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。根据h(n)= hd(n)w(n)求所需设计的FIR滤波器单位脉冲响应； 求H(ejω)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。 窗函数的傅式变换W(ejω)的主瓣决定了H(ejω)过渡带宽。W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度.。

四．课程设计的步骤 1.语音信号的采集：

（1）在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。

采样频率即取样频率,指每秒钟取得声音样本的次数。它的采样频率越高,声音的质量也就越好,但人耳的分辨率很有限,所以太高的频率就分辨不出好坏来。对于高于48KHz的采样频率人耳已无法辨别出来了。由采样定理知，采样频率要大于信号最高频率的两倍，否则会产生混叠，人的声音频率在300Hz-3400Hz所以取采样频率为8000Hz。这次实验是利用MATLAB来录制的.采样位数可以理解为声卡处理声音的解析度。这个数值越大，解析度就越高，录制和回放的声音就越真实。电脑中的声音文件是用数字0和1来表示的。所以在电脑上录音的本质就是把模拟声音信号转换成数字信号。反之，在播放时则是把数字信号还原成模拟声音信号输出。实验中wavread函数返回的Nbits为16位，声音效果还是很好的。

%我的录音程序

fs=8000;%采样频率 fprintf('开始录音');

y=wavrecord(2*fs,fs,'double')%录制声音2秒 fprintf('录音结束');

wavplay(y,fs);%播放录音

wavwrite(y,fs,'d:/shengyin');%声音的存储 wavread('d:/shengyin');%读取声音 figure(1);

plot(abs(y)); %画出声音采样后的时域波形 grid on;%添加网格

title('语音信号的时域波形');%给图形加注标签 z=fft(y,8192);%对信号做8192点FFT变换 f=fs*(0:4095)/8192;%采样点数 figure(2);

plot(f,abs(z(1:4096)));%绘制原始语音信号的幅频响应图

6

axis([0 4000 0 20]);

title('FFT后声音的频谱特性');%给图形加注标签 grid on;%添加网格

%录制刚开始时，常会出现实际发出声音落后录制动作半拍的现象，现在排除对这些无效点 fs=8000;

size=wavread('d:/shengyin','size')%声音的数据长度 z=wavread('d:/shengyin',[2000 16000]);%读取声音 wavplay(z,fs);%播放声音

wavwrite(z,fs,'d:/shengyin1');%声音的存储 wavread('d:/shengyin1');%读取声音 figure(3);

plot(abs(z)); %画出声音采样后的时域波形 grid on;%添加网格

title('排除无效点后语音信号的时域波形');%给图形加注标签说明 z=fft(y,8192); %对信号做fs点FFT变换 figure(4);

f=fs*(0:4095)/8192;%采样点数

plot(f,abs(z(1:4096))); %绘制语音信号的幅频响应图 axis([0 4000 0 20]);

title('排除无效点经FFT后声音的频谱特性');%给图形加注标签说明 grid on;%添加网格

语音信号的时域波形0.060.050.040.030.020.0100200040006000800010000120001400016000

7

排除无效点后语音信号的时域波形0.060.050.040.030.020.010050001000015000

FFT后声音的频谱特性20181614121082005001000150020002500300035004000

排除无效点经FFT后声音的频谱特性20181614121082005001000150020002500300035004000

2. 取8000个和16000个数据进行频谱分析，得到幅度和相位谱，比较二者异同并分析原因。程序如下

x1=wavread('D:\\shengyin',8000); %读取语音信号的数据，赋给变量x1 sound(x1,8000); %播放语音信号

x2=wavread('D:\\shengyin',16000); %读取语音信号的数据，赋给变量x2 sound(x1,8000); %播放语音信号 y1=fft(x1,8192); %对信号做8192点FFT变换 y2=fft(x2,8192);

8

f=fs*(0:4095)/8192; figure(1)

subplot(1,2,1);

plot(x1)

title('取8000数据的语音信号'); xlabel('时间'); ylabel('幅值'); subplot(1,2,2);

plot(x2)

title('取16000数据的语音信号'); xlabel('时间'); ylabel('幅值'); figure(2)

subplot(1,2,1);

plot(angle(y1)); title('相位图') subplot(1,2,2); plot(angle(y2)); title('相位图') figure(3)

subplot(1,2,1)

plot(f,abs(y1(1:4096))); axis([0 4000 0 20]); title('语音信号频谱') xlabel('频率'); ylabel('幅值'); subplot(1,2,2)

plot(f,abs(y2(1:4096))); axis([0 4000 0 20]); title('语音信号频谱')

xlabel('频率'); ylabel('幅值');

取8000数据的语音信号取16000数据的语音信号0.050.050.0450.0450.040.040.0350.0350.030.03值值幅0.025幅0.0250.020.020.0150.0150.010.010.0050.0050050001000000500010000时间时间

9

相位图43210-1-2-3-443210-1-2-3-4相位图05000100000500010000

8000个数据语音信号频谱2018161412201816141216000个数据语音信号频谱幅值108200100200频率300400幅值108200100200频率300400

分析异同点

从两者的时域图看出取8000点数据的语音信号被截断了，而取16000点数据的声音则完好地呈现出来。从频谱图看出幅值大小与fft选用点数N有关，N越大，幅值就相对大一点。

3. 设计数字滤波器和对信号滤波 窗函数法设计FIR滤波器： （1）窗函数设计低通滤波器 程序设计如下： clear;

fs=8000;%采样频率 fprintf('开始录音');

y=wavrecord(2*fs,fs,'double')%录制声音2秒 fprintf('录音结束'); wavplay(y,fs);%播放录音

wavwrite(y,fs,'d:/shengyin3');%声音的存储 [z1,fs,bits]=wavread('d:/shengyin3') Y1=fft(z1,8192);

10

fp=1000;fc=1200;As=100;Ap=1;Fs=7000; wc=2*pi*fc/Fs; wp=2*pi*fp/Fs;% wdel=wc-wp;

beta=10.056;%根据阻带衰减选凯泽窗 N0=ceil(6.4*2*pi/wdel);

N=N0+mod(N0+1,2); %---确保N为奇数 wn= kaiser(N+1,beta); ws=(wp+wc)/2/pi; b=fir1(N,ws,wn); figure(1); freqz(b,1);

x=fftfilt(b,z1); X=fft(x,8192); figure(2);

f=fs*(0:4095)/8192; subplot(3,2,1);

plot(f,abs(Y1(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波前信号频谱');

subplot(3,2,2);plot(f,abs(X(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波后信号频谱'); subplot(3,2,3); plot(angle(Y1));

title('滤波前相位谱'); subplot(3,2,4); plot(angle(X));

title('滤波后相位谱'); subplot(3,2,5); plot(z1);

title('滤波前信号波形'); subplot(3,2,6);plot(x); title('滤波后信号波形');

sound(x,fs,bits); wavwrite(x,fs,'d:/shengyin4.wav');%声音的存储

11

滤波前信号频谱10振幅滤波后信号频谱10振幅50500100020003000频率滤波前相位谱40000100020003000频率滤波后相位谱400050-50.10-0.150-50.10-0.105000滤波前信号波形1000005000滤波后信号波形1000000.511.5x 102400.511.5x 1024

可以看出FIR低通滤波器把1000赫兹以上的频率成分基本滤去，达到信号处理效果。在实际系统中，一个输入信号可以分解为多个正弦信号的叠加，为了使得输出信号不会产生相位失真，必须要求它所包含的这些正弦信号通过系统的时间是一样的。从系统的频率响应来看，就是要求它的相频特性是一条直线。在FIR滤波器的线性相位的性质，这就避免相位失真。 （2）窗函数设计高通滤波器 程序设计如下：

12

clear;close all

wavwrite(y,fs,'d:/shengyin3.Wav');%声音的存储 [z1,fs,bits]=wavread('d:/shengyin3')%读取声音 Y1=fft(z1,8192);

fp=2800;fc=3000;As=100;Ap=1;Fs=7000; wc=2*pi*fc/Fs; wp=2*pi*fp/Fs; wdel=wc-wp; beta=10.056;

N1=ceil(6.4*2*pi/wdel);

N=N1+mod(N1+1,2); %---高通滤波器要确保N为奇数 wn= kaiser(N,beta); ws=(wp+wc)/2/pi;

b=fir1(N-1,ws,'high',wn); figure(1); freqz(b,1);

x=fftfilt(b,z1); X=fft(x,8192); figure(2);

f=fs*(0:4095)/8192; subplot(3,2,1);

plot(f,abs(Y1(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波前信号频谱');

subplot(3,2,2);plot(f,abs(X(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波后信号频谱'); subplot(3,2,3); plot(angle(Y1));

title('滤波前相位谱'); subplot(3,2,4); plot(angle(X));

title('滤波后相位谱'); subplot(3,2,5); plot(z1);

title('滤波前信号波形'); subplot(3,2,6);plot(x); title('滤波后信号波形'); sound(x,fs,bits);

wavwrite(x,fs,'d:/shengyin5.wav');%声音的存储 图形分析如下：

13

滤波前信号频谱1010滤波后信号频谱振幅50振幅0100020003000频率滤波前相位谱4000500100020003000频率滤波后相位谱400050-50.10-0.150-50.010-0.0105000滤波前信号波形1000005000滤波后信号波形1000000.511.5x 102400.511.5x 1024

（3）窗函数设计带通滤波器

[z1,fs,bits]=wavread('d:/shengyin3') Y1=fft(z1,8192);

fp1=1200 ;fp2=3000 ;fc1=1000 ;fc2=3200 ;As=100 ;Ap=1 ;Fs=7000 ; wp1=2*pi*fp1/Fs; wc1=2*pi*fc1/Fs; wp2=2*pi*fp2/Fs; wc2=2*pi*fc2/Fs; wdel=wp1-wc1; beta=10.056;

N1=ceil(6.4*2*pi/wdel);

N=N1+mod(N1+1,2); %---确保N为奇数

14

ws =[(wp1+wc1)/2/pi,(wp2+wc2)/2/pi]; wn= kaiser(N+1,beta); b=fir1(N,ws,wn); figure(1); freqz(b,1)

x=fftfilt(b,z1); X=fft(x,8192);

f=fs*(0:4095)/8192; figure(2);

f=fs*(0:4095)/8192; subplot(3,2,1);

plot(f,abs(Y1(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波前信号频谱');

subplot(3,2,2);plot(f,abs(X(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波后信号频谱'); subplot(3,2,3); plot(angle(Y1));

title('滤波前相位谱'); subplot(3,2,4); plot(angle(X));

title('滤波后相位谱'); subplot(3,2,5); plot(z1);

title('滤波前信号波形'); subplot(3,2,6);plot(x); title('滤波后信号波形'); sound(x,fs,bits);

wavwrite(x,fs,'d:/shengyin8.wav');

15

100)Bd( 0edutinga-100M-20000.10.20.30.40.50.60.70.80.91Normalized Frequency ( rad/sample)5000)se0egred(-5000 esahP-10000-1500000.10.20.30.40.50.60.70.80.91Normalized Frequency ( rad/sample)

滤波前信号频谱滤波后信号频谱1010幅幅振5振5001000200030004000001000200030004000频率频率滤波前相位谱滤波后相位谱5500-50500010000-50500010000滤波前信号波形滤波后信号波形0.10.100-0.100.511.52-0.100.511.52x 104x 104

双线性变换法设计IIR低通滤波器 ①选用butter 程序设计如下： clear;close all

[z1,fs,bits]=wavread('d:/shengyin3') Y1=fft(z1,8192);

16

fp=1000;fc=1200;As=100;Ap=1;Fs=7000; wc=2*fc/Fs; wp=2*fp/Fs;

[N,ws]=buttord(wc,wp,Ap,As); [b,a]=butter(N,ws); figure(1);

[h,w]=freqz(b,a);

plot(w/pi,abs(h));grid; xlabel('\\omega/\\pi'); ylabel('振幅（幅值）');

title('Butterworth型数字低通滤波器的幅频响应'); x=filter(b,a,z1); X=fft(x,8192); figure(2);

f=fs*(0:4095)/8192; subplot(3,2,1);

plot(f,abs(Y1(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波前信号频谱');

subplot(3,2,2);plot(f,abs(X(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波后信号频谱'); subplot(3,2,3); plot(angle(Y1));

title('滤波前相位谱'); subplot(3,2,4); plot(angle(X));

title('滤波后相位谱'); subplot(3,2,5); plot(z1);

title('滤波前信号波形'); subplot(3,2,6);plot(x); title('滤波后信号波形'); sound(x,fs,bits);

wavwrite(x,fs,'d:/shengyin6.wav');图形分析

17

滤波前信号频谱1010滤波后信号频谱振幅50振幅0100020003000频率滤波前相位谱4000500100020003000频率滤波后相位谱400050-50.10-0.150-50.10-0.105000滤波前信号波形1000005000滤波后信号波形1000000.511.5x 102400.511.5x 1024

（5）双线性变换法设计高通滤波器

②选用cheby1 程序设计如下：

[z1,fs,bits]=wavread('d:/shengyin3') Y1=fft(z1,8192);

fp=3000;fc=2800;As=100;Ap=1;Fs=7000; wc=2*fc/Fs; wp=2*fp/Fs;

[N,ws]=cheb1ord(wc,wp,Ap,As); [b,a]=cheby1(N,Ap,wp,'high'); figure(1);

[h,w]=freqz(b,a);

plot(w/pi,abs(h));grid; xlabel('\\omega/\\pi'); ylabel('振幅（幅值）');

title('切比雪夫型数字高通滤波器的幅频响应'); x=filter(b,a,z1); X=fft(x,8192);

18

f=fs*(0:4095)/8192; figure(2);

f=fs*(0:4095)/8192; subplot(3,2,1);

plot(f,abs(Y1(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波前信号频谱');

subplot(3,2,2);plot(f,abs(X(1:4096))); axis([0 4000 0 10]); xlabel('频率'); ylabel('振幅');

title('滤波后信号频谱'); subplot(3,2,3); plot(angle(Y1));

title('滤波前相位谱'); subplot(3,2,4); plot(angle(X));

title('滤波后相位谱'); subplot(3,2,5); plot(z1);

title('滤波前信号波形'); subplot(3,2,6);plot(x); title('滤波后信号波形'); sound(x,fs,bits);

wavwrite(x,fs,'d:/shengyin7.wav');

19

滤波前信号频谱1010滤波后信号频谱振幅50振幅0100020003000频率滤波前相位谱4000500100020003000频率滤波后相位谱400050-50.10-0.150-50.010-0.0105000滤波前信号波形1000005000滤波后信号波形1000000.511.5x 102400.511.5x 1024

以上几种滤波，我们都可以从信号滤波前后的波形图以及频谱图上看出变化。当然，也可以用sound()函数来播放滤波后的语音，从听觉上直接感受语音信号的变化，但由于人耳听力的，有些情况下我们很难听出来的。

五．设计报告思考

1在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。如以上实验过程中，采用双线性变化法设计的butter和cheby1数字滤波器，从图中可以看到这种非线性关系。而窗函数设计法则保持原有的线性相位。

2.分析实际信号频谱过程中，分别用N=8092点N=4096来分析幅度谱。由于采样频率f =8000Hz，故Nyquist频率为4000Hz。整个频谱图是以Nyquist 频率为轴对称的。因此利用fft对信号作频谱分析，只要考察0～Nyquist频率（采样频率一半）范围的幅频特性。比较可见，幅值大小与fft选用点数N有关，只要选的N合适，则不影响结果。但对点数N与信号的相位联系不大明白。 3双线性变换的主要优点：双线性变换不会出现由于高频部分超过折叠频率而产生混叠现象。双线性变换法的缺点：Ω和ω之间的非线性关系得到非线性相位的数字滤波器，使数字滤波器的频响偏移模拟滤波器的频响。窗函数法相位响应有严格的线性,不存在稳定性问题, 设计简单。

4．从性能上来说，IIR滤波器较低的阶数获得高的选择性，所用的存贮单元少，所以经济而效率高。但是这个高效率是以相位的非线性为代价的。选择性越好，则相位非线性越严重。相反，FIR滤波器却可以得到严格的线性相位，然而由于FIR滤波器传输函数的极点固定在原点，所以只能用较高的阶数达到好的选择性；对于同样的滤波器设计指标，FIR滤波器所要求的阶数可以比IIR滤波器高5~10倍，结果，成本较高，信号延时也较大；如果按相同的选择性和相同的线性要求来说，则IIR滤波器就必须加全通网络进行相位较正，同样要大增加滤

20

波器的节数和复杂性。从结构上看IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定。另外，在这种结构中，由于运算过程中对序列的舍入处理，这种有限字长效应有时会引入寄生振荡。相反，FIR滤波器主要采用非递归结构，不论在理论上还是在实际的有限精度运算中都不存在稳定性问题，运算误差也较小。此外，FIR滤波器冲激响应是有限长的，可以采用快速付里叶变换算法，而IIR滤波器不可以，所以运算速度可以快得多。 六．结论

在之前数字信号与处理的学习以及完成课后的作业的过程中，已经使用过MATLAB，对其有了一些基础的了解和认识。通过这次练习是我进一步了解了信号的产生、采样及频谱分析的方法。 以及其中产生信号和绘制信号的基本命令和一些基础编程语言。让我感受到只有在了解课本知识的前提下，才能更好的应用这个工具；并提高对所学知识的综合应用能力，在实践上初步实现对数字信号的处理。并且熟练的应用MATLAB也可以很好的加深我对课程的理解，方便我的思维。这次设计使我了解了MATLAB的使用方法，学会分析滤波器的优劣和性能，提高了分析和动手实践能力。同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助。

七．参考文献：

1.《数字信号处理》丁玉美，高西全等编著，西安：西安电子科技大学出版社 2.《数字信号处理教程》程佩清著，北京：清华大学出版社

3.《数字信号处理——理论、算法与实现（第二版）》胡广书编著，北京：电子工业出版社

4.《数字信号处理（第二版）学习指导书》高西全，丁玉美编著，西安：清华大学出版社

5.《数字信号处理实验指导书（MATLAB版）》孙洪，余翔宇等译，北京：电子工业出版社

6.《基于MATLAB的系统分析与设计——信号处理》楼顺天，李博菡编著，西安：

21