数字信号处理上机实习报告.doc

数字信号处理上机实习报告

———————————————————————————————— 作者：
———————————————————————————————— 日期：

28

个人收集整理 勿做商业用途
专题一 离散卷积的计算

一、实验内容
设线性时不变（LTI）系统的冲激响应为h(n),输入序列为x（n)
1、h（n)=(0.8)n，0≤n≤4； x（n)=u（n)—u（n-4)
2、h（n)=（0.8）nu（n)， x（n）=u(n)—u（n-4）
3、h(n)=(0。8)nu（n）, x(n)=u（n)
求以上三种情况下系统的输出y（n）。
二、实验目的
1、掌握离散卷积计算机实现.
2、进一步对离散信号卷积算法的理解.
三、原理及算法概要
算法：把冲激响应h(n)与输入序列x（n)分别输入到程序中,然后调用离散卷积函数y=conv（x。，h)即可得到所要求的结果.
原理:
离散卷积定义为

当序列为有限长时，则

四．理论计算
1、h（n)=(0。8)n，0≤n≤4; x（n）=u（n）—u(n—4)

（a) 当n〈0 时，y(n）=0
(b） 当时,（0。8）n
(c） 当时，（0.8)n
(d） 当n〉7时，y（n）=0
理论结果与上图实验结果图中所示吻合。

2、h(n)=(0.8)nu(n）, x(n）=u(n)—u(n-4)

(a) 当n<0 时,y(n）=0
(b) 当时，(0。8)n
（c） 当时，（0。8)n
（d) 当时，（0.8）n
（e) 当n>23时，y（n）=0
在卷积序列中,因为n趋于无穷，在实验中无法实现，而0。8的20次幂基本接近于0，所以这里只取到n为20就可以了。由以上可以看出，理论结果与上图实验结果图中所示吻合。

3、h(n）=（0.8)nu（n）， x(n)=u(n）

(a) 当n〈0 时,y（n)=0
(b) 当时，(0。8）n
(c) 当时，(0.8)n
(d) 当n〉140时,y（n）=0
在卷积序列中，因为n趋于无穷，在实验中无法实现，而0。8的70次幂的叠加基本接近于最大值,所以这里只取到n为70就可以了。由以上可以看出，理论结果与上图实验结果图中所示吻合。

五．程序
x1=[1 1 1 1 ］；nx1=0：3；
h1=[1 0。8 0.64 0。8^3 0。8^4］;nh1=0:4;
y1=conv（x1,h1)；
subplot（3，3,1)；stem(nx1，x1）;title（'序列x1');
xlabel（'n'）；ylabel（’x1(n)'）;
subplot(3,3,2)；stem(nh1，h1）；title（’序列h1’)；
xlabel('n'）；ylabel（’h1(n)’)；
subplot（3,3,3)；stem(y1）;title（'序列y1');
xlabel（'n'）；ylabel（'y1（n）’)；

x2=[1 1 1 1]；nx2=0：3；
nh2=0：1:20；
h2=（0。8)。^nh2;
y2=conv（x2,h2)；
subplot（3，3，4）;stem(x2)；title（’序列x2'）;
xlabel（’n’);ylabel(’x2（n）');
subplot(3,3，5);stem(h2);title（'序列h2')；
xlabel（'n’);ylabel('h2(n）’）；
subplot（3，3,6);stem(y2)；title（'序列y2');
xlabel（’n'）；ylabel（’y2(n)')

nx3=0:1：20；
x3=1.^nx3;
nh3=0:1：20;
h3=(0.8).^nh3；
y3=conv（x3,h3)；
subplot（3,3，7);stem（nx3，x3）；title(’序列x3')；
xlabel('n');ylabel('x3（n)’）;
subplot（3,3,8)；stem(nh3，h3）；title(’序列h3')；
xlabel('n'）；ylabel（’h3（n）’）；
subplot(3，3，9）；stem(y3)；title('序列y3’）；
xlabel（’n'）;ylabel（'y3(n）'本文为互联网收集,请勿用作商业用途个人收集整理，勿做商业用途
六．程序运行结果

七．结果分析
有限长序列的离散卷积计算结果与理论值一致，而存在无限长序列做卷积时,由于在程序处理时是用比较长有限长序列代替的，所以与理论值基本相同.
八.专题实习心得
该专题主要进行的是有matlab实现离散卷积的计算，分为三个类型：冲激响应h(n）与输入序列x(n）都为有限长，一个序列为有限长一个序列为无限长和两个序列均为无限长。这部分实习内容不难,是原来做过的，主要是为了拾拣起原来所学的知识。第一种类型是最简单,也是最基本的,直接调用函数y=conv(x.,h)即可。在第二种和第三种类型的计算中遇到了一些困难，在输入序列时，由于存在无限长的序列，不知道该怎么输入，查过了相关的题目才弄明白。通过本专题的实习，让我对数字信号处理中的一些基础知识有了一些回顾，对原来所学过的知识也熟悉了一些。

专题二 离散傅里叶变换及其应用

一、 实验内容
设有离散序列 x(n）=cos（0。48πn)+cos(0.52πn）
分析下列三种情况下的幅频特性。
（1） 采集数据长度N=16，分析16点的频谱,并画出幅频特性。
（2） 采集数据长度N=16，并补零到64点，分析其频谱,并画出幅频特性。
（3） 采集数据长度N=64，分析46点的频谱,并画出幅频特性。
观察三幅不同的幅频特性图，分析和比较它们的特点及形成原因。
2、实现序列的内插和抽取所对应的离散傅里叶变换。（选做）

求和对应的128点傅里叶变换，比较三个结果所得到的幅频特性,分析其差别产生的原因。
二、实验目的
1、了解DFT及FFT的性质和特点
2、利用FFT算法计算信号的频谱。
三、关键算法
算法1：读入离散序列x（n)=cos（0.48πn)+cos（0.52πn)，采集长度为N=16的数据，调用matlab中的函数fft（x,16）与fft(x,64)对其作离散傅里叶变换得到16点、64点的频谱。采集数据长度为N=64,，调用matlab中的函数fft(x，46）对其作离散傅里叶变换得到46点的频谱。
算法2：读入离散序列信号，每隔N=4采集一个数值，得到新的抽取离散序列，每隔K=0.25采集一个数值得到新的内插离散序列。分别对其做离散傅氏变换,得到对应的频谱图。

原理概要：
当抽样数N=2M时,以下为算法蝶形图.

一般规律如下:
1、 当N=2M时，则要进行M次分解，即进行M级蝶形单元的计算
2、按自然顺序输入，输出是码位倒置。
3、每一级包含N/2个基本蝶形运算
4、第L级有2L—1个蝶群，蝶群间隔为N/2L-1
如果是Matlab实现的话，可用以下两种方法计算信号频谱
1、调用库函数为:fft(），直接计算X（k）
2、进行矩阵运算
四．程序
程序1：
n=0：1：15；
x1=cos（0。48＊3。14＊n）+cos(0.52＊3。14＊n）;
g1=abs(fft（x1，16））；
subplot（3,2，1)；stem(x1）;title（'x1')；
subplot（3，2，2)；stem（g1);title（'g1’）；

n2=0：1:15；
x2=cos(0.48*3.14＊n2）+cos(0。52＊3。14＊n2）;
x2=［x2 zeros(1,48）］;
g2=abs(fft（x2,64）);
subplot（3，2，3)；stem(x2)；title（'x2’）;
subplot（3，2,4）；stem（g2）；title（'g2')；

n3=0：1：64；
x3=cos（0。48＊3.14＊n3)+cos(0.52＊3.14＊n3）;
g3=abs（fft(x3,46)）;
subplot(3,2，5）；stem（x3）；title(’x3’）;
subplot(3，2，6）；stem（g3)；title(’g3'）;

程序2：
x=［];
for n=0:127
x=［x (cos（1/36＊pi*n)+cos（1。5/36＊pi*n））］
end
X=abs(fft（x)）；
subplot（3,2,1）;stem（x）;title（’序列x'）；
xlabel（'n’)；ylabel(’x（n）')；
subplot（3,2，2);stem(X)；title（’序列X');
xlabel(’\omega’);ylabel('X')；

x1=［];
for n=0：4：508
x1=［x1 （cos（1/36*pi＊n）+cos(1。5/36*pi＊n））]
end
X1=abs（fft（x1)）；
subplot(3，2，3)；stem(x1）;title（'序列x1');
xlabel('n’);ylabel('x1(n)’）；
subplot(3，2,4);stem（X1）；title(’序列X1'）；
xlabel（'\omega')；ylabel(’X1’）;

x2=[］；
for n=0:0.25：31。75
x2=［x2 (cos（1/36*pi*n）+cos（1.5/36＊pi＊n)）]
end
X2=abs(fft（x2）);
subplot（3，2，5）;stem(x2）；title(’序列x2');
xlabel（’n');ylabel（'x2（n）’)；
subplot（3,2,6);stem（X2)；title('序列X2’）;
xlabel（’\omega'）;ylabel（’X2’)；
本文为互联网收集，请勿用作商业用途本文为互联网收集,请勿用作商业用途
五．程序运行结果
结果1：

结果2:

六．结果分析
　N 点DFT的频谱分辨率是2 π/N.一节指出可以通过补零观察到更多的频点,但是这并不意味着补零能够提高真正的频谱分辨率。这是因为x［n] 实际上是x(t） 采样的主值序列，而将x[n]补零得到的x'［n] 周期延拓之后与原来的序列并不相同，也不是x（t） 的采样。因此是不同离散信号的频谱。对于补零至M点的x’的DFT，只能说它的分辨率2 π/M仅具有计算上的意义,并不是真正的、物理意义上的频谱。频谱分辨率的提高只能通过提高采样频率实现.
七。专题实习心得
离散傅里叶变换是一种快速算法,由于有限长序列在其频域也可离散化为有限长序列，因此离散傅里叶变换在数字信号处理中是非常有用的.DFT是重要的变换，在分析有限长序列的有用工具、信号处理的理论上有重要意义、运算方法上起核心作用，谱分析、卷积、相关都可以通DFT在计算机上实现。DFT解决了两个问题：一是离散与量化，二是快速运算.通过编程实践体会到了时域、频域信号的对应关系,也对采样频率的含义有了深刻的认识,同时也加深了对采样信号频谱周期性的理解.

八．思考题
1、直接计算DFT与用FFT计算,理论上的速度差别应是多少?
整个DFT运算总共需要4N 2次实数乘法和2N(2N—1）次实数加法。整个FFT运算总共需要（N2/2)+（N/2）=N（N+1)/2≈N2/2次复数乘法和N（N/2-1)+N=N2/2次复数加法。由此可见，理论上用FFT计算运算速度是DFT的二倍。
2、 什么是高密度频谱及高分辨率频谱？实验1中至少应采集几个样本才能区分两个频率分量？
DFT的频谱分辨率一定时,在尾部补零可以得到DFT的高密度频谱，可以细化当前分辨率下的频谱，但不能改变DFT的分辨率。高分辨率是指对信号中两个靠得较近的频谱分量的识别能力比较高，在采样频率不变时,采样点数N越大，频谱分辨率越高。在实验1中至少要采集8个点才能区分两个频率分量。

专题三 IIR滤波器的设计
一、实验内容
1、设计一个Butterworth数字低通滤波器，设计指标如下:
通带截止频率：0.2π，幅度衰减不大于1分贝
阻带截止频率：0.3π,幅度衰减大于15分贝
2、让不同频率的正弦波通过滤波器，验证滤波器性能。
3、分析不同滤波器的特点和结果。
4、编程设计实现IIR滤波器.
二、实验目的
掌握不同IIR滤波器的性质、特点。
通过实验学习如何设计各种常用的IIR滤波器，以便在实际工作中能根据具体情况使用IIR滤波器。
三、原理及算法概要
算法：输入通带截止频率Wp，阻带截止频率Ws,通带衰减Rp，阻带衰减Rs，通过这些数值调用［N Wn］=buttord（Wp,Ws,Rp，Rs） 函数计算巴特沃斯数字滤波器的阶数N和截止频率wn，再根据阶数N通过函数［b，a］=butter(N,Wn），即可得到所要的巴特沃斯滤波器。设计一个正弦波信号，再调用函数A=filter（b，a，I）让正弦波信号通过滤波器,得到滤波信号.
原理概要：
1、滤波器类型
Butterworth滤波器
Butterworth滤波器的特点是在通带内的频率特性是平坦的，并且随着频率的增加而衰减.Butterworth滤波器又是最简单的滤波器。
N阶低通Butterworth滤波器的幅度平方函数为:
Chebyshev滤波器
Chebyshev滤波器的在通带内的响应是等纹波的，而在阻带内是单调下降的，或在通带内是单调下降的,在阻带内是等纹波的特性
其幅度平方函数为
其中VN是Chebyshev多项式.
椭圆滤波器
椭圆滤波器在通带和阻带内都是等纹波振荡。
椭圆滤波器的特性函数为:
其中UN是N阶雅可比函数。

模拟滤波器设计完成后，就可以进行典型滤波器设计的第二个步骤,通过冲激响应不变法和双线性变换转变为相应的数字滤波器.
2、变换方法
（1）冲激响应不变法
冲激响应不变法的基本原理是从滤波器的冲激响应出发，对模拟滤波器冲激响应h（t）进行取样，所得到的离散序列h（nT)作为数字滤波器的单位取样响应。
H(z）是由H（s）通过下式的对应关系得到.

（2）双线性变换是在所得到满足性能指标要求的模拟滤波器的基础上,通过变换
,从而得到相应的数字滤波器。
四．程序
Wp=0.2;
Ws=0。3;
Rp=1；
Rs=15；
［N Wn]=buttord（Wp，Ws，Rp,Rs) %用于计算巴特沃斯数字滤波器的阶数N和截止频率wn

［b，a]=butter（N,Wn)； %计算N阶巴特沃斯数字滤波器系统函数分子、分母多项式的系数向量b、a，设计所需的低通滤波器

［h，omega]=freqz(b，a,512);％返回量h包含了离散系统频响 ，调用中若N默认,默认值为512。

plot（omega/pi，20＊log10(abs（h）));grid;
xlabel（’\omega/\pi’);ylabel（’Gain,dB’);
title(’IIR Butterworth Lowpass Filter’）；

Wp=0.2；
Ws=0。3；
Rp=1;
Rs=15;
[N1，Wn1]=buttord（Wp，Ws，Rp,Rs）；％用于确定阶次
[b,a］=butter（N,Wn）；%用于直接设计巴特沃兹数字滤波器，即为IIR滤波器
%freqz（b,a)；
t=1:300
I=sin（0。1＊pi*t)+sin(0。4*pi＊t）;%设计正弦波
plot(I)；
figure；
A=filter(b,a,I）;%正弦波通过滤波器
plot(A）；

本文为互联网收集，请勿用作商业用途本文为互联网收集，请勿用作商业用途
五．程序运行结果
N =

6

Wn =

0.2329

N1 =
6
Wn1 =

0。2329

六．结果分析
Butterworth滤波器在通带内的频率特性是平坦的，并且随着频率的增加而衰减。正弦信号在经过IIR滤波器滤波后,高频信号被滤除，低频信号被保留了下来。
七。专题实习心得
通过本专题的学习，熟悉和巩固了Butterworth数字低通滤波器的设计方法和原理，实现滤波器设计的有关经典算法，更重要的是熟练掌握使用MATLAB语言设计各种要求的数字滤波器。这一部分的内容相对来说是有些难度了,做起来花费的精力也多了一些,不过对数字信号处理内容的掌握上又加深了一层。
八．思考题
1、比较双线性变换与冲激响应不变法的优缺点.
冲激响应不变法的特点：
模仿模拟滤波器的单位抽样相应 ,时域模仿性好产生频率响应的混叠失真。保持系统的相位特性换不改变系统的因果性和稳定性；
因为存在频率响应的混叠效应，只适用于限带的模拟滤波器，高通和带阻滤波器不宜采用冲激响应不变法,冲激响应不变法是一种时域模仿，缺点是产生频率响应的混叠失真。冲激响应不变法产生混叠失真的原因是映射是多对一的映射。双线性变换的特点是一对一的映射,消除了多值变换性，也就消除了频谱混叠.
1、 和频率变换相比，数字变换有什么优点？
数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变

3、实验中求取相应模拟滤波器原形时,为什么要对模拟频率Ωc进行归一化？
用归一化巴特沃思低通滤波器为原型滤波器，一旦归一化低通滤波器的系统函数确定后，其它巴特沃思低通滤波、高通、带通、带阻滤波器的传递函数都可以通过变换法从归一化低通原型的传递函数 Ha(S) 得到。归一化原型滤波器是指截止频率Ωc 已经归一化成Ωc’=1的低通滤波器。对于截止频率为某个Ωc的低通滤波器，则令 S/Ωc 代替归一化原型滤波器系统函数中的 S∧ ，即 S∧ S/Ωc 对于其他高通、带通、带阻滤波器，可应用后面讨论到的频带变换法，由其变换得出.

专题四 用窗函数设计FIR 滤波器
一、实验内容
选取合适窗函数设计一个线性相位FIR低通滤波器，使它满足如下性能指标：
通带截止频率：ωp=0.5π,通带截止频率处的衰减不大于3分贝；
阻带截止频率：ωs=0。66π,阻带衰减不小于40分贝.
二、实验目的
1、掌握用窗函数法设计FIR滤波器的原理和方法。
2、熟悉线性相位滤波器特性。
3、了解各种窗函数对滤波器特性的影响。
三、原理及算法概要
算法：通过其通带截止频率ωp与阻带截止频率ωs算出其过渡带的宽度与滤波器的长度，从而得到理想滤波器的截止频率,根据所要求的理想滤波器,得到hd(n）。由于其通带截止频率处的衰减不大于3分贝与阻带衰减不小于40分贝，我选择最接近的汉宁窗，最后调用函数h=hd.*win 及freqz(h，1，512)得到实际汉宁窗的响应和实际滤波器的幅度响应。

原理概要：
1、利用窗函数法设计FIR 滤波器
FIR滤波器的最大特点是其相位特性可以设计为严格的线性,而其幅值可以任意设置，这样输出波形就不会相位失真.
理想低通滤波器的单位取样响应hd（n)是无限长的，所以要用一个有限长的因果序列
h(n)进行逼近，最有效的方法是截断hd（n)，即用有限长的窗函数w（n）来截取
hd（n），表示为h（n)=hd（n）w（n）
为获得线性相位的FIR滤波器，h（n)必须满足中心对称条件，序列h（n)应有一定的延迟α，且α=（N—1）/2
频率响应逼近hd（ejw）的FIR滤波器，最简单的窗函数为矩形窗
1 n<（N—1）/2
W(n）=
0 n〉（N—1)/2
加窗后的频谱 加窗后使实际频响偏离理想频响，影响主要有两个方面:
（1） 通带和阻带之间存在过渡带，过渡带宽度取决于窗函数频响的主瓣宽度。
（2） 通带和阻带区间有纹波，这是由窗函数的旁瓣引起的，旁瓣越多，纹波越多.
增加窗函数的宽度N，其主瓣宽度减小，但不改变旁瓣的相对值。
为了改善滤波器的性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；
旁瓣衰减尽可能大，数量尽可能大,从而改善纹波状况，使实际频响H（ejω）更好地逼近理想频响Hd（ejω）

除了矩形窗外，一般还可以采用以下几种窗函数
Battlet窗 ：
汉宁窗：
海明窗
布来克曼窗
其中:WR是矩形窗的频谱

四．程序
wp =0。5*pi；
ws=0。66＊pi；
wdelta =ws-wp; %过渡带宽度
N=ceil（8*pi/wdelta) %滤波器长度
if rem(N,2)==0
N=N+1;
end
Nw =N；
wc =（wp+ws）/2； %理想低通滤波器的截止频率
n =0： N—1；
alpha =(N—1)/2;
m =n—alpha+0。00001；
hd =sin（wc*m）。/(pi*m); ％一个响应
win =(hanning（Nw）)’； ％汉宁窗
h=hd。*win; %实际汉宁窗的响应
freqz（h，1，512); ％实际滤波器的幅度响应

五．程序运行结果
N =

50

六．结果分析
上图中第一个图为FIR 滤波器的相频特性，第二个图为幅频特性；根据滤波器特性：通带截止频率处的衰减不大于3分贝；阻带衰减不小于40分贝。选用汉宁窗.FIR滤波器的最大特点是其相位特性可以设计为严格的线性，而其幅值可以任意设置,这样输出波形就不会相位失真。理想低通滤波器的单位取样响应hd（n）是无限长的,所以要用一个有限长的因果序列
h（n）进行逼近,最有效的方法是截断hd（n），即用有限长的窗函数w(n）来截取
hd（n)，表示为h(n）=hd（n）w(n)
为获得线性相位的FIR滤波器，h（n）必须满足中心对称条件,序列h（n）应有一定的延迟α,且α=（N-1)/2
频率响应逼近hd(ejw）的FIR滤波器，最简单的窗函数为矩形窗
1 n<(N-1)/2
W（n）=
0 n>（N—1）/2
加窗后的频谱 加窗后使实际频响偏离理想频响,影响主要有两个方面：
（3） 通带和阻带之间存在过渡带，过渡带宽度取决于窗函数频响的主瓣宽度。
（4） 通带和阻带区间有纹波,这是由窗函数的旁瓣引起的，旁瓣越多,纹波越多.
增加窗函数的宽度N，其主瓣宽度减小，但不改变旁瓣的相对值。
为了改善滤波器的性能,要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带;
旁瓣衰减尽可能大，数量尽可能大，从而改善纹波状况，使实际频响H（ejω)更好地逼近理想频响Hd（ejω）

七。专题实习心得
FIR有限冲击响应数字滤波器由于设计灵活，在滤波效果和响应时间之间能较好的折衷,因此在数字信号处理领域应用广泛。用窗函数法设计FIR数字滤波器在截断点依旧取模拟样本冲激响应的全值更利于理解和仿真.但是，窗函数法设计FIR数字滤波器的频率特性与模拟样本频率特性不完全一样,且具有线性相移的特点.窗函数法设计FIR数字滤波器是傅立叶变换的典型运用.通过学习有益于加深对数字滤波器的理解，并提高使用傅立叶变换分析解决实际问题的能力。

八:思考题
1、 在这个实验中，你选择哪一种窗函数？试说明原因.
选用汉宁窗。根据所给定滤波器特性：通带截止频率处的衰减不大于3分贝;阻带衰减不小于40分贝。参照书本上的各滤波器的特性。
2、用窗函数法设计FIR 滤波器时，滤波器的过渡带宽度和阻带衰减和哪些因素有关？
窗函数的主瓣宽度决定滤波器过渡带宽度，旁瓣幅值决定阻带衰减。选择窗函数应使其频谱满足：主瓣宽度尽可能的窄，以使过渡带尽可能陡；最大旁瓣相对于主瓣尽可能小，使能量尽可能集中在主瓣内，这样可以使肩峰和波动减小.

综合
一、实验内容
录制一段自己的语音信号，时间为10s左右，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法或双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的语音信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比,分析信号的变化；回放语音信号；最后，用MATLAB设计一信号处理系统界面。
二．算法
调用函数function pushbutton1_Callback(hObject, eventdata, handles）实现一个信号处理系统界面。选择左键时，用双线性变换法设计滤波器来对信号进行处理，选择右键时,用窗函数法设计滤波器来对信号进行处理。读取语音信号,对语音信号进行f=8000的频率进行采样,调用函数y1=fft(x1，2048）对所采集的点做2048点FFT变换。先设计butterworth模拟滤波器,再用双线性变换法实现模拟滤波器到数字滤波器的转换.最后调用函数f1=filter（bz，az,x2)对加了噪声的语音信号进行滤波，得到滤波后的频谱图。

三．程序
function pushbutton1_Callback（hObject, eventdata， handles)
% hObject handle to pushbutton1 （see GCBO）
% eventdata reserved — to be defined in a future version of MATLAB
％ handles structure with handles and user data （see GUIDATA)
fs=8000； ％语音信号采样频率为8000
x1=wavread（'F:\实习\信息工程实习\第五题\语音.wav');
t=（0：length（x1)—1）/8000;
y1=fft(x1，2048）; %对信号做2048点FFT变换
f=fs*（0：1023)/2048;
figure(1）
subplot(2，2,1）
plot（t，x1) ％做原始信号的时域波形
grid on；axis tight；
title('原始语音信号');
xlabel(’time(s)')；
ylabel(’幅度');
subplot(2，2,2)
plot（f，abs（y1(1：1024))) ％做原始信号的FFT频谱
grid on;axis tight;
title（’原始语音信号的FFT频谱'）
xlabel('Hz’)；
ylabel（'幅度’);
%双线性变换法设计的巴特沃斯滤波器
A1=0。05;A2=0.10；
d=[A1＊cos(2*pi*3800*t）+A2＊sin（2*pi*3600＊t）]';
x2=x1+d;
wp=0.8＊pi;
ws=0.85*pi;
Rp=1;
Rs=15；
Fs=8000；
Ts=1/Fs;
wp1=2/Ts*tan(wp/2); ％将模拟指标转换为数字指标
ws1=2/Ts*tan(ws/2）;
[N，Wn］=buttord（wp1,ws1,Rp,Rs,’s’）； %选择滤波器最小阶数
［Z，P,K］=buttap（N); %创建butterworth模拟滤波器
[Bap,Aap］=zp2tf（Z,P,K）;
[b，a]=lp2lp(Bap，Aap,Wn）；
[bz，az］=bilinear(b,a，Fs）; %用双线性法实现模拟到数字的转换
［H，W］=freqz（bz,az)； ％绘制频率响应曲线
subplot（2，2,3)
plot（W＊Fs/(2*pi),abs（H)）
grid on;axis tight;
xlabel(’频率（Hz）’)
ylabel（’频率响应’）
title('Butterworth'）
f1=filter（bz，az,x2)；
figure（2)
subplot(2，2,1）
plot（t，x2） ％画出滤波前的时域图
grid on;axis tight;
title('滤波前的时域波形’）；
subplot（2,2，2)
plot(t，f1）; ％画出滤波后的时域图
grid on；axis tight;
title(’滤波后的时域波形’）;
y3=fft(f1，2048);
y2=fft（x2，2048);
subplot(2，2,3）；
plot(f,abs（y2（1:1024)））; %画出滤波前的频谱图
grid on;axis tight;
title（’滤波前的频谱’)
xlabel('Hz')；
ylabel（'幅度'）；
subplot（2，2，4)
plot(f，abs(y3(1:1024))）； %画出滤波后的频谱图
grid on;axis tight；
title（'滤波后的频谱')
xlabel('Hz'）；
ylabel（'幅度'）；

function pushbutton2_Callback(hObject， eventdata, handles）
% hObject handle to pushbutton2 （see GCBO)
% eventdata reserved — to be defined in a future version of MATLAB
% handles structure with handles and user data （see GUIDATA）
fs=8000； ％语音信号采样频率为8000
x1=wavread(’F：\实习\信息工程实习\第五题\语音。wav’）;
t=(0：length（x1)-1)/8000;
y1=fft（x1,2048)； %对信号做2048点FFT变换
f=fs＊(0：1023）/2048；
figure(1)
subplot（2,1,1)
plot（t,x1) ％做原始信号的时域波形
grid on;axis tight;
title(’原始语音信号’）;
xlabel（'time(s)’）；
ylabel（'幅度’);
subplot（2,1，2)
plot(f，abs（y1（1:1024)）） %做原始信号的FFT频谱
grid on;axis tight；
title（'原始语音信号的FFT频谱’）
xlabel（'Hz'）；
ylabel(’幅度'）;
%窗函数设计滤波器
t=（0:length(x1)-1)/8000；
f=fs＊(0：2047)/4096；
A1=0.05；A2=0。10；
d=[A1＊cos(2＊pi＊3600＊t)+A2*sin（2*pi＊3800＊t)]’；
x2=x1+d；
wp=0.8*pi;
ws=0.85＊pi；
wdelta=ws—wp;
N=ceil(6。6＊pi/wdelta）; %取整
wn=（0。8+0。85）＊pi/2；
［bz，az]=fir1(N,wn/pi，hamming(N+1))； ％选择窗函数,并归一化截止频率
figure(2)
freqz(bz,az);
grid on;axis tight；
f2=filter(bz，az,x2);
figure（3)
subplot(2,2，1）
plot(t,x2);
grid on;axis tight;
title（'滤波前的时域波形'）；
subplot(2，2，2)
plot(t,f2);
grid on；axis tight;
title（’滤波后的时域波形');
y3=fft(f2，4096);
f=fs＊（0:2047)/4096;
figure(8)
y2=fft（x2，4096);
subplot（2，2,3）；
plot(f，abs（y2(1：2048)）);
grid on；axis tight;
title（'滤波前的频谱’）
xlabel（’Hz’）；
ylabel（'幅度');
subplot(2,2，4)
plot（f，abs(y3（1:2048)））;
grid on;axis tight；
title(’滤波后的频谱’）
xlabel（’Hz’）；
ylabel(’幅度');

四．实验结果 个人收集整理，勿做商业用途文档为个人收集整理，来源于网络
实验结果：

选择左边的按键：

选择右边的按键:

汉明窗：

五．结果分析
因为录音是自己用电脑自己说话来录音的，所以从语音的频谱可以看出，几乎没有高频信号，之前所设计的滤波器都是低通滤波,为了再使用低通滤波器,所以人为加了高频信号，通过低通滤波器之后，滤波后的时域波形对滤波前的时域波形看起来相对光滑一些，滤波后的频谱图，滤除了添加的高频噪声信号，保留了低频信号，语音效果也改善了不少。通过最后的频谱图可以看出,与之前语音信号的频谱图几乎一样。而通过比较双线性变换法与窗函数设计的滤波器可以看出，它们的滤波效果相似。
六．专题实习体会
本部分的实习比较难,虽然读了不少相关的程序，查了很多资料，做起来感觉很不容易.这一阶段的实习让我体会到了能读懂程序和自己编写程序的差距是很大的,理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，经过大量的亲身实践后才能够得到更深的掌握,最终达到融会贯通的效果。
本次数字信号处理实习体会：
数据信号处理是一种很有用的课程，但课程本身理论性强、公式推导较多、概念比较抽象，通过上机的实践可加深对课本理论知识的理解和认识。通过学习可以设计出良好的算法,高效地组织处理数据.这门课在学习的时候做过类似的程序设计，想着应该不会太难，可是拿到题目了以后，才发现我们所掌握的知识和实际问题有些联系不上，才发现理论与实际相结合说起来容易,要做起来需要很多的练习实践才能达到。由于时间过得有点久了,很多知识都生疏了,但通过一段时间的学习后，把原来学的东西都拣起来了不少，看起程序来也