数字信号处理期末综合实验报告

数字信号处理综合实验报告

实验题目：基于Matlab 的语音信号去噪及仿真

专业名称：

学 号：

姓 名：

日 期：

报告内容：

一、 实验原理

1、去噪的原理

1.1 采样定理

在进行模拟/数字信号的转换过程中，当采样频率fs.max 大于信号中，最高频率fmax 的2倍时，即：fs.max>=2fmax,则采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5～10倍；采样定理又称奈奎斯特定理。 1924年奈奎斯特(Nyquist)就推导出在理想低通信道的最高大码元传输速率的公式: 理想低通信道的最高大码元传输速率=2W*log2 N (其中W 是理想低通信道的带宽,N 是电平强度) 为什么把采样频率设为8kHz? 在数字通信中，根据采样定理, 最小采样频率为语音信号最高频率的2倍

频带为F 的连续信号 f (t ) 可用一系列离散的采样值f (t 1), f (t 1±Δt ) ，f (t 1±2Δt ) ，... 来表示, 只要这些采样点的时间间隔Δt ≤1/2F ，便可根据各采

样值完全恢复原来的信号f (t ) 。 这是时域采样定理的一种表述方式。 时域采样定理的另一种表述方式是：当时间信号函数f (t ) 的最高频率分量为fM 时, f (t ) 的值可由一系列采样间隔小于或等于1/2fM 的采样值来确定, 即采样点的重复频率f ≥2fM 。图为模拟信号和采样样本的示意图。

时域采样定理是采样误差理论、随机变量采样理论和多变量采样理论的基础。对于时间上受限制的连续信号f (t ) （即当│t │>T 时, f (t )=0,

这里T =T 2-T 1是信号的持续时间），若其频谱为F （ω）, 则可在频域上用一系列

离散的采样值

（1-1）

采样值来表示, 只要这些采样点的频率间隔

（1-2）

。

1.2 采样频率 采样频率，也称为采样速度或者采样率，定义了每秒从连续信号中提取并组成离散信号的采样个数，它用赫兹（Hz ）来表示。采样频率的倒数是采样周期或者叫作采样时间，它是采样之间的时间间隔。通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。

采样频率只能用于周期性采样的采样器，对于非周期性采样的采样器没有规则限制。 采样频率的常用的表示符号是 fs。 通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。采样频率与声音频率之间有一定的关系，根据采样定理，只有采样频率高于声音信号最高频率的两倍时，才能把数字信号表示的声音还原成为原来的声音。这就是说采样频率是衡量声卡采集、记录和还原声音文件的质量标准。

采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。每增加一个采样位数相当于力度范围增加了6dB 。采样位数越多则捕捉到的信号越精确。对于采样率来说你可以想象它类似于一个照相机，44.1kHz 意味着音频流进入计算机时计算机每秒会对其拍照达441000次。显然采样率越高，计算机摄取的图片越多，对于原始音频的还原也越加精确

2、数字滤波器设计的基本原理

1.1 FIR数字滤波器的设计 FIR ：有限脉冲响应滤波器。有限说明其脉冲响应是有限的。与IIR 相比，它具有线性相位、容易设计的优点。这也就说明，IIR 滤波器具有相位不线性，不容易设计的缺点。而另一方面，IIR 却拥有FIR 所不具有的缺点，那就是设计同样参数的滤波器，FIR 比IIR 需要更多的参数。这也就说明，要增加DSP 的计算量。DSP 需要更多的计算时间，对DSP 的实时性有影响。FIR 滤波器的设计比较简单，就是要设计一个数字滤波器去逼近一个理想的低通滤波器。通常这个理想的低通滤波器在频域上是一个矩形窗。根据傅里叶变换我们可以知道，此函数在时域上是一个采样函数。通常此函数的表达式为：

sa （n ）＝sin （n ）/n（1-3）

但是这个采样序列是无限的，计算机是无法对它进行计算的。故我们需要对此采样函数进行截断处理。也就是加一个窗函数。就是传说中的加窗。也就是把这个时域采样序列去乘一个窗函数，就把这个无限的时域采样序列截成了有限个序列值。但是加窗后对此采样序列的频域也产生了影响：此时的频域便不在是一个理想的矩形窗，而是成了一个有过渡带，阻带有波动的低通滤波器。通常根据所加的窗函数的不同，对采样信号加窗后，在频域所得的低通滤波器的阻带衰减也不同。通常我们就是根据此阻带衰减去选择一个合适的窗函数。如矩形窗、汉宁窗、汉明窗、BLACKMAN 窗、凯撒窗等。

1.2 IIR数字滤波器的设计

对于数字高通、带通滤波器的设计，通用方法为双线性变换法。可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。具体设计步骤如下：

（1）确定所需类型数字滤波器的技术指标。

（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频

率，转换公式为Ω=2/T tan(0.5ω) （1-4）

(3)将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。

（4）设计模拟低通滤波器。

（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。

（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。

我们知道，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。为了克服之一缺点，可以采用双线性变换法。

二、 实验步骤

其大概流程框图可如下表示：（图2-1）

图2-1

1. 滤波器的设计。

利用窗函数法设计FIR 滤波器的步骤。如下：

（1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N 。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。

（2）构造希望逼近的频率响应函数。

（3）计算h(n).。

（4）加窗得到设计结果。

接下来，我们根据语音信号的特点给出有关滤波器的技术指标：

低通滤波器的性能指标：

fp=1000Hz，fc=1200Hz，As=50db ,Ap=1dB

高通滤波器的性能指标：

fp=3500Hz，fc=4000Hz，As=50dB，Ap=1dB

在Matlab 中，可以利用函数fir1设计FIR 滤波器，利用Matlab 中的函数freqz 画出各步步器的频率响应。

MATLAB 信号处理工具箱函数cheblap,cheblord 和cheeby1是切比雪夫I 型滤波器设计函数。我们用到的是cheeby1函数，其调用格式如下：

[B,A]=cheby1(N,Rp,wpo,’ftypr ’)

[B,A]=cheby1(N,Rp,wpo,’ftypr ’, ’s ’)

利用模拟滤波器设计IIR 数字低通滤波器的步骤：如下：

（1）确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。

（2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。

（3）按照模拟低通滤波器的技术指标设计及过渡模拟低通滤波器。

（4）用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。

MATLAB 信号处理工具箱函数cheblap,cheblord 和cheeby1是切比雪夫I 型滤波器设计函数。我们用到的是cheeby1函数，其调用格式如下：

[B,A]=cheby1(N,Rp,wpo,’ftypr ’)

[B,A]=cheby1(N,Rp,wpo,’ftypr ’, ’s ’)

函数butter,cheby1和ellip 设计IIR 滤波器时都是默认的双线性变换法，所以 在设计滤波器时只需要代入相应的实现函数即可。

2. 语音信号的录制。

单击自己的电脑开始程序，选择所有程序，接着选择附件，再选择录音。自己录入语音信号，然后保存在MATLAB 文件夹里面，命名为“chenghaijie1.wav ”。

3. 在MATLAB 平台上读入语音信号、绘制频谱图并回放原始语音信号。 利用MATLAB 中的wavread 命令来读入（采集）语音信号，将它赋值给某一向量。 [y,fs,bits]=wavread(' [N1 N2]);用于读取语音，采样值放在向量y 中，fs 表示采样频率(Hz)，bits 表示采样位数。[N1 N2]表示读取从N1点到N2点的值（若只有一个N 的点则表示读取前N 点的采样值。

4. 利用matlab 函数randn 编程加入一段随机噪音信号，再利用设计的FIR 和IIR 滤波器去噪，分别绘制去噪后的频谱图、回放语音信号与原始信号的频谱图、原始语音信号比较，并且比较两种滤波器的优缺点和得出语音信号的频段。

三、 实验结果

下面我们将给出设计FIR 数字滤波器的主要程序和图像。

图2—2 FIR 低通滤波器

FIR 高通滤波程序见附录2；

FIR 高通滤波图像：（图2—3）

FIR 高通滤波器

1

0.8

0.6

0.4

0.2

图2—3 FIR 高通滤波器

下面我们将给出IIR 数字滤波器的主要程序。

IIR 低通滤波器程序见附录3；

IIR 低通滤波器图像：（图2—4）

IIR 低通滤波器

1.4

1.2

1

0.8

0.6

0.4

0.2

0图2—4 IIR低通滤波器

IIR 滤波器高通程序见附录4；

IIR 滤波器高通图像：（图2—5）

IIR 高通滤波器

图2—5 IIR高通滤波器

原始语音信号的回放及频谱分析程序：见附录5；

原始语音信号及其频谱分析图像、运行程序可以听到原始信号的回放：（图2—6）

图2—6 原始语音信号及其频谱分析图像

利用randn 函数把一段随机噪音信号加入原始语音信号的信号处理过程程序：见附录6； 加噪后语音信号的时域波形、频谱图：

（图2—7）

（图2—7）加噪后语音信号的时域波形、频谱图

利用FIR 低通滤波器法去噪程序：见附录7；

FIR 滤波器去噪前与去噪后语音信号的时域波形、频谱图：（图2—8）

（图2—8）FIR 滤波器去噪前与去噪后语音信号的时域波形、频谱图 利用IIR 低通滤波器法去噪程序：见附录8；

IIR 滤波器去噪前与去噪后语音信号的时域波形、频谱图：（图2—9）

（图2—9）IIR 滤波器去噪前与去噪后语音信号的时域波形、频谱图

四、 实验分析

1. 由（图2—8）FIR 滤波器去噪前与去噪后语音信号的时域波形、频谱图和（图2—9）IIR 滤波器去噪前与去噪后语音信号的时域波形、频谱图可看出原始语音信号和加噪语音信号时域波形和频谱图的区别。加噪后的语音信号的时域波形比原始语音信号要模糊得多，频谱图则是在频率5000Hz 以后出现了明显的变化。

再通过滤波前的信号波形和频谱图的对比，可以明显看出滤波后的波形开始变得清晰了，有点接近原始信号的波形图了。滤波后信号的频谱图也在5000Hz 以后开始逐渐接近原始语音信号的频谱图。

再从对语音信号的回放，人耳可以明显辨别出加噪后的语音信号比较浑浊，还有很明显嘎吱嘎吱的杂音在里面。滤波后，语音信号较加噪后的信号有了明显的改善，基本可以听清楚了，而且杂音也没有那么强烈，但是声音依然没有原始语音信号那么清晰脆耳。

2. 结合去噪后的频谱图对比可知FIR 滤波器和IIR 滤波器的优缺点

IIR 数字滤波器采用递归型结构，即结构上带有反馈环路。IIR 滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。

（1）IIR 数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。FIR 滤波器则要求较低。

（2）IIR 滤波器运算误差大，有可能出现极限环振荡，FIR 相比之下运算误差较小，不会出现极限环振荡。

（3）IIR 幅频特性精度很高，不是线性相位的，可以应用于对相位信息不敏感的音频信号上；

（4）与FIR 滤波器的设计不同，IIR 滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。在MATLAB 下设计不同类型IIR 滤波器均有与之对应的函数用于阶数的选择。

（5）IIR 单位响应为无限脉冲序列FIR 单位响应为有限的

（6）FIR 幅频特性精度较之于iir 低，但是线性相位，就是不同频率分量的信号经过FIR 滤波器后他们的时间差不变。这是很好的性质。

（7）IIR 滤波器有噪声反馈，而且噪声较大，FIR 滤波器噪声较小。

FIR 幅频特性精度较之于IIR 低，但是线性相位，就是不同频率分量的信号经过FIR 滤波器后他们的时间差不变，这是很好的性质。

3. 通过对语音信号做频谱分析和通过比较加噪前后，语音的频谱和语音回放，能明显的感觉到加入噪声后回放的声音与原始的语音信号有很大的不同, 前者随较尖锐的干扰啸叫声。从含噪语音信号的频谱图中可以看出含噪声的语音信号频谱, 在整个频域范围内分是布均匀。其实，这正是干扰所造成的。通过滤波前后的对比，低通滤波后效果最好，高通滤波后的效果最差。由此可见，语音信号主要分布在低频段，而噪声主要分布在高频段。

五、 实验总结

基于Matlab 的语音信号去噪及仿真课题的特色在于它将语音信号看作一个向

量，于是就把语音数字化了。那么，就可以完全利用数字信号处理的知识来解决语音及加噪处理问题。我们可以像给一般信号做频谱分析一样，来对语音信号做频谱分析，也可以较容易的用数字滤波器来对语音进行滤波处理。通过比较加噪前后，语音的频谱和语音回放，能明显的感觉到加入噪声后回放的声音与原始的语音信号有很大的不同, 前者随较尖锐的干扰啸叫声。从含噪语音信号的频谱图中可以看出含噪声的语音信号频谱, 在整个频域范围内分是布均匀。其实，这正是干扰所造成的。通过两种滤波器滤波前后的对比，可知FIR 滤波器和IIR 滤波器的优缺点，低通滤波后效果最好，高通滤波后的效果最差。由此可见，语音信号主要分布在低频段，而噪声主要分布在高频段。

六、 参考文献

[1] Sanjit K.Mitra，数字信号处理实验指导书，电子工业出版社，2005 年 1月

[2] 胡航，语音信号处理，哈尔滨工业大学出版社，2000 年 5 月

[3］姚天任. 数字语音处理[M].武汉：华中科技大学出版社，2005 年 3 月

七、 附录

附录 1

Ft=8000;%设置抽样频率为8000Hz

Fp=1000;%设置通带边界为1000Hz

Fs=1200;%设置阻带边界为1200Hz

wp=2*Fp/Ft;%计算归一化通带边界角频率

ws=2*Fs/Ft;%计算归一化阻带边界角频率

rp=1;%设置通带纹波为1dB

rs=50;%设置阻带衰减为50dB

p=1-10.^(-rp/20);%通带阻带波纹

s=10.^(-rs/20);

fpts=[wp ws];

mag=[1 0];

dev=[p s];

[n21,wn21,beta,ftype]=kaiserord(fpts,mag,dev);%kaiserord求阶数截止频率

b21=fir1(n21,wn21,Kaiser(n21+1,beta));%由fir1设计滤波器

[h,w]=freqz(b21,1);%得到频率响应

plot(w/pi,abs(h));%绘制FIR 低通滤波器

title('FIR低通滤波器');%标题

附录2

Ft=8001;%设置抽样频率为8001Hz

Fp=4000;%设置通带边界为4000Hz

Fs=3500;%设置阻带边界为3500Hz

wp=2*Fp/Ft;%计算归一化通带边界角频率

ws=2*Fs/Ft;%计算归一化阻带边界角频率

rp=1;%设置通带纹波为1dB

rs=50;%设置阻带衰减为50dB

p=1-10.^(-rp/20);%通带阻带波纹

s=10.^(-rs/20);

fpts=[ws wp];

mag=[0 1];

dev=[p s];

[n23,wn23,beta,ftype]=kaiserord(fpts,mag,dev);%kaiserord求阶数截止频率

b23=fir1(n23,wn23,'high',Kaiser(n23+1,beta));%由fir1设计滤波器

[h,w]=freqz(b23,1);%得到频率响应

plot(w*12000*0.5/pi,abs(h));%绘制FIR 高通滤波器

title('FIR高通滤波器');%标题

axis([3000 6000 0 1.2]);%确定横纵坐标的范围

附录3

Ft=8000;%设置抽样频率为8000Hz

Fp=1000;%设置通带边界为1000Hz

Fs=1200;%设置阻带边界为1200Hz

wp=2*pi*Fp/Ft;%计算归一化通带边界角频率

ws=2*pi*Fs/Ft;%计算归一化阻带边界角频率

fp=2*Ft*tan(wp/2);%计算通带边界频率

fs=2*Fs*tan(wp/2);%计算阻带边界频率

[n11,wn11]=buttord(wp,ws,1,50,'s');%求低通滤波器的阶数和截止频率

[b11,a11]=butter(n11,wn11,'s');%利用butter 求S 域的频率响应的参数

[num11,den11]=bilinear(b11,a11,0.5);%双线性变换实现S 域到Z 域的变换

[h,w]=freqz(num11,den11);%根据参数求出频率响应

plot(w*8000*0.5/pi,abs(h));%绘制IIR 低通滤波器

title('IIR低通滤波器');%标题

legend('用butter 设计');%添加图例的标注

grid;%使绘制的图有网格

附录4

Ft=8001;%设置抽样频率为8001Hz

Fp=4000;%设置通带边界为4000Hz

Fs=3500;%设置阻带边界为3500Hz

wp1=tan(pi*Fp/Ft);%通带到低通滤波器参数转换

ws1=tan(pi*Fs/Ft);%阻带到低通滤波器参数转换

wp=1;%计算归一化通带边界角频率

ws=wp1*wp/ws1;%计算归一化阻带边界角频率

[n13,wn13]=cheb1ord(wp,ws,1,50,'s');%求模拟的低通滤波器阶数和截止频

[b13,a13]=cheby1(n13,1,wn13,'s');%利用cheby1求S 域的频率响应的参数

[num,den]=lp2hp(b13,a13,wn13); %将S 域低通参数转为高通的

[num13,den13]=bilinear(num,den,0.5); %利用双线性变换实现S 域到Z 域转

[h,w]=freqz(num13,den13);%得到频率响应

plot(w*21000*0.5/pi,abs(h));%绘制IIR 高通滤波器

title('IIR高通滤波器');%标题

legend('用cheby1设计');%添加图例的标注

附录5

[x,fs,bits]=wavread('chenghaijie1.wav');%利用MATLAB 中的wavread 命令来读入（采集）语音信号，将它赋值给某一向量

sound(x,fs,bits);%播放语音信号

X=fft(x,100000);%对信号做50000点FFT 变换

magX=abs(X);%求幅度

angX=angle(X);%求相位

plot(x);title('原始信号波形');%绘制原始信号的波形

plot(X); title('原始语音信号采样后的频谱图')%绘制原始语音信号采样后的频谱图

附录6

[y,fs,bits]=wavread('chenghaijie1.wav');%利用MATLAB 中的wavread 命令来读入（采集）语音信号，将它赋值给某一向量

sound(y,fs)%播放语音信号

n=length(y)%计算向量y 的长度

y_p=fft(y,n);%对向量y 进行fft 变换

f=fs*(0:n/2-1)/n;%对应点的频率

figure(1)%建立第一幅图

subplot(2,1,1);%分配绘制的位置为第一块位置

plot(y);%绘制原始语音信号采样后时域波形

title('原始语音信号采样后时域波形');%标题

xlabel('时间轴')%标注X 轴

ylabel('幅值 A')%标注Y 轴

subplot(2,1,2);%分配绘制的位置为第二块位置

plot(f,abs(y_p(1:n/2)));%绘制原始语音信号采样后的频谱图

title('原始语音信号采样后的频谱图');

xlabel('频率Hz');%标注X 轴

ylabel('频率幅值');%标注Y 轴

L=length(y)%计算向量y 的长度

noise=0.1*randn(L,2);%设置加入的随机噪声函数

y_z=y+noise;%将原始信号与随机噪声函数相叠加

sound(y_z,fs)%试听加入随机噪声函数的信号

n=length(y);%计算向量y 的长度

y_zp=fft(y_z,L);%对加入正弦噪声向量y 进行fft 变换

f=fs*(0:L/2-1)/L;%计算对应点的频率

figure(2)%建立第二幅图

subplot(2,1,1);%分配绘制的位置为第一块位置

plot(y_z);%绘制加噪语音信号时域波形

title('加噪语音信号时域波形');%标题

xlabel('时间轴')%标注X 轴

ylabel('幅值 A')%标注Y 轴

subplot(2,1,2);%分配绘制的位置为第二块位置

plot(f,abs(y_zp(1:L/2)));%绘制加噪语音信号频谱图

title('加噪语音信号频谱图');%标题

xlabel('频率Hz');%标注X 轴

ylabel('频率幅值');%标注Y 轴

附录7

[y,fs,bits]=wavread('chenghaijie1.wav');%利用MATLAB 中的wavread 命令来读入（采集）语音信号，将它赋值给某一向量

sound(y,fs)%播放语音信号

n=length(y)%计算向量y 的长度

y_p=fft(y,n);%对向量y 进行fft 变换

f=fs*(0:n/2-1)/n;%对应点的频率

figure(1)%建立第一幅图

subplot(2,1,1);%分配绘制的位置为第一块位置

plot(y);%绘制原始语音信号采样后时域波形

title('原始语音信号采样后时域波形');%标题

xlabel('时间轴')%标注X 轴

ylabel('幅值 A')%标注Y 轴

subplot(2,1,2);%分配绘制的位置为第二块位置

plot(f,abs(y_p(1:n/2)));%绘制原始语音信号采样后的频谱图

title('原始语音信号采样后的频谱图');

xlabel('频率Hz');%标注X 轴

ylabel('频率幅值');%标注Y 轴

L=length(y)%计算向量y 的长度

noise=0.1*randn(L,2);%设置加入的随机噪声函数

y_z=y+noise;%将原始信号与随机噪声函数相叠加

sound(y_z,fs)%试听加入随机噪声函数的信号

n=length(y);%计算向量y 的长度

y_zp=fft(y_z,L);%对加入正弦噪声向量y 进行fft 变换

f=fs*(0:L/2-1)/L;%计算对应点的频率

figure(2)%建立第二幅图

subplot(2,1,1);%分配绘制的位置为第一块位置

plot(y_z);%绘制加噪语音信号时域波形

title('加噪语音信号时域波形');%标题

xlabel('时间轴')%标注X 轴

ylabel('幅值 A')%标注Y 轴

subplot(2,1,2);%分配绘制的位置为第二块位置

plot(f,abs(y_zp(1:L/2)));%绘制加噪语音信号频谱图

title('加噪语音信号频谱图');%标题

xlabel('频率Hz');%标注X 轴

ylabel('频率幅值');%标注Y 轴

Ft=5000;%设置抽样频率为5000Hz

Fp=1000;%设置通带边界为1000Hz

Fs=1200;%设置阻带边界为1200Hz

wp=2*Fp/Ft;%计算归一化通带边界角频率

ws=2*Fs/Ft;%计算归一化阻带边界角频率

rp=1;%设置通带纹波为1dB

rs=50;%设置阻带衰减为50dB

p=1-10.^(-rp/20);%通带阻带波纹

s=10.^(-rs/20);

fpts=[wp ws];

mag=[1 0];

dev=[p s];

[n21,wn21,beta,ftype]=kaiserord(fpts,mag,dev);%kaiserord求阶数截止频率

b21=fir1(n21,wn21,Kaiser(n21+1,beta));%由fir1设计滤波器

[h,w]=freqz(b21,1);%得到频率响应

plot(w/pi,abs(h));%绘制FIR 低通滤波器

title('FIR低通滤波器');%标题

x=fftfilt(b21,y_z);%重叠相加法实现线性卷积的计算

X=fft(x,n);%计算x 的fft

figure(3);%建立第三幅图

subplot(2,2,1);plot(f,abs(y_zp(1:n/2)));%绘制滤波前信号的频谱 title('滤波前信号的频谱');%标题

subplot(2,2,2);plot(f,abs(X(1:n/2)));%绘制滤波后信号的频谱

title('滤波后信号的频谱');%标题

subplot(2,2,3);plot(y_z);%绘制滤波前信号的时域波形

title('滤波前信号的时域波形')%标题

subplot(2,2,4);plot(x);%绘制滤波后信号的时域波形

title('滤波后信号的时域波形')%标题

sound(x,fs,bits)%重新听取滤噪后的信号

附录8

Ft=8000;%设置抽样频率为8000Hz

Fp=1000;%设置通带边界为1000Hz

Fs=1200;%设置阻带边界为1200Hz

wp=2*pi*Fp/Ft;%计算归一化通带边界角频率

ws=2*pi*Fs/Ft;%计算归一化阻带边界角频率

fp=2*Ft*tan(wp/2);%计算通带边界频率

fs=2*Fs*tan(wp/2);%计算阻带边界频率

[n11,wn11]=buttord(wp,ws,1,50,'s');%求低通滤波器的阶数和截止频率

[b11,a11]=butter(n11,wn11,'s');%利用butter 求S 域的频率响应的参数

[num11,den11]=bilinear(b11,a11,0.5);%双线性变换实现S 域到Z 域的变换

[h,w]=freqz(num11,den11);%根据参数求出频率响应

plot(w*8000*0.5/pi,abs(h));%绘制IIR 低通滤波器

title('IIR低通滤波器');%标题

legend('用butter 设计');%添加图例的标注

grid;%使绘制的图有网格

[y,fs,nbits]=wavread ('chenghaijie1.wav');%利用MATLAB 中的wavread 命令来读入（采集）语音信号，将它赋值给某一向量

n = length (y) ;%求出语音信号的长度

noise=0.01*randn(n,2);%随机函数产生噪声 s=y+noise;%语音信号加入噪声

S=fft(s); %傅里叶变换

z11=filter(num11,den11,s);%滤噪声 sound(z11);%听取滤噪声后的信号

m11=fft(z11);%求滤波后的信号

subplot(2,2,1);%分配绘制的位置为第一块位置 plot(abs(S),'g');%绘制滤波前信号的频谱 title('滤波前信号的频谱');%标题

grid;%网格

subplot(2,2,2);%分配绘制的位置为第二块位置 plot(abs(m11),'r');%绘制滤波后信号的频谱 title('滤波后信号的频谱');%标题

grid;%网格

subplot(2,2,3);%分配绘制的位置为第三块位置 plot(s);%绘制滤波前信号的波形

title('滤波前信号的波形');%标题

grid;%网格

subplot(2,2,4);%分配绘制的位置为第四块位置 plot(z11);%绘制滤波后的信号波形

title('滤波后的信号波形');%标题