数字信号处理实验报告

    姓名：孔德权

            学院：电信学院

            学号：101045113

            指导老师：韩萍

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

1、实验目的

（1）  熟悉用双线性变换法设计IIR数字滤波器的原理与方法。

（2）  掌握数字滤波器的设计方法。

（3）  通过观察对实际心电图信号的滤波过程，获得数字滤波的感性认识。

2、实验原理

   基于模拟滤波器数学模型--微分方程的积分，利用对积分的数值逼近可以得到相应的数字滤波器。假设模拟滤波器的系统函数

对模拟滤波器取样后逼近后，得到的数字滤波器系统函数

可见，数字滤波器可通过下述变换关系有模拟滤波器得到

上述s平面到z平面之间的变换是双线性变换关系，故称之为“双线性”变换。它对高阶系统同样成立。

3、实验内容

（1）  用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计参数为:在低通频带内频率低于时,最大衰减小于1dB；在阻带内[，]频率区间上，最小衰减大于15dB。

（2）  以为采样间隔，打印出数字滤波器在频率区间[0，]上的幅频响应特性曲线。

（3）  用所设计的滤波器对实际心电图信号采样序列进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。

4、实验步骤

（1）  复习有关巴特沃斯模拟滤波器的设计和用双线性变换法设计IIR数字滤波器的内容，即变换关系。根据教材中p194页例题中给定的设计指标，调用MATLAB信号处理工具箱函数buttord和butter设计出相应的数字滤波器H(Z).

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

[附：人体心电图采样信号在测量过程中往往受到工业高频干扰，所以，必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面的序列就是一个实际心电图信号采样序列样本，其中存在高频干扰，实验时，将其作为输入信号，滤除其中的干扰成分。

（3）  在计算机上运行仿真程序，并调用绘图子程序，完成实验内容（2）和（3）。

5、实验过程

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);

title('滤波器的扶贫响应曲线');

xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-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,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,0];

figure(2);

subplot(2,1,1);

stem(xn,'.');

title('心电图信号采样序列x(n)');

yn=filter(b,a,xn);

subplot(2,1,2)

stem(yn,'.');

title('滤波后的心电图信号');

 6、实验结果

 

  7、思考题

用双线性变换法设计数字滤波器过程中,变换公式中T的取值对设计结果有无影响？为什么？

答: T的取值对设计结果没有影响，因为不管是设计模拟滤波器还是由模滤波器转换为数字滤波器的过程中系数2/T均被约分掉。

实验四   FIR数字滤波器的设计及滤波实现

1、实验目的

(1)掌握用窗函数法设计FIR数字滤波器的原理和方法。

(2)熟悉线性相位FIR数字滤波器特性。

(3)了解各种窗函数对滤波特性的影响。

2、实验原理

如果所希望的滤波器的理想频率响应函数为，则其对应的单位脉冲响应为

窗函数设计方法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，且是因果的，所以用窗函数将截断，并进行加权处理，得到

就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为

式中，N为所选窗函数的长度。

这里需要注意的是：设计出的滤波器性能取决于窗函数的类型和窗口的长度N的取值，设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度（各类窗函数的频率特性见教材的第三章）。

选定窗函数和长度N后，求出单位脉冲响应和。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整次幂，以便用FFT计算，如果观察细节，补领点数增多即可。如果不满足设计要求，则要重新选择窗函数类型和长度N,再次验算，直至满足要求。

如果要求线性相位特性，则还必须满足

根据上式中的正、负和长度N的奇、偶性又将线性相位FIR滤波器分为四类，要根据所设计的滤波特性正确选择其中一类。例如要设计线性相位低通特性，可选择一类，而不能选一类。

3、实验步骤

（1）    复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握

设计步骤。

（2）          编写程序

(a)       编写能产生矩形窗、升余弦窗的窗函数子程序。

(b)   编写主程序。主程序框图如下（仅供参考），其中幅度特性要求用dB表示。

设

画图时，用打印幅度特性。

第点对应的频率。为使曲线

包络更接近的幅频特性，DFT变换

的区间要选得大一些。例如窗口长度N=33

时，可通过在末尾补零的方法，使长

度变为64，再进行64点的DFT，则可得到

更精确的幅度衰减特性曲线。

4、上机实验内容

（1）  用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率，窗口长度N=15,33。要求在两种窗口长度情况下，分别求出，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽。总结窗口长度N对滤波特性的影响。

      设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数，即

其中，

（2）  N=33, ，用两种窗函数设计线性相位低通滤波器。绘制相应的幅频特性曲线，观察3dB和20dB带宽以及阻带最小衰减，比较两种窗函数的滤波特性的影响。

5、实验过程

clear all;

N=15;

wp=0.25*pi;

window=hanning(N);

b=fir1(N-1,wp/pi,window);

figure(1);

freqz(b,1,1024);

title('衰减特性(db)');

figure(2);

subplot(2,1,1);

stem(b,'.');

title('实际低通滤波器');

subplot(2,1,2);

y=freqz(b,1,1024);

plot(abs(y));

title('幅频特性');

N=33;

wp=0.25*pi;

window=hanning(N);

b=fir1(N-1,wp/pi,window);

figure(3);

freqz(b,1,1024);

title('衰减特性(db)');

figure(4);

subplot(2,1,1);

stem(b,'.');

title('实际低通滤波器');

subplot(2,1,2);

y=freqz(b,1,1024);

plot(abs(y));

title('幅频特性');

N=15;

wp=0.25*pi;

window=hamming(N);

b=fir1(N-1,wp/pi,window);

figure(5);

freqz(b,1,1024);

title('衰减特性(db)');

figure(6);

subplot(2,1,1);

stem(b,'.');

title('实际低通滤波器');

subplot(2,1,2);

y=freqz(b,1,1024);

plot(abs(y));

title('幅频特性');

N=33;

wp=0.25*pi;

window=hamming(N);

b=fir1(N-1,wp/pi,window);

figure(7);

freqz(b,1,1024);

title('衰减特性（db）');

figure(8);

subplot(2,1,1);

stem(b,'.');

title('实际低通滤波器');

subplot(2,1,2);

y=freqz(b,1,1024);

plot(abs(y));

title('幅频特性');

N=15;

wp=0.25*pi;

window=blackman(N);

b=fir1(N-1,wp/pi,window);

figure(9);

freqz(b,1,1024);

title('衰减特性(db)');

figure(10);

subplot(2,1,1);

stem(b,'.');

title('实际低通滤波器');

subplot(2,1,2);

y=freqz(b,1,1024);

plot(abs(y));

title('幅频特性');

N=33;

wp=0.25*pi;

window=blackman(N);

b=fir1(N-1,wp/pi,window);

figure(11);

freqz(b,1,1024);

title('衰减特性(db)');

figure(12);

subplot(2,1,1);

stem(b,'.');

title('实际低通滤波器');

subplot(2,1,2);

y=freqz(b,1,1024);

plot(abs(y));

title('幅频特性');

N=15;

wp=0.25*pi;

window=boxcar(N);

b=fir1(N-1,wp/pi,window);

figure(13);

freqz(b,1,1024);

title('衰减特性(db)');

figure(14);

subplot(2,1,1);

stem(b,'.');

title('实际低通滤波器');

subplot(2,1,2);

y=freqz(b,1,1024);

plot(abs(y));

title('幅频特性');

N=33;

wp=0.25*pi;

window=boxcar(N);

b=fir1(N-1,wp/pi,window);

figure(15);

freqz(b,1,1024);

title('衰减特性(db)');

figure(16);

subplot(2,1,1);

stem(b,'.');

title('实际低通滤波器');

subplot(2,1,2);

y=freqz(b,1,1024);

plot(abs(y));

title('幅频特性');

6、实验结果

（1）汉宁窗：

窗口长度N=15

窗口长度N=33

（2）汉明窗：

窗口长度N=15

窗口长度N=33

（3） 布莱克曼窗：

窗口长度N=15

窗口长度N=33

（4） 矩形窗：

窗口长度N=15

窗口长度N=33

7、思考题

（1）  如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器？写出步骤.

答：步骤如下：

① 首先确定模拟带通滤波器的技术指标。

② 确定归一化低通技术要求

③ 设计归一化低通

      clear all;

N=15;                        %N为窗口长度

wp=0.25*pi;                 %wp为截止频率

window=hanning(N);

b=fir1(N-1,wp/pi,window);%频率归一化

figure(1);

freqz(b,1,1024);           %衰减特性显示

（2）  如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率和 ，试求理想带通的单位脉冲响应

答:通过和求出这样便可以利用公式=sin[Wc(n-a)]/[pi*(n-a)]求出。

 

第二篇：地震子波 数字信号实验报告

地震子波波形显示及一维地震合成记录

姓名：杨肖迪       学号：050422009040     专业：地球信息科学与技术  2009级

实验目的

1.   认识地震子波（以雷克子波为例），对子波的波形有直观的认识。

2.     利用褶积公式合成一维地震记录。

实验步骤

1.     雷克子波

  （零相位子波）

   （最小相位子波）

其中代表子波的中心频率，代表子波宽度,随着的增大，子波能量后移，当=7时，最小相位子波可视为混合相位子波，这里 = 25 Hz，= 4;

零相位子波与最小相位子波是同一家族的子波。对于给定的振幅谱，只要配上不同的相位谱，便得到不同的子波。但零相位子波的能量聚集在首部，即从开始时就具有最大能量，无短暂积累过程。而最小相位子波的能量聚集在序列首部，最大能量有短暂积累过程。

2.   根据公式编程实现零相位子波、最小相位子波的波形显示；

零相位子波：

时间序列长度为200，时间采样间隔为0.001s，

生成雷克子波的程序为：

for i=1:200

w1(i)=exp((-4*pi^2*f^2/r^2)*i^2*dt^2)*cos(2*pi*f*i*dt);

end

其中子波中心频率f=25HZ，子波宽度r=4，子波周期为1/f=0.04s，对应图中40位置，即t=40*1.001=0.04s。能量聚集在首部，开始时就具有最大能量，无积累过程。其相位和振幅为：

当振幅最大时，相位为零，即此时波的振幅为实数，达到最大值。

最小相位子波：

时间序列长度为200，时间采样间隔为0.001s，

生成雷克子波程序为：

for i=1:200

        w2(i)=exp((-4*pi^2*f^2/r^2)*i^2*dt^2)*sin(2*pi*f*i*dt);

end

其中子波中心频率f=25HZ，子波宽度r=4，子波周期为1/f=0.04s，对应图中40位置，即t=40*1.001=0.04s。能量聚集在序列首部，是最小能量延迟的。

其相位和振幅为：

从中可以看出，当振幅最大时，相位不为零，是非零相位的，相对零相位子波来说，最大能量是延迟的。

当 ，时最小相位子波可视为混合相位子波，如图所示：

混合相位子波的能量聚集在序列中部，是混合能量延迟的。

3.     设计反射系数(n=500)，n为地层深度，其中，为第一层介质深度；，为第二层介质深度；，为第三层介质深度；，为第四层介质深度；，为第五层介质深度；其它为0。

其地质模型为：

地震波在介质中传播，当到达介质分界面时，发生反射和透射，反射波被检波器接受，生成地震记录。反射系数表示地震波在两层介质分界面的能量重新分配，如r(100)=1.0，表示地震波入射到分界面时，只有一种波，反射纵波（或反射横波）。反射系数不为1.0时，表示当地震波入射到分界面时，产生两种反射波。反射系数为正，表示反射波相位与入射波相位相差2π；反射系数为负，表示反射波相位与入射波相位相差π，存在半波损失。

4.     应用褶积公式合成一维地震记录，并图形显示；

应用褶积公式求f（n）的程序为：

for i=1:n+m-1

    for j=1:m

        if i-j>0&i-j<=200

            a1(j)=r(j)*w1(i-j);

            a2(j)=r(j)*w2(i-j);

            f1(i)=f1(i)+a1(j);

            f2(i)=f2(i)+a2(j);

        end

    end

end

其中f(n)长度为n+m-1，n为w(n)雷克子波时间序列长度，m为r(m)反射系数序列的长度。线性褶积f(n)的长度为雷克子波长度与地质模型长度之和减一，即f(n)=n+m-1。f1为零相位子波w1(n)与反射系数r(n)的摺积；f2为最小相位子波w2(n)与反射系数r(n)的摺积。

实验结果

一维地震记录图形显示为：

零相位子波：

由记录中可以看到，当n=100、200、300、400、500时反射波的能量是最大的。因为当r(m)=1.0时，即j=100时，i=100时，w(i-j)=1.0，是能量最大的，即w(0)=1.0。同理，当n=200、300、400、500时能量也是最大的。

最小相位子波：

由地震记录中可以看出，当n=100、200、300、400、500时，能量不是最大的，而是延迟短暂时间后出现的，因为当r(m)=1.0时，即j=100时，但i=100时，w(i-j)不是最大能量的，即w(0)≠1.0，最大值是延迟出现。同理，当n=200、300、400、500时，能量也不是最大的，而是要延迟出现。也说明最小相位的子波相比零相位子波来说是能量延迟的。

实验分析

根据所学知识对实验结果进行分析；

地震子波由震源激发，在地层中传播，因为在沉积地层中，每层介质的物理性质不相同，从而使得地震波的传播速度也不相同。当地震波传播到两层介质的分界面时，会发生反射，由于每层介质的反射系数不同，所以反射波的能量也不相同，检波器接收到不同时刻的、不同能量的反射波，形成一个地震记录。

由合成地震记录中可以看出，最小相位子波相对零相位子波来说是相位滞后的，能量延迟的，但两者为同一家族的子波。合成地震记录中横坐标为时间，纵坐标为振幅。每一时刻的值由m个值的和组成，m为反射系数r(n)的长度，整个地震记录由m+n-1个时刻的值组成。对于零相位的地震记录来说，当r(m)=1.0时，即j=100时，i=100时，w(i-j)=1.0，是能量最大的，即w(0)=1.0。同理，当n=200,、300、400、500时，能量也是最大的。对于最小相位的地震记录来说，当r(m)=1.0时，即j=100时，但i=100时，w(i-j)不是最大能量的，即最大能量不是在w(0)出现，而是延迟出现。同理，当n=200、300、400、500时，能量也不是最大的，而是要延迟出现。

附：源程序代码

clear

f=25;  %地震子波中心频率

r=4;   %子波宽度

n=200;  %时间序列长度

dt=0.001;  %时间域采样间隔

for i=1:200

    w1(i)=exp((-4*pi^2*f^2/r^2)*i^2*dt^2)*cos(2*pi*f*i*dt);%零相位子波

    w2(i)=exp((-4*pi^2*f^2/r^2)*i^2*dt^2)*sin(2*pi*f*i*dt);%最小相位子波

end

m=500;%地层深度

for i=1:m

    r(i)=0;%地层性质´

end

r(100)=1.0;%第一层介质反射系数

r(200)=-0.7;%第二层介质反射系数

r(300)=0.5;%第三层介质反射系数

r(400)=0.4;%第四层介质反射系数

r(500)=0.6;%第五层介质反射系数

for i=1:n+m-1 %线性褶积长度

    f1(i)=0;

    f2(i)=0;

end

for i=1:n+m-1%地震记录长度

    for j=1:m %地层深度

        if i-j>0&i-j<=200 %满足雷克子波的时间序列长度

            a1(j)=r(j)*w1(i-j);

            a2(j)=r(j)*w2(i-j);

            f1(i)=f1(i)+a1(j);

            f2(i)=f2(i)+a2(j);

        end

    end

end

subplot(2,2,1),plot(w1,'k'),axis([0,200,-1,1]),title('零相位子波')

subplot(2,2,2),plot(w2,'k'),axis([0,200,-1,1]),title('最小相位子波')

subplot(2,2,3),plot(f1,'k'),title('零相位子波地震记录')

subplot(2,2,4),plot(f2,'k'),title('最小相位子波地震记录')