数字信号处理报告

《数字信号处理》

课

程

设

计

报

告

设计一信号的基本运算

【一】设计目的

熟悉信号的基本运算，通过仿真结果加深理解。

【二】设计内容

1.现有一段以采样频率8192Hz采样得到的语音数据x(k)，为了模仿出回音的效果，可以对数据进行如下处理：x(k)=x(k)+a*x(k-d)，其中d为时延，a为时延信号的衰减幅度。根据上述原理，进行仿真实验，并演示回声的效果。(参数：时延d=0.4秒，衰减幅度a=0.5)

2.对上述语音信号进行时间反转x(-k)、上采样x(2k)、下采样x(k/2)操作，并演示这些运算的效果。

【三】设计过程，调试，分析。

数字信号处理的主要研究对象是语音信号和图像信号,语音信号的研究可以从时域和

频域两个方面来进行。其中时域的分析处理有两种方法: 一种是进行语音信号分析, 这属于线性处理的范畴,主要是通过信号的加减、时移、倍乘、卷积、求相关函数等来实现; 另一种是生成和变换成各种调制信号, 这属于非线性的范畴,主要是对信号平均累加器的动

态范围进行压缩扩张, 用门限方法对噪声的抑制。对频域分析处理, 即对信号的频率特性在频谱中加以分析研究, 这拓展了信号分析的范围, 是对不确定信号分析的主要方法。在实际应用中, 信号的时域频分析经常同时进行。

1 语音信号的采集与时频分析

1. 1 语音信号的采集

利用 PC 机上的声卡和WINDOWS操作系统可以进行数字信号的采集。将话筒输入计算机的语音输入插口上, 启动录音机。按下录音按钮, 接着对话筒说话/大家好0 ,说完后停止录音, 屏幕左侧将显示所录声音的长度。点击放音按钮,可以实现所录音的重现。以文件名/ dajiahao0 保存入 g:\MATLAB\ work中。可以看到, 文件存储器的后缀默认wav。

在 Matlab 软件平台下可以利用函数wavread对语音信号进行采样,得到了声音数据变量x,同时把x的采样频率Fs= 8000Hz和数据位Nbits= 8Bit放进了MATALB的工作空间。其程序如下

Fs= 8000Hz; % 给出抽样频率

Nbits= 8Bit; % 放入数据位

[ x, fs, Nbits] = wavread(- g: \ Matlab \ wor k \ chaotian.wav ) ; % 把语音信号进行加载入 Matlab 仿真软件平台中

sound( x, fs) ; % 对加载的语音信号进行回放

stem( x) ; title(‘ 语音信号的时域波形’) ;

1. 2 语音信号的时域分析

因录音采样频率无法准确。故用余弦序列表示信号的迟延，反转，上采样下采样。

迟延：

题目要求时延d为0.4s.取点数k为0.4*Fs。产生一组新序列将原序列依次推后k点。

【四】程序清单与运行结果

程序清单：

y=load('E:\erf1s1t0');

fs=8192;

n=length(y);

m=n+round(0.4*fs);

q=(m-1)/fs;

l=round(0.4*fs);

z=zeros(1,m);

x=z;

for i=1:n x(i)=y(i);end %原信号

for i=1:n z(l+i)=y(i);end %信号的延拖

o=x+0.4*z; %回声信号

sound(o) %声音

t=0:1/fs:q ;

subplot(3,1,1);

plot(t,x)

title('原信号');

subplot(3,1,2);

plot(t,z)

title('时延信号');

subplot(3,1,3);

plot(t,o)

title('回声信号');

运行结果：

翻转：

K为序列长度，将序列第i个点值依次与第k-i点值进行调换产生新序列，即为翻转序列。

采样：

下采样，也就是对信号的抽取。其实，上采样和下采样都是对数字信号进行重采，重采的采样率与原来获得该数字信号的采样率比较，大于原信号的称为上采样，小于的则称为下采样。上采样的实质也就是内插或插值。

程序清单：

%对上述语音信号进行时间反转x(-k)、上采样x(2k)、下采样x(k/2)操作

data=load('E:\erf1s1t0');

x=data; %原声信号

d=0.4; %d为时延

Fs=8192; %Fs为采样频率

k=0.4*8192;

a=0.5; %a为时延信号的衰减幅度

z=linspace(0,0,k); %linspace为线性间隔的向量

y=[z x'];

x=[x' z];

f=x+a*y; %回声信号

f1=fliplr(x); %翻转信号

f2=downsample(x,2); %下采样信号

f3=upsample(x,2); %上采样信号

subplot(4,1,1);

plot(x);title('原声信号')

subplot(4,1,2);

plot(f1);title('翻转信号')

subplot(4,1,3);

plot(f2);title('下采样')

subplot(4,1,4);

plot(f3);title('上采样')

sound(f1)

运行结果：

设计二正余弦信号的谱分析

【一】设计目的

用DFT实现对正余弦信号的谱分析；

观察DFT长度和窗函数长度对频谱的影响；

对DFT 进行谱分析中的误差现象获得感性认识。 【二】设计内容

1.对一个频率为10Hz ，采样频率为64Hz 的32点余弦序列进行谱分析，画出其频谱图；若将频率改为11Hz ，其他参数不变，重新画出该序列的频谱图，观察频率泄漏现象，分析原因。

2.考察DFT 的长度对双频率信号频谱分析的影响。设待分析的信号为

150)

2sin()2sin(5.0)(21≤≤+=n n f n f n x ππ （1.2）

令两个长度为16的正余弦序列的数字频率为22.01=f 及34.02=f 。取N 为四个不同值16，32，64，128。画出四个DFT 幅频图，分析DFT 长度对频谱分辨率的影响。

3.在上题中若把两个正弦波的频率取得较近，令22.01=f ，25.02=f ，试问怎样选择FFT 参数才能在频谱分析中分辨出这两个分量？ 【三】设计过程，调试，分析。

连续信号的频谱分析在工程上有着广泛的应用，计算机的第一步就是把连续信号离散化，一是采样，二是截断。由此会产生频率混叠误差、栅栏效应和截断误差。

（一） 对一个频率为10Hz ，和11Hz ，采样频率为64Hz 的32点余弦序列用Matlab 计算它的DFT

程序为

N=32;Fs=64;f=10; n=0:N-1;

x=cos(2*pi*10*n/64) y=cos(2*pi*11*n/64) subplot(2,1,1) X=abs(fft(x,N)) stem(n,X,'.');

xlabel('n');title('10HZ 32点变化幅频曲线') subplot(2,1,2) Y=abs(fft(y,N)) stem(n,Y,'.');

xlabel('n');title('11HZ 32点变化幅频曲线')

运行结果：

通过图可看出：频率为10Hz的余弦曲线DFT只有两个点不等于零，位于k=5和k=27处，k=5对应于频率10Hz，k=27对应于频率54Hz(也就是-10Hz)。这样DFT确实正确的分辨了余弦信号的频率。但是这样的理想结果是恰巧得到的，此时我们借去了五个完整的余弦周期

（f*N/Fs=5）.

将频率改为11Hz，采样频率和窗长度依然为32点，计算图像可看出：频谱图上k=5和k=27处都有较大的峰值，而其它的点上幅度不再为零。这两个峰值对应的频率为10Hz和12Hz，所以，信号的峰值位于两者之间，本来是单一的11Hz频率的能量会分不到许多DFT频率上，这种现象叫频率泄露，来源于截断效应。

（二）对双信号的频谱分析程序为：

16点，32点，64点，128点的频谱分析程序：

N=16;n=0:15;

f1=0.22;f2=0.34;

x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n);

subplot(4,2,1),stem(n,x,'.');

%xlable('n'); ylable('x1(n)');

title('余弦序列');

X=abs(fft(x,N));

subplot(4,2,2)

k=0:N-1;

stem(k,X,'.')

%xlable('k');ylable('X(k)')

string=[num2str(N),'点FFT幅频曲线']; title(string);

N=32;n=0:15;

f1=0.22;f2=0.34;

x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n); subplot(4,2,3),stem(n,x,'.')

%xlable('n'); ylable('x1(n)');

title('余弦序列');

X=abs(fft(x,N));

subplot(4,2,4)

k=0:N-1;

stem(k,X,'.')

%xlable('k');ylable('X(k)')

string=[num2str(N),'点FFT幅频曲线']; title(string);

N=64;n=0:15;

f1=0.22;f2=0.34;

x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n); subplot(4,2,5),stem(n,x,'.')

%xlable('n'); ylable('x1(n)');

title('余弦序列');

X=abs(fft(x,N));

subplot(4,2,6)

k=0:N-1;

stem(k,X,'.')

%xlable('k');ylable('X(k)')

string=[num2str(64),'点FFT幅频曲线']; title(string);

N=128;n=0:15;

f1=0.22;f2=0.34;

x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n); subplot(4,2,7),stem(n,x,'.')

%xlable('n'); ylable('x1(n)');

title('余弦序列');

X=abs(fft(x,N));

subplot(4,2,8)

k=0:N-1;

stem(k,X,'.')

%xlable('k');ylable('X(k)')

string=[num2str(128),'点FFT幅频曲线'];

title(string); 运行结果:

分析：DFT 样本值就是其DTFT 在相应位置的采样。

在图一中很难看出两个峰值，因此要提高它的分辨率，故把N 增大，逐渐可以看出它有两个峰值，将k 换算成数字频率f=w/2*pi=k/N.这样可确定峰值的位置大体在f=0.21和0.35之附近，与信号的给定频率有一定的误差，这也是截断和泄露带来的问题，在这图上还可以看到一些较小的峰，这是很难判断是输入信号固有的，还是由泄露引起的。

说明了增加DFT 长度N 减小了相邻样本间的频率间距，提高频谱的视在分辨率，因而可以提高样本未知的测定精度。

（三） 把两个正弦波的频率取得较近，另fr=[0.22,0.25]应怎样选择FFT 参数

要能分清两个频率，分辨率至少应达到f=0.03.因为此处的数字频率是对采样频率Fs 进行归一化后的，几fr 最大为1，因此总的样本数至少要达到1/0.03=33

加窗以后可以使频谱函数更加光滑，便于分辨峰值位置和准确的数值，为了提高实际分辨率，应该尽量增加信号的长度n 及DF 长度N,当受到条件限制不能提高n ，则单独提高N 可以提高视在分辨率。

当22.01=f ，25.02=f 时，16点，32点，64点，128点的频谱运行结果：

（四）谱分析中的误差现象

①混叠误差。产生混叠误差的原因是：信号的离散化是通过抽样实现的，而抽样频率再高总是有限的。除带限信号外，如果信号的最高频率Ω趋于无穷，则实际器件无法满足抽样定理。而抽样过程如果不满足抽样定理，就会产生频谱的混叠，即混叠误差。要减少或避免混叠误差，应该提高抽样频率，以设法满足抽样定理，或者采用抗混叠滤波这样的信号预处理措施。

②栅栏效应。对于非周期信号来说，理论上应具有连续的频谱，但数字谱分析是用DFT 来近似的，是用频谱的抽样值逼近连续频谱值，只能观察到有限个频谱值，每一个间隔中的频谱都观察不到。如同通过“栅栏”观察景物一样，一部分被“栅栏”所阻挡，看不见，把这种现象称为栅栏效应。连续时间信号只要采用数字谱分析的方法，就必定产生栅栏效应，栅栏效应只能减小而无法避免。能够感受的频谱最小间隔值，称为频谱分辨率，而频率分辨率与信号截断长度成反比，且频率分辨率越小，插值效果越好。因此即使连续信号是有限时间的，也应选择一个大的截断长度来改善栅栏效应。

③截断误差。截断误差是由于对信号进行截断，把无限长的信号限定为有限长，即令有限区间以外的函数值均为零值的近似处理而产生的，这种处理相当于用一个矩形(窗)信号乘待分析的连续时间信号。必然会引起吉布斯效应（波动），也会把窗函数的频谱引入信号频谱 ,造成混叠，所以需要考虑其误差的问题。

设计三数字滤波器的设计及实现

【一】设计目的

1.熟悉IIR数字滤波器和FIR数字滤波器的设计原理和方法；

2.学会调用MATLAB信号处理工具箱中的滤波器设计函数设计各种IIR和FIR数字滤

波器，学会根据滤波要求确定滤波器指标参数；

3.掌握用IIR和FIR数字滤波器的MATLAB实现方法，并能绘制滤波器的幅频特性、

相频特性；

4.通过观察滤波器的输入、输出信号的时域波形及其频谱，建立数字滤波的概念。【二】设计内容

1.调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，观

察st的时域波形和幅频特性曲线；

2.通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号

的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率，假定要求滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB；

3.编程调用MATLAB滤波器设计函数分别设计这三个数字滤波器，并绘图显示其幅频特

性曲线；

4.用所设计的三个滤波器分别对复合信号st进行滤波，分离出st中的三路不同载波

频率的调幅信号，并绘图显示滤波后信号的时域波形和频谱，观察分离效果。【三】设计过程，调试，分析。

（一）IIR数字滤波设计方法

数字滤波是数字信号处理的一种重要算法，广泛用于对信号的过滤、检测与参数的估计等信号处理中。数字滤波器按其单位脉冲响应的长度可分为有限脉冲响应（FIR）滤波器和无限脉冲响应（IIR）滤波器两类。

数字滤波器按其通频带分为低通滤波器，高通滤波器，带通滤波器和带阻滤波器。

在本次实验中旨在用椭圆滤波器实现IIR低通，带通，高通滤波。

IIR 滤波器基本概念

一个 N 阶递归型数字滤波器的差分方程

H ( Z) =

ME br z- r

r = 0

N

1+K =E akz- k

从以上的系统函数可知, 设计 IIR 滤波器的任务就是通过计算寻求一个因果、物理上可实现的系统函数 H ( Z) , 使其频率响应 H ( e jw

) 满足所希望得到的频域指标, 即符合给定的通带截止频率、阻带截止、通带衰减和阻带衰减。 数字滤波器的设计过程

数字滤波器的一般设计过程可分为以下四个步骤: 1) 按照实际的需要, 确定滤波器的性能要求或指标;

2) 用一个因果稳定的离散线性时不变系统,去逼近这一性能指标; 3) 用有限精度的运算实现所设计的系统; 4) 通过模拟, 验证所设计的系统是否符合性能要求 （二）

1：调用信号产生函数mstg 产生由三路抑制载波调幅信号相加构成的复合信号st,画出st 的时域波形和幅频特性，三路信号在时域混叠无法在时域分离，但在频域分离，所以可以通过滤波的方法在频域分离。

抑制载波单频调幅信号的数学表达式为 []))(2cos())(2cos(2

1

)2cos()2cos()(000t f f t f f t f t f t s c c c ++-=

=ππππ 这三路调幅信号的载波频率分别为250Hz 、500Hz 和1000Hz

2：通过观察st 的幅频特性，分别可以确定三路信号的低通，带通，高通滤波器。的通带截止频率和阻带截止频率，要求滤波器的通带最大衰减为0.1db,阻带最小衰减为60db 分离混合信号st 的滤波器指标参数选取如下： 载波为250Hz 的调幅信号，用低通滤波器分离。指标为：

通带截止频率fp=300HZ,阻带截止频率fs=400Hz,通带最大衰减ap=0.1db,阻带最小衰减as=60db

对载波频率为500Hz 的调幅信号，用带通滤波器分离，指标为：

通带截止频率fl=400HZ,fu=500Hz,阻带截止频率fl=600Hz,fu=700Hz,通带最大衰减ap=0.1db,阻带最小衰减as=60db

对载波频率为1000Hz 的调幅信号，用高通滤波器分离，其指标为：

通带截止频率fp=900HZ,阻带截止频率fs=800Hz,通带最大衰减ap=0.1db,阻带最小衰减，60db

3：用椭圆滤波器实现 椭圆滤波器特点：

幅值响应在通带和阻带内都是等波纹的。

对于给定的阶数和给定的波纹要求，椭圆滤波器能获得较其它滤波器更窄的过渡带宽，就这点而言，椭圆滤波器是最优的。

通带和阻带内波纹固定时，阶数越高，过渡带越窄 阶数固定，通带和阻带波纹越小，过渡带越宽

【四】 程序清单与运行结果

程序清单：

%1. 调用信号产生函数mstg 产生由三路抑制载波调幅信号相加构成的复合信号st ，观察%st 的时域

% 波形和幅频特性曲线； function st=mstg

%产生信号序列st ，并显示st 的时域波形和频谱

%st=mstg返回三路调幅信号相加形成的混合信号，长度N=800

N=800; %信号长度N=800

Fs=10000; %采样频率Fs=10kHz,Tp为采样时间

T=1/Fs;

Tp=N*T;

t=0:T:(N-1)*T;

k=0:N-1;

f=k/Tp;

fc1=Fs/10; %第一路调幅信号载波频率fc1=1000Hz

fm1=fc1/10; %第一路调幅信号的调制信号频率fm1=100Hz

fc2=Fs/20; %第二路调幅信号载波频率fc2=500Hz

fm2=fc2/10; %第二路调幅信号的调制信号频率fm2=50Hz

fc3=Fs/40; %第三路调幅信号载波频率fc3=250Hz

fm3=fc3/10; %第三路调幅信号的调制信号频率fm3=25Hz

xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第一路调幅信号

xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第二路调幅信号

xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第三路调幅信号

st=xt1+xt2+xt3; %三路信号相加，得到复合信号

fxt=fft(st,N); %计算信号st的频谱

%以下为绘图命令

figure(1)

subplot(2,1,1)

plot(t,st);

grid;

xlabel('t/s');

ylabel('s(t)');

axis([0,Tp,min(st),max(st)]);

title('(a)s(t)的波形')

subplot(2,1,2)

stem(f,abs(fxt)/max(abs(fxt)),'.');

grid;

title('(b)s(t)的频谱')

axis([0,Fs/8,0,1.2]);

xlabel('f/Hz');

ylabel('幅度')

%低通滤波器

fp=320;fs=350;

wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B 和A

y1t=filter(B,A,st); %滤波器软件实现

y1=fft(y1t);

% 绘图部分

figure(2);

subplot(3,3,1);

[H,W]=freqz(B,A,800);

plot(W*Fs/2/pi,abs(H));

xlabel('Hz');

ylabel('H(w)');

title('低通滤波器');

axis([0,2000,0,1.2]);

grid;

subplot(3,3,2);

plot(t,y1t);

xlabel('t');

ylabel('y(t)');

title('分离出的250Hz的波形');

subplot(3,3,3);

stem(f,abs(y1)/max(abs(y1)),'.');

xlabel('Hz');

ylabel('|H|');

title('频谱');

axis([0,1200,0,1]);

grid;

%带通滤波器

fpl=420;fpu=580;fsl=300;fsu=700;

wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;

[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B 和A

y2t=filter(B,A,st); %滤波器软件实现

y2=fft(y2t);

%绘图部分

figure(2);subplot(3,3,4);

[H,W]=freqz(B,A,800);

plot(W*Fs/2/pi,abs(H));

xlabel('Hz');

ylabel('H(w)');

title('带通滤波器');

axis([0,2000,0,1.2]);

grid;

subplot(3,3,5);

plot(t,y2t);

xlabel('t');

ylabel('y(t)');

title('分离出的500Hz的波形');

subplot(3,3,6);

stem(f,abs(y2)/max(abs(y2)),'.');

xlabel('Hz');

ylabel('|H|');

title('频谱');

axis([0,1200,0,1]);

grid;

%高通滤波器

fp=870;fs=780;

wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp,'high'); %调用ellip计算椭圆带通DF系统函数系数向量B和A

y3t=filter(B,A,st); %滤波器软件实现

y3=fft(y3t);

%绘图部分

figure(2);subplot(3,3,7);

[H,W]=freqz(B,A,800);

plot(W*Fs/2/pi,abs(H));

xlabel('Hz');

ylabel('H(w)');

title('高通滤波器');

axis([0,2000,0,1.2]);

grid;

subplot(3,3,8);

plot(t,y3t);

xlabel('t');

ylabel('y(t)');

title('分离出的1000Hz的波形');

subplot(339);

stem(f,abs(y3)/max(abs(y3)),'.');

xlabel('Hz');

ylabel('|H|');

title('频谱');

axis([0,1200,0,1]);

grid;

运行结果：

设计四语音信号滤波处理

【一】设计目的

1.了解语音信号的产生、采集，能绘制语音信号的频率响应曲线及频谱图；

2.学会用MATLAB对语音信号进行分析和处理；

3.掌握用滤波器去除语音信号噪声的方法，观察去噪前后的语音信号。

【二】设计内容

1.利用Windows下的录音机录制一段自己的话音，时间在1s内。然后在Matlab软件

平台下，利用wavread函数对语音信号进行采样，记住采样频率和采样点数；

2.画出语音信号的时域波形，对采样后的语音信号进行快速傅立叶变换，得到信号的

频谱特性，画出采样信号的时域波形和频谱图；

3.根据对语音信号谱分析的结果，确定滤除噪声所需滤波器的技术指标，设计合适的

数字滤波器，并画出滤波器的频域响应；

4. 用所设计的滤波器对采集的信号进行滤波，在同一个窗口画出滤波前后信号的时域

波形和频谱，并对滤波前后的信号进行对比，分析信号的变化 ； 5. 回放语音信号；sound(x)

6. 为使编制的程序操作方便，设计一个信号处理用户界面。在所设计的系统界面上可

以选择滤波器的类型，输入滤波器的参数，显示滤波器的频率响应，选择信号等。

【三】设计原理 1. 语音信号的采集

在MATLAB 软件平台下，利用函数wavread （ ）对语音信号采集，并记录采样频率和采样点数。将语音信号转换成计算机能够运算的有限长序列。

wavread 函数的调用格式如下：

y=wavread(file)

读取file 所规定的wav 文件，返回采样值放在向量y 中。

[y,fs]=wavread(file)

采样值放在向量y 中，fs 表示采样频率（Hz ）。 2. 用FFT 作谱分析

FFT 即快速傅立叶变换，它是从DFT 运算中发展起来的，利用系数nk

N W 的对称性和周期性减少运算量。长度为N 的序列直接计算DFT 需要2

N 次复乘和(1)N N -次复加，而用FFT 进行运算一般需要

2log 2

N

N 次复乘和2log N N 次复加，

从而使DFT 的运算大大简化。

用FFT 对连续信号进行谱分析的步骤如下：

图3.1 连续信号谱分析过程

引入前置低通滤波器LPF 是为了消除或减少时域连续信号转换成序列时可能出现的频谱混叠现象。)(n w 表示时域有限的窗函数。 3. 设计滤波器去除语音信号的噪声

通过wavread（）函数将语音信号读入，通过频率采样及fft ( )产生信号，并对之加噪，通过窗函数法设计滤波器滤掉该语音信号的噪声，并对比滤波前后的语音波形和频谱。

【四】程序清单与运行结果：

程序清单：

Fs =20000; %采样频率

[x,Fs]=wavread('E:\sound.wav');

figure(1)

subplot(2,2,1)

plot(x,'b'); %做原始语音信号的时域图形

title('原语音信号');grid on;

%sound(x,Fs);

subplot(2,2,2)

plot(abs(fft(x)));%画出原始语音信号的频谱图

title('原语音信号的频谱');grid on;

y=rand(size(x))/5;

y=x+y;

subplot(2,2,3);

plot(y); %做噪声语音信号的时域图形

title('加噪后语音信号');

xlabel('时间n');ylabel('音量n');

%sound(y,Fs);

subplot(2,2,4);

y5=fft(y,length(y)); %做length(y)点的FFT

plot(abs(y5));%画出噪声语音信号的频谱图

title('加噪后语音信号的频谱');grid on;

fs=20000;

fc2=5000;

wc2=2*pi*fc2/fs;

wp2=2*pi*4800/fs;ws2=2*pi*5000/fs;

N2=ceil(6.6*pi/(ws2-wp2))+1;%确定窗的长度

Window=(hamming(N2));

b2=fir1(N2-1,1/pi,0,Window);%利用窗函数法设计FIR滤波器

[H,W]=freqz(b2,1,512);%数字滤波器频率响应

H_db=20*log10(abs(H));

figure(2);

subplot(3,1,1)

plot(W,H_db);

title('低通滤波器');grid on;

y_high = filter(b2,1,y);%对信号进行高通滤波

subplot(3,1,2)

plot(y_high);title('信号经过低通滤波器(时域)'); subplot(3,1,3)

plot(abs(fft(y_high)));

title('信号经过低通滤波器(频域)');grid on; sound(y_high,Fs);

运行结果：