《数字信号处理》实验报告

《数字信号处理》
实验报告
年级：2011级
班级：信通4班
姓名：朱明贵
学号：111100443
老师：李娟
福州大学
2013 年11 月
实验一快速傅里叶变换(FFT)及其应用
一、实验目的
1.在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。

2.熟悉应用FFT对典型信号进行频谱分析的方法。

3.了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。

4.熟悉应用FFT实现两个序列的线性卷积和相关的方法。

二、实验类型
演示型
三、实验仪器
装有MATLAB语言的计算机
四、实验原理
在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换(DFT)。

这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT定义为：
反变换为：
有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。

FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。

它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。

常用的FFT 是以2为基数的，其长度。

它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。

(一)在运用DFT进行频谱分析的过程中可能的产生三种误差
1．混叠
序列的频谱是被采样信号频谱的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。

避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。

2．泄漏
实际中我们往往用截短的序列来近似很长的甚至是无限长的序列，这样可以使用较短的DFT 来对信号进行频谱分析，这种截短等价于给原信号序列乘以一个矩形窗函数，也相当于在频域将信号的频谱和矩形窗函数的频谱卷积，所得的频谱是原序列频谱的扩展。

泄漏不能与混叠完全分开，因为泄漏导致频谱的扩展，从而造成混叠。

为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减至最小。

3．栅栏效应
DFT 是对单位圆上Z 变换的均匀采样，所以它不可能将频谱视为一个连续函数，就一定意义上看，用DFT 来观察频谱就好像通过一个栅栏来观看一个图景一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住，不能别我们观察到。

减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值，从而变动DFT 的点数，这一方法实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或谷点暴露出来。

(二)用FFT 计算线性卷积
用FFT 可以实现两个序列的圆周卷积。

在一定的条件下，可以使圆周卷积等于线性卷积。

一般情况，设两个序列的长度分别为N 1和N 2，要使圆周卷积等于线性卷积的充要条件是FFT 的长度N ≥N 1＋N 2-1，对于长度不足N 的两个序列，分别将他们补零延长到N 。

当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。

有两种方法：
1． 重叠相加法
将长序列分成与短序列相仿的片段，分别用FFT 对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。

2．重叠保留法
这种方法在长序列分段时，段与段之间保留有互相重叠的部分，在构成总的卷积输出时只需将各段线性卷积部分直接连接起来，省掉了输出段的直接相加。

五、实验内容和要求
1、一个连续信号含两个频率分量，经采样得
1,,1,0])125.0(2cos[]125.02sin[)(-=∆++⋅=N n n f n n x Λππ
已知16=N ，f ∆分别为1/16和1/64，观察其频谱；当128=N 时，f ∆不变，其结果有何不同，为什么？
代码：
N=16;
n=0:N-1;
Df=1/16;
Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);
Xk=fft(Xn,N);
subplot(245);
stem(n,abs(Xk));
xlabel('n');ylabel('X(k)');title('N=16,Df=1/16,频谱图'); subplot(241);
stem(n,abs(Xn));
xlabel('n');ylabel('X(n)');title('N=16,Df=1/16,时序图'); Df=1/64;
Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);
Xk=fft(Xn,N);
subplot(246);
stem(n,abs(Xk));
xlabel('n');ylabel('X(k)');title('N=16,Df=1/64,频谱图'); subplot(242);
stem(n,abs(Xn));
xlabel('n');ylabel('X(n)');title('N=16,Df=1/64,时序图');
N=128;
n=0:N-1;
Df=1/16;
Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);
Xk=fft(Xn,N);
subplot(247);
stem(n,abs(Xk));
xlabel('n');ylabel('X(k)');title('N=128,Df=1/16,频谱图'); subplot(243);
stem(n,abs(Xn));
xlabel('n');ylabel('X(n)');title('N=128,Df=1/16,时序图'); Df=1/64;
Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);
Xk=fft(Xn,N);
subplot(248);
stem(n,abs(Xk));
xlabel('n');ylabel('X(k)');title('N=128,Df=1/64,频谱图'); subplot(244);
stem(n,abs(Xn));
xlabel('n');ylabel('X(n)');title('N=128,Df=1/64,时序图');
2、用FFT 分别实现的110,8.0)(≤≤=n n x n
和)6()()(--=n u n u n h 的线性卷积和10点、20点圆周卷积，记录其波形，并说明他们之间的关系。

代码： function [ y ] = circonv( x1,x2,N )
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
if length(x1)>N
error('N should higher than or equal to the length of x1') end
if length(x2)>N
error('')
end
x1=[x1,zeros(1,N-length(x1))];
x2=[x2,zeros(1,N-length(x2))];
m=0:1:N-1;
x2=x2(mod(-m,N)+1);
H=zeros(N,N);
for n=1:1:N
H(n,:)=cirshift(x2,n-1,N);
y=x1*H';
function [ y] = cirshift( x,m,N )
%UNTITLED4 Summary of this function goes here % Detailed explanation goes here
if length(x)>N
error('')
end
x=[x zeros(1,N-length(x))];
n=0:1:N-1;
n=mod(n-m,N);
y=x(n+1);
n=0:11;
x=0.8.^n;
h=[1,1,1,1,1,1];
N=17;
x=[x,zeros(1,5)];
h=[h,zeros(1,11)];
X=fft(x);
H=fft(h);
Y=X.*H;
y=ifft(Y);
subplot(321);stem(x);xlabel('n');ylabel('x'); title('序列x(n)');
subplot(323);stem(h);xlabel('n');ylabel('h'); title('序列h(n)');
subplot(325);stem(y);xlabel('n');ylabel('y'); title('序列y(n),线性卷积');
n=0:11;
N=20;
x=0.8.^n;
h=[1,1,1,1,1,1];
y=circonv(x,h,N);
x=[x,zeros(1,5)];
h=[h,zeros(1,11)];
subplot(322);stem(x);xlabel('n');ylabel('x'); title('序列x(n)');
subplot(324);stem(h);xlabel('n');ylabel('h'); title('序列h(n)');
subplot(326);stem(y);xlabel('n');ylabel('y'); title('序列y(n)，圆周卷积');
频谱图：
说明他们之间的关系？
周期卷积是线性卷积的周期延拓。

六、思考题
1.实验中的信号序列)(n x c 和)(n x d ，在单位圆上的Z 变换频谱)(ωj c e
X 和)(ωj d e X 会相同吗？如果不同，说出哪一个低频分量更多一些，为什么？
答：不同；其中，Xd(n)的低频分量较多，由图形可以看出，在低频处，Xd(n)的取值较多，呈递减趋势。

2. 对一个有限长序列进行DFT 等价于将该序列周期延拓后进行DFS 展开，因为DFS 也只是取其中一个周期来计算，所以FFT 在一定条件下也可以用以分析周期信号序列。

如果实正弦信号1.0,)2sin(=f fn π用16点FFT 来做DFS 运算，得到的频谱是信号本身的真实谱吗？为什么？
答：可以把有限长非周期序列假设为一无限长周期序列的一个主直周期，即对有限长非周期序列进行周期延拓，延拓后的序列完全可以采用DFS 进行处理，即采用复指数。

实验二 IIR 数字滤波器的设计
一、实验目的
1. 掌握双线性变换法及脉冲响应不变法设计IIR 数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR 数字滤波器的MATLAB 编程。

2. 观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点。

3. 熟悉Butterworth 滤波器、Chebyshev 滤波器和椭圆滤波器的频率特性。

二、实验类型
设计型
三、实验仪器
装有MATLAB 语言的计算机
四、实验原理
1． 脉冲响应不变法
用数字滤波器的单位脉冲响应序列)(n h 模仿模拟滤波器的冲激响应)(t h a ,让)(n h 正好等于)(t h a 的采样值，即)()(nT h n h a =，其中T 为采样间隔，如果以)(s H a 及)(z H 分别表示)(t h a 的拉式变换及)(n h 的Z 变换，则
)2(1)(m T
j s H T z H m a e z sT ∑∞-∞==+=π 2．双线性变换法
S 平面与z 平面之间满足以下映射关系：
);(,2121,11211
ωωσj re z j s s T s T z z z T s =+=-+
=+-⋅=-- s 平面的虚轴单值地映射于z 平面的单位圆上，s 平面的左半平面完全映射到z 平面的单位圆内。

双线性变换不存在混叠问题。

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

以低通数字滤波器为例，将设计步骤归纳如下：
1. 确定数字滤波器的性能指标：通带临界频率c f 、阻带临界频率r f 、通带波动δ、阻带内的最小衰减At 、采样周期T 、采样频率s f ；
2. 确定相应的数字角频率，T f c c πω2=；T f r r πω2=；
3. 计算经过预畸的相应模拟低通原型的频率，
)2(2c c tg T ω=Ω，)2
(2r r tg T ω=Ω； 4. 根据Ωc 和Ωr 计算模拟低通原型滤波器的阶数N ，并求得低通原型的传递函数
)(s H a ；
5. 用上面的双线性变换公式代入)(s H a ，求出所设计的传递函数)(z H ；
6. 分析滤波器特性，检查其指标是否满足要求。

五、实验内容和要求
1．用脉冲响应不变法设计一个巴特沃斯数字低通滤波器,要求在0-0.2π内衰耗不大于3dB,在0.6 π –π内衰耗不小于60dB ，采样频率Fs=500 Hz.
代码：
Omegap=0.2*pi*500;
Omegas=0.6*pi*500;
rp=3;
as=60;
fs=500;
wp=Omegap/fs;
ws=Omegas/fs;
[n,Omegac]=buttord(Omegap,Omegas,rp,as,'s');
[b,a]=butter(n,Omegac,'s');
[bz,az]=impinvar(b,a,fs);
w0=[wp,ws];
hx=freqz(bz,az,w0);
[h,w]=freqz(bz,az);
plot(w*fs/(2*pi),abs(h));
grid
xlabel('频率hz')
ylabel('频率响应幅度')
频谱图：
2．分别用脉冲响应不变法和双线性变换法分别设计一个巴特沃斯滤波器数字低通滤波器,使其特性逼近一个巴特沃斯模拟滤波器的性能指标如下：通带截止频率为2π×2000（rad/s）,阻带截止频率为2π×3000（rad/s），通带衰耗不大于3dB,阻带衰耗不小于15dB，采样频率Fs=100000 Hz，观察记录所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。

比较这两种方法的优缺点。

代码（1）：
omegap=2*pi*2000;
omegas=2*pi*3000;
rp=3;
as=15;
fs=10000;
wp=omegap/fs;
ws=omegas/fs;
[n,omegac]=buttord(omegap,omegas,rp,as,'s');
[b,a]=butter(n,omegac,'s');
[bz,az]=impinvar(b,a,fs);
w0=[wp,ws];
hx=freqz(bz,az,w0);
[h,w]=freqz(bz,az);
plot(w*fs/(2*pi),abs(h));
grid
xlabel('频率hz')
ylabel('频率响应幅度')
屏幕图（1）：
代码（2）：
omegap=2*pi*2000;
omegas=2*pi*3000;
rp=3;
as=15;
fs=10000;
wp=omegap/fs;
ws=omegas/fs;
omegap1=2*fs*tan(wp/2);
omegas1=2*fs*tan(ws/2);
[n,omegac]=buttord(omegap1,omegas1,rp,as,'s'); [b,a]=butter(n,omegac,'s');
[bz,az]=bilinear(b,a,fs);
w0=[wp,ws];
hx=freqz(bz,az,w0);
[h,w]=freqz(bz,az);
plot(w*fs/(2*pi),abs(h));
grid
xlabel('频率hz')
ylabel('频率响应幅度')
频谱图（2）：
比较优缺点：
脉冲响应不变法的优点：1，模拟频率到数字频率的转换时线性的。

2，数字滤波器单位脉冲响应的数字表示近似原型的模拟滤波器单位脉冲响应，因此时域特性逼近好
缺点： 会产生频谱混叠现象，只适合带限滤波器 双线性变换法优点： 克服多值映射得关系，可以消除频率的混叠 缺点： 是非线性的，在高频处有较大的失真。

3． 编写滤波器仿真程序，完成对实际采集的心电图信号序列)(n x （具体数据见下面）
的总响应序列)(n y ，可直接调用MATLAB filter 函数实现仿真。

附：人体心电图采样信号在测量过程中往往受到工业高频干扰，所以，必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。

下面的序列就是一个实际心电图信号采样序列样本)(n x ，其中存在高频干扰，实验时，将其作为输入信号，滤除其中的干扰成分。

⎪⎪⎭
⎪
⎪⎬⎫
⎪⎪⎩⎪⎪⎨⎧--------------------------------=0,2,2,2,2,0,0,2,2,0,0,0,4,2,0,0,0,0,0,4,6,6,6,10,12,12,8,4,2,4,32,66,90,84,60,38,16,0,8,12,6,2'6,6,4,4,6,6,4,2,4,6,4,0,2,4)(n x
代码：
clear all;
wp=0.2*pi;
ws=0.3*pi;
rp=1;
rs=15;
[n,wn]=buttord(wp/pi,ws/pi,rp,rs);
[b,a]=butter(n,wn);
N=0.5/(0.02);
figure(1);
freqz(b,a,N);
grid;
xlabel('频率hz');
ylabel('频率响应幅度');
xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38, -60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2, -4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
figure(2);
subplot(2,1,1);
stem(xn,'.');
title('心电图信号采样序列');
yn=filter(b,a,xn);
subplot(2,1,2)
stem(yn,'.');
title('滤波后的心电图信号');
频谱图：
六、思考题
1.双线性变换法中Ω和ω之间的关系是非线性的，在实验中你注意到这种非线性关系了
吗？从那几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？
答：在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。

如以上实验过
程中，采用双线性变化法设计的butter和cheby1数字滤波器，从图中可以看到这
种非线性关系。

2.能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？
答：IIR数字滤波器的设计实际上是求解滤波器的系数和，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。

如果在S平面
上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。

但是
它的缺点是，存在频率混迭效应，故只适用于阻带的模拟滤波器。

实验三 FIR数字滤波器的设计
一、实验目的
1.掌握用窗函数法，频率采样法及优化设计法设计FIR滤波器的原理及方法，熟悉相应的MA TLAB编程。

2.熟悉线性相位FIR滤波器的幅频特性和相频特性。

3.了解各种不同窗函数对滤波器性能的影响。

二、实验类型
设计型
三、实验仪器
装有MATLAB语言的计算机
四、实验原理
线性相位实系数FIR滤波器按其N值奇偶和h(n)的奇偶对称性分为四种：
1．h(n)为偶对称，N为奇数
H(e jω)的幅值关于ω=0，π，2π成偶对称。

2．h(n)为偶对称，N为偶数
H(e jω)的幅值关于ω=π成奇对称，不适合作高通。

3．h(n)为奇对称，N为奇数
H(e jω)的幅值关于ω=0，π，2π成奇对称，不适合作高通和低通。

4．h(n)为奇对称，N为偶数
H(e jω)ω=0、2π＝0，不适合作低通。

FIR滤波器的常用设计方法：
(一) 窗口法
窗函数法设计线性相位FIR滤波器步骤
1. 确定数字滤波器的性能要求：临界频率p ω、st ω，滤波器单位脉冲响应长度N 。

2. 根据性能要求，合理选择单位脉冲响应)(n h 的奇偶对称性，从而确定理想频率响
应)(ω
j d e
H 的幅频特性和相频特性。

3. 求理想单位脉冲响应)(n h d ，在实际计算中，可对)(ω
j d e
H 按M(M 远大于N)点
等距离采样，并对其求IDFT 得)(n h M ，用)(n h M 代替)(n h d 。

4. 选择适当的窗函数)(n w ，根据)()()(n w n h n h d =求所需设计的FIR 滤波器单位脉
冲响应。

5. 求)(ω
j e
H ，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N ，
重复上述设计过程，以得到满意的结果。

窗函数的傅氏变换)(ω
j e W 的主瓣决定了)(ωj e H 的过渡带宽。

)(ωj e W 的旁瓣大小和
多少决定了)(ω
j e
H 在通带和阻带范围内波动的幅度。

常用的几种窗函数有：
(1) 矩形窗 w(n)=R N (n)； (2)
Hanning 窗
；
(3) Hamming 窗
；
(4) Blackmen 窗
；
(5) Kaiser 窗。

式中I o (x)为零阶贝塞尔函数。

（二）频率采样法
频率采样法是从频域出发，将给定的理想频率响应)(ω
j d e
H 加以等间隔采样
然后以此)(k H d 作为实际FIR 数字滤波器的频率特性的采样值)(k H ，即令
由)(k H 通过IDFT 可得有限长序列)(n h
将上式代入到Z 变换中去可得
其中Φ(ω)是内插函数
五、实验内容及步骤
1. N=50，编程并画出矩形窗、汉宁窗、海明窗和布莱克曼窗的时域波形和归一化的幅度谱，并比较各自的主要特点。

代码：
wvtool(boxcar(50),hanning(50),hamming(50),blackman(50)));grid; 频谱图：
并比较各自的主要特点?
矩形窗优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象；汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗．但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降；海明窗的第一旁瓣衰减为一42dB ．海明窗的频谱也是由3个矩形时窗的频谱合成，但其旁瓣衰减速度为20dB/（10oct ），这比汉宁窗衰减速度慢；布莱克曼窗的幅度函数主要由五部分组成，他们的位移都不同，其幅度也是不同的WRg （w ）使旁瓣再进一步抵消。

旁瓣峰值幅度进一步增加，其幅度谱主瓣宽度是矩形窗的3倍。

设计程序时用backman 函数调用。

2. N=15，带通滤波器的两个通带边界分别是πω
3.01=，πω5.02=。

用汉宁（Hanning ）
窗设计此线性相位带通滤波器，观察它的实际3dB 和20dB 带宽。

N=45，重复这一设计，观察幅频和相位特性的变化，注意长度N 变化的影响。

N=15汉宁代码：
clc;clear all
n=15;
w1=0.3;w2=0.5; wn=[w1,w2];
b2=fir1(n,wn,hanning(n+1)); freqz(b2,1);
title('汉宁窗，N=15'); N=15汉宁频谱：
N=45汉宁代码：
clc;clear all
n=45;
w1=0.3;w2=0.5;
wn=[w1,w2];
b2=fir1(n,wn,hanning(n+1));
freqz(b2,1);
title('汉宁窗，N=45');
N=45汉宁频谱：
3.分别改用矩形窗和Blackman窗，设计2中的带通滤波器，观察并记录窗函数对滤波器幅频特性的影响，比较三种窗的特点。

N=15矩形代码：
clc;clear all
n=15;
w1=0.3;w2=0.5;
wn=[w1,w2];
b2=fir1(n,wn,boxcar(n+1));
freqz(b2,1);
title('矩形窗，N=15');
N=15矩形频谱：
N=45矩形代码：
clc;clear all
n=45;
w1=0.3;w2=0.5;
wn=[w1,w2];
b2=fir1(n,wn,boxcar(n+1)); freqz(b2,1);
title('矩形窗，N=45');
N=45矩形频谱：
N=15Blackman代码：
clc;clear all
n=15;
w1=0.3;w2=0.5;
wn=[w1,w2];
b2=fir1(n,wn,Blackman(n+1)); freqz(b2,1);
title('Blackman窗，N=15'); N=15Blackman频谱：
N=45Blackman代码：
clc;clear all
n=45;
w1=0.3;w2=0.5;
wn=[w1,w2];
b2=fir1(n,wn,Blackman(n+1)); freqz(b2,1);
title('Blackman窗，N=45');
N=45Blackman 频谱：
比较三种窗的特点?
矩形窗优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象；汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗．但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降；布莱克曼窗的幅度函数主要由五部分组成，他们的位移都不同，其幅度也是不同的WRg （w ）使旁瓣再进一步抵消。

旁瓣峰值幅度进一步增加，其幅度谱主瓣宽度是矩形窗的3倍。

设计程序时用backman 函数调用
六、思考题
1. 定性地说明用本实验程序设计的FIR 滤波器的3dB 截止频率在什么位置？它等于理想频率响应)(ωj d e H 的截止频率吗？
答：3dB 截止频率在最高幅度下降一半是所对应的频率。

因为设计的时候总会有误差，所以并不等于理想Hd(ejw)的截止频率。

2.如果没有给定h(n)的长度N ，而是给定了通带边缘截止频率p ω和阻带临界频率st ω，以及相应的衰减，能根据这些条件用窗函数法设计线性相位FIR 低通滤波器吗？
答：可以，可以用雷米兹交替法来设计。