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

实验题目：基于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(t1),f(t1±Δt)，

f(t1±2Δt)，...来表示,只要这些采样点的时间间隔Δt≤1/2F，便可根据各采

样值完全恢复原来的信号f(t)。 这是时域采样定理的一种表述方式。 时域采样定理的另一种表述方式是：当时间信号函数f(t)的最高频率分量为

fM时,f(t)的值可由一系列采样间隔小于或等于1/2fM的采样值来确定,即采样点的重复频率f≥2fM。图为模拟信号和采样样本的示意图。

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

T=T2-T1是信号的持续时间〕，假设其频谱为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数字滤波器的主要程序和图像。 FIR低通滤波器程序见附录1； FIR低通滤波器图像：〔图2—2〕

FIR低通滤波器1.41.210.80.60.40.2000.10.20.30.40.50.60.70.80.91

图2—2 FIR低通滤波器 FIR高通滤波程序见附录2； FIR高通滤波图像：〔图2—3〕

FIR高通滤波器10.80.60.40.203000350040004500500055006000

图2—3 FIR高通滤波器

下面我们将给出IIR数字滤波器的主要程序。 IIR低通滤波器程序见附录3； IIR低通滤波器图像：〔图2—4〕

IIR低通滤波器1.4用butter设计1.210.80.60.40.2005001000150020002500300035004000图2—4 IIR低通滤波器

IIR滤波器高通程序见附录4； IIR滤波器高通图像：〔图2—5〕

IIR高通滤波器1用cheby1设计0.80.60.40.20020004000600080001000012000

图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('滤波后的信号波形');%标题