《测试信号分析与处理》(附实验结果).doc

《测试信号分析与处理》实验指导书
实验一差分方程、卷积、z变换
一、实验目的
通过该实验熟悉 matlab软件的基本操作指令，掌握matlab软件的使用方法，掌握数字信号处理中的基本原理、方法以及matlab函数的调用。

二、实验设备
1、微型计算机1台；
2、matlab软件1套
三、实验原理
Matlab 软件是由mathworks公司于1984年推出的一套科学计算软件，分为总包和若干个工具箱，其中包含用于信号分析与处理的sptool工具箱和用于滤波器设计的fdatool工具箱。

它具有强大的矩阵计算和数据可视化能力，是广泛应用于信号分析与处理中的功能强大且使用简单方便的成熟软件。

Matlab软件中已有大量的关于数字信号处理的运算函数可供调用，本实验主要是针对数字信号处理中的差分方程、卷积、z变换等基本运算的matlab函数的熟悉和应用。

差分方程（difference equation)可用来描述线性时不变、因果数字滤波器。

用x表示滤波器的输入，用y表示滤波器的输出。

a0y[n]+a1y[n-1]+…+a N y[n-N]=b0x[n]+b1x[n-1]+…+b M x[n-M] (1)
ak,bk 为权系数，称为滤波器系数。

N为所需过去输出的个数，M 为所需输入的个数卷积是滤波器另一种实现方法。

y[n]= ∑x[k] h[n-k] = x[n]*h[n] (2) 等式定义了数字卷积，*是卷积运算符。

输出y[n] 取决于输入x[n] 和系统的脉冲响应h[n]。

传输函数H(z)是滤波器的第三种实现方法。

H(z)=输出/输入= Y(z)/X(z) （3）即分别对滤波器的输入和输出信号求z变换，二者的比值就是数字滤波器的传输函数。

序列x[n]的z变换定义为
X (z)=∑x[n]z-n (4) 把序列x[n] 的z 变换记为Z{x[n]} = X(z)。

由X(z) 计算x[n] 进行z 的逆变换x[n] = Z-1{X(z)}。

Z 变换是Z-1的幂级数，只有当此级数收敛，Z 变换才有意义，而且同一个Z 变换等式，收敛域不同，可以代表不同序列的Z 变换函数。

这三种数字滤波器的表示方法之间可以进行相互转换。

四、实验步骤
1、熟悉matlab软件基本操作指令。

读懂下列matlab程序指令，键入程序
并运行，观察运行结果。

Conv.m% 计算两个序列的线性卷积；
%-----------------------------------------------------------------
clear;
N=5;
M=6;
L=N+M-1;
x=[1,2,3,4,5];
h=[6,2,3,6,4,2];
y=conv(x,h);
nx=0:N-1;
nh=0:M-1;
ny=0:L-1;
subplot(231);
stem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on;
subplot(232);
stem(nh,h,'.k');xlabel('n');ylabel('h(n)');grid on;
subplot(233);
stem(ny,y,'.k');xlabel('n');ylabel('y(n)');grid on;
filter.m;%求一个离散系统的输出；
clear;
x=ones(100);
t=1:100;
b=[.001836,.007344,.011016,.007374,.001836];
a=[1,-3.0544,3.8291,-2.2925,.55075];
y=filter(b,a,x);
clear;
impz .m% 计算滤波器的冲击响应
b=[.001836,.007344,.011016,.007374,.001836];
a=[1,-3.0544,3.8291,-2.2925,.55075];
[h,t]=impz(b,a,40);
subplot(221)
stem(t,h,'.');grid on;
ylabel('h(n)')
xlabel('n')
filter.m% 计算滤波器的阶跃响应
x=ones(100);t=1:100;
y=filter(b,a,x);
subplot(222)
plot(t,x,'g.',t,y,'k-');grid on;
ylabel('x(n) and y(n)')
xlabel('n')
2、编程求出下列问题的解
1）、滤波器的差分方程为：y[n]=x[n]-0.8x[n-1]-0.5y[n-1]
求出此滤波器脉冲响应和阶跃响应的前十个采样值。

2）、系统的脉冲响应为h[n]=e-n(u[n]-u[n-3]),用卷积求系统的阶跃响应。

五、实验讨论和分析
1、差分方程、卷积、z变换和傅里叶变换之间如何进行转换？
2、边界效应是如何产生的？它对信号的滤波效果有何影响？
实验二数字滤波器综合设计
一、实验目的
通过该设计实验掌数字滤波器设计的一般步骤，掌握利用matlab 软件设计数字滤波器的方法，熟悉sptool工具箱的使用方法。

二、实验设备
1、微型计算机1台；
2、matlab软件1套
三、实验原理
一）、滤波器的形状及重要参数
理想滤波器的形状是矩形，图 1 给出非理想滤波器。

图1
通带：增益高的频率范围，信号可以通过，称为滤波器的通带。

阻带：增益低的频率范围，滤波器对信号有衰减或阻塞作用，称滤波器的阻带。

滤波器截止频率：增益为最大值的0.707倍时所对应的频率为滤波器截止频率增益通常用分贝（dB）表示。

增益（dB）= 20log(增益）
增益为0.707 时对应-3dB，因此截止频率常被称为-3dB。

滤波器的带宽：对于低通滤波器宽带是从0 ~ - 3dB
对于高通滤波器宽带是从- 3dB~采样频率的一半
对于带通滤波器带宽是截止频率之间的频率距离
二）加窗低通FIR 滤波器的设计
1. 在过渡带宽度的中间，选择通带边缘频率(Hz):
f1=所要求的通带边缘频率+(过渡带宽度)/2
2. 计算Ω1=2πf1/fs，并将此值代入理想低通滤波器的脉冲响应h1[n] 中：
h1[n] = sin(nΩ1)/nπ
3. 从表中选择满足阻带衰减及其他滤波器要求的窗函数，用表中N 的公式
计算所需要的非零项数目。

选择奇数项，这样脉冲响应可以完全对称，避免了滤波器产生相位失真，对于|n|≤(N-1)/2，计算窗函数w[n]。

4. 对于|n|≤(N-1)/2，从式h[n]=h1[n]w[n]计算（有限）脉冲响应，对于其他
n 值h[n]=0，此脉冲响应是非因果的。

5. 将脉冲响应右移(N-1)/2，确保第一个非零值在n=0处，使此低通滤波器
为因果的。

三）、设计低通巴特沃斯滤波器：
1) 确定待求通带边缘频率fp1 Hz 、待求阻带边缘频率fs1 Hz 和待求阻
带衰减- 20logδsdB(或待求阻带增益20logδsdB)。

通带边缘频率对应–3dB增益。

2) 用式Ω=2πf/fs 把由Hz 表示的待求边缘频率转成由弧度表示的数
字频率，得到Ωp1 和Ωs1 。

3) 计算预扭曲模拟频率以避免双线性变化带来的失真。

由ω=2fs tan(Ω/2)
求得ωp1和ωs1，单位是弧度/秒。

4) 由已给定的阻带衰减- 20logδs(或增益- 20logδs)确定阻带边缘增益
δs 。

5) 计算所需滤波器的阶数n 取整数。

6）把ωp1代入n 阶模拟巴特沃斯滤波器传输函数H(s)中，并对H(s) 进行双线性变换得到n 阶数字传输函数H(z)。

滤波器实现所需的差分方程可直接从传输函数H(s) 求出。

四）、低通切比雪夫Ⅰ型滤波器的设计：
1）确定待求的通带与阻带边缘频率fp1 和fs1 、待求的通带边缘增益20log(1- δp) 和待求的阻带衰减-20logδs（或待求的阻带增益20logδs ）。

2）用公式Ω=2πf/fs 将待求的边缘频率转换为数字频率（用弧度表示），得到Ωp1 和Ωs1 。

3）对数字频率采用预扭曲以避免双线性变换引起的误差。

由ω=2fs tan(Ω/2) 得到ωp1和ωs1，单位是弧度/秒。

4）由指定的通带边缘增益20log(1- δp) ，确定通带边缘增益1- δp 。

计算参数ε。

5）由指定的衰减-20logδs（或增益20logδs），确定阻带边缘增益δs 。

6）计算所需的阶数n。

7）将ωp1 和δp 代入n 阶模拟切比雪夫Ⅰ型滤波器的传输函数H(s)，并对其进行双线性变换，得到n 阶数字滤波器传输函数H(z)。

实现滤波器所需的差分方程可由传输函数H(z) 直接得到。

四、实验步骤
1、任选第9、10章后滤波器设计题各2题，利用matlab编程完成滤
波器的设计，并画出滤波器的脉冲响应、幅度响应和相位响应图。

五、实验讨论和分析
1、设计得到的滤波器与设计要求有无差别？如果有，请分析误差产生的原因。

2、FIR滤波器与IIR滤波器的优缺点分别是什么？针对具体信号进行滤波时，
如何选择？
实验一
clear;
N=5;
M=6;
L=N+M-1;
x=[1,2,3,4,5];
h=[6,2,3,6,4,2];
y=conv(x,h);
nx=0:N-1;
nh=0:M-1;
ny=0:L-1;
subplot(231);
stem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on; subplot(232);
stem(nh,h,'.k');xlabel('n');ylabel('h(n)');grid on; subplot(233);
stem(ny,y,'.k');xlabel('n');ylabel('y(n)');grid on;
clear;
x=ones(100);
t=1:100;
b=[.001836,.007344,.011016,.007374,.001836]; a=[1,-3.0544,3.8291,-2.2925,.55075];
y=filter(b,a,x);
clear;
b=[.001836,.007344,.011016,.007374,.001836];
a=[1,-3.0544,3.8291,-2.2925,.55075];
[h,t]=impz(b,a,40);
subplot(221)
x=ones(100);t=1:100;
y=filter(b,a,x);
subplot(222)
plot(t,x,'g.',t,y,'k-');grid on;
ylabel('x(n) and y(n)')
xlabel('n')
2、编程求出下列问题的解
1）、滤波器的差分方程为：y[n]=x[n]-0.8x[n-1]-0.5y[n-1] 求出此滤波器脉冲响应和阶跃响应的前十个采样值。

syms x
x0=2;
f=x^3-3*x-1;
eps=1e-6;
maxcnt=1000;
fx=diff(f,x);
x1=x0;
cnt=1;
while cnt<=maxcnt
x2=x1-subs(f/fx,x,x1);
if abs(x1-x2)<eps
break;
end
[cnt,x1,x2] %迭代次数，迭代前，迭代后
x1=x2;
cnt=cnt+1;
end
subs(f,x,x2)
输入matlab软件的输出：
syms x
x0=2;
f=x^3-3*x-1;
eps=1e-6;
maxcnt=1000;
fx=diff(f,x);
x1=x0;
cnt=1;
while cnt<=maxcnt
x2=x1-subs(f/fx,x,x1);
if abs(x1-x2)<eps
break;
end
[cnt,x1,x2] %迭代次数，迭代前，迭代后
x1=x2;
cnt=cnt+1;
end
subs(f,x,x2)
ans =
1.0000
2.0000 1.8889
ans =
2.0000 1.8889 1.8795
ans =
3.0000 1.8795 1.8794
ans =
8.8818e-016
2）、系统的脉冲响应为h[n]=e-n(u[n]-u[n-3]),用卷积求系统的阶跃响应。

M=3;
N=10;
L=N+M-1;
x=[1,1,1,1,1,1,1,1,1,1];
h=[1,exp(-1),exp(-2)];
y=conv(x,h);
nx=0:N-1;
nh=0:M-1;
ny=0:L-1;
subplot(131);
stem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on; subplot(132);
stem(nh,h,'.k');xlabel('n');ylabel('h(n)');grid on; subplot(133);
stem(ny,y,'.k');xlabel('n');ylabel('y(n)');grid on;
实验二
习题9.15
f1=4000;%信号频率Hz
f2=5000;%信号频率Hz
f3=6000;%信号频率Hz
fs=12000;%采样频率Hz
N=32;%采样点数
t=(0:N-1)/fs;%采样时间
x1=sin(2*pi*f1*t);%信号采样值x2=sin(2*pi*f2*t);%信号采样值x3=sin(2*pi*f3*t);%信号采样值
x=x1+x2+x3;
y=filter(h,1,x);
f1=3000+250;
fs=12000;
w=2*f1/fs;
n=3.32*fs/500;
h=makelp(n,w,'hanning'); [mag,phase,w]=dtft(h);
plot(t,x,'g',t,y,'k-')
习题9.20：
h=makelp(19,0.625,'hamming'); [mag,phase,w] = dtft(h); plotdtft(mag,phase,w,2)
stem( 0:18,h,'.')
ylabel('h(n)');
xlabel('n');
习题9.23：
h=bandfilt(59,0.31875,0.68125,1,'hanning'); [mag,phase,w]=dtft(h);
plotdtft(mag,phase,w,2);
stem(0:116,h,'.');
ylabel('h(n)');
xlabel('n');
习题10.12
n=buttord(0.25,0.375,3,44); [b,a]=butter(n,0.25); [mag,phase,w] = dtft(b,a); plotdtft(mag,phase,w,1); [h,t]=impz(b,a,40); subplot(111)
stem(t,h,'.');grid on;
ylabel('h(n)')
xlabel('n')
习题10.15
n = cheb1ord(0.64,0.72,3,21.94); [b,a]=cheby1(n,3,0.64); [mag,phase,w] = dtft(b,a); plotdtft(mag,phase,w,1);
[h,t]=impz(b,a,40);
>> subplot(111)
stem(t,h,'.');grid on;
ylabel('h(n)')；
xlabel('n')；。