数字信号处理实验.doc

实验二 DFT用于频谱分析
(一)、在运用DFT进行频谱分析的过程中可能产生三种误差：
  (1) 混叠
    序列的频谱时被采样信号的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。
避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。
(2) 泄漏
实际中我们往往用截短的序列来近似很长的甚至是无限长的序列，这样可以使用较短的DFT来对信号进行频谱分析，这种截短等价于给原信号序列乘以一个矩形窗函数，也相当于在频域将信号的频谱和矩形窗函数的频谱卷积，所得的频谱是原序列频谱的扩展。
泄漏不能与混叠完全分开，因为泄漏导致频谱的扩展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减至最小。 
DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数，就一定意义上看，用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住，不能别我们观察到。
    减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值，从而变动DFT的点数，这一方法实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或谷点暴露出来。
用FFT可以实现两个序列的圆周卷积。在一定的条件下，可以使圆周卷积等于线性卷积。一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于线性卷积的充要条件是FFT的长度
N≥N1＋N2
对于长度不足N的两个序列，分别将他们补零延长到N。
当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。有两种方法：
重叠相加法。将长序列分成与短序列相仿的片段，分别用FFT对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。
重叠保留法。这种方法在长序列分段时，段与段之间保留有互相重叠的部分，在构成总的卷积输出时只需将各段线性卷积部分直接连接起来，省掉了输出段的直接相加。
(三)、用周期图法(平滑周期图的平均法)对随机信号作谱分析
实际中许多信号往往既不具有有限能量，由非周期性的。
无限能量信号的基本概念是随机过程，也就是说无限能量信号是一随机信号。周期图法是随机信号作谱分析的一种方法，它特别适用于用FFT直接计算功率谱的估值。
将长度为N的实平稳随机序列的样本x(n)再次分割成K段，每段长度为L,即L=N/K。每段序列仍可表示为：
xi(n)=x(n+(i-1)L)，0≤n≤L-1，1≤i≤K
但是这里在计算周期图之前，先用窗函数w(n)给每段序列xi(n)加权，K个修正的周期图定义为

其中U表示窗口序列的能量，它等于：

在此情况下，功率谱估计量可表示为：

三、实验内容及步骤
    实验中用到的信号序列：

a)  Gaussian序列

b) 衰减正弦序列

c) 三角波序列

d)反三角波序列

上机实验内容
(1)观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。
(2)观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。
(3)观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列xc(n)和xd(n)的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。在xc(n)和xd(n)末尾补零，用N=16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFT频谱还有相同之处吗？这些变化说明了什么？
(4)一个连续信号含两个频率分量，经采样得
x(n)=sin2π*0.125n+cos2π*(0.125+Δf)n n=0,1……,N-1 已知N=16，Δf分别为1/16和1/64，观察其频谱；当N=128时，Δf不变，其结果有何不同，为什么？
(5)用FFT分别实现xa(n)（p＝8，q＝2）和 xb(n)（a＝0.1，f＝0.0625）的16点圆周卷积和线性卷积。
(6)产生一512点的随机序列xe(n)，并用xc(n)和xe(n)作线性卷积，观察卷积前后xe(n)频谱的变化。要求将xe(n)分成8段，分别采用重叠相加法和重叠保留法。
四，实验程序及结果图像
1, 高斯序列
N=16;
n=0:1:15;
p=8;
q=2;
%q=4;
%q=8;
a=0.1;
f=0.0625;
xa=exp(-((n-p).^2)./q);
xb=exp(-a*n).*sin(2*pi*f*n);
figure(1)
stem(n, xa,'.');
title('xa(n)序列')
xlabel('n')
ylabel('xa(n)')
grid on

[H, w] = freqz(xa, 1, [], 'whole', 1);
Hamplitude = abs(H);
Hphase = angle(H);
Hphase = unwrap(Hphase);
figure(2)
subplot(2, 1, 1)
plot(w, Hamplitude)
title('幅频响应')
xlabel('w/(2*pi)')
ylabel('|H(exp(jw))|')
grid on
subplot(2, 1, 2)
plot(w, Hphase)
title('相频响应')
xlabel('w/(2*pi)')
ylabel('fai(H(exp(jw)))')
grid on

(1) P=8，q=2时的幅频特性与相频特性曲线

（2） p=8，q=4时的幅频特性与相频特性曲线

（3）P=8,q=8时的幅频特性与相频特性曲线

(4) p=13,q=8时的幅频特性与相频特性曲线

（5） p=14,q=8时的幅频特性与相频特性曲线

(6) 思考题回答
高斯序列可以作为低通滤波器使用，通过调整p，q的值，可以调整滤波器参数，其中，p决定时域中心位置，q决定点数与幅度，频域上，通过取不同的q值，可以得到不同的通带增益与过渡带宽度。
p固定时，高斯序列的中心位置固定，q的变化会影响高斯序列时域上的点数与幅度，q越大，点数越多，幅度越大，频域上，q越大，通带增益越大，过渡带越窄。
q固定为8时，p为13时发生明显泄露，此时，频谱发生混叠。
(2)观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。

2，衰减正弦序列

（1） 实验程序
n=0:1:15;
a=0.1;
f1=0.0625;f2=0.04375;f3=0.05625;
xb1=exp(-a*n).*sin(2*pi*f1*n);

figure
subplot(3,2,1)
stem(n, xb1,'.');

title('f=0.0625的时域特性')
xlabel('n')
ylabel('xb1(n)')
grid on

[H, w] = freqz(xb1, 1, [], 'whole', 1);
Hamplitude = abs(H);
subplot(3,2,2)
plot(w, Hamplitude)
title('f=0.0625的幅频响应')

xlabel('w/(2*pi)')

ylabel('|H(exp(jw))|')
grid on

xb2=exp(-a*n).*sin(2*pi*f2*n);
subplot(3,2,3)
stem(n, xb2,'.');

title('f=0.04375的时域特性')
xlabel('n')
ylabel('xb2(n)')
grid on

[H, w] = freqz(xb2, 1, [], 'whole', 1);
Hamplitude = abs(H);
subplot(3,2,4)
plot(w, Hamplitude)

title('f=0.04375的幅频响应')
xlabel('w/(2*pi)')
ylabel('|H(exp(jw))|')
grid on

xb3=exp(-a*n).*sin(2*pi*f3*n);
subplot(3,2,5)
stem(n, xb3,'.');
title('f=0.05625的时域特性')

xlabel('n')
ylabel('xb3(n)')
grid on

[H, w] = freqz(xb3, 1, [], 'whole', 1);
Hamplitude = abs(H);
subplot(3,2,6)
plot(w, Hamplitude)

title('f=0.05625的幅频响应')
xlabel('w/(2*pi)')
ylabel('|H(exp(jw))|')
grid on

（2）实验结果图形显示

（3） 思考题回答
衰减正弦序列可以作为带通滤波器，通过改变f 的值，可以控制其中心频率。
当f 改变时，频谱的形状和谱峰出现的位置会发生混叠和泄漏现象，因为f 改变会导致采样间隔增大，从而导致频谱混叠的产生。可以通过减小采样间隔或增大衰减来使混叠引起的误差在允许的范围内。
3，三角序列与反三角序列
clc
clear all;
close all;

n1=0:1:3;
xc1=n1+1;
n2=4:7;
xc2=8-n2;
xc=[xc1,xc2];
n = [n1,n2];
figure
stem(n,xc);
xlabel('n');
ylabel('xc');
title('三角序列');

n1=0:1:3;
xd1=4-n1;
n2=4:7;
xd2=n2-3;
xd=[xd1,xd2];
n = [n1,n2];
figure
stem(n,xd);
xlabel('n');
ylabel('xd');
title('反三角序列');

N = 8; %N点的FFT
[H1,w1] = freqz(xc,1, 256, 'whole', 1);
Hamplitude1 = abs(H1);
figure
plot(2*w1, Hamplitude1)
title('xc幅频响应')
xlabel('w/pi')
ylabel('|H(exp(jw))|')
grid on
[H2,w2] = freqz(xd,1, 256, 'whole', 1);
Hamplitude2 = abs(H2);
figure
plot(2*w2, Hamplitude2)
title('xd幅频响应')
xlabel('w/pi')
ylabel('|H(exp(jw))|')
grid on
[H3, w3] = freqz(xc, 1, N, 'whole', 1);
Hamplitude3 = abs(H3);
figure
subplot(2, 1, 1)
h3 = stem(2*w3, Hamplitude3, '*');
title('xc幅频响应进行N点FFT')
xlabel('n')
ylabel('|H(exp(jw))|')
grid on
[H4, w4] = freqz(xd, 1, N, 'whole', 1);
Hamplitude4 = abs(H4);
subplot(2, 1, 2)
h4 = stem(2*w4, Hamplitude4, '*');
title('xd幅频响应进行N点FFT')
xlabel('n')
ylabel('|H(exp(jw))|')
grid on
（1） N=8时的的结果图形显示

（2） N=16时的的结果图形显示

(3) 思考题回答
三角波与反三角波序列在时域上互补，8点的FFT变换相同，补0后，两信号的幅频特性没有变化，但FFT变换不同了。这些变换说明，信号的幅频特性不随其补零延拓而变化，但是FFT变换的结果依赖点数N的取值。
4,连续信号采样
clc
clear all;
close all;
N = 128;
%N=16；
f1 = 1/16;
n = 0:N-1;
xn = sin(2*pi*0.125.*n)+ cos(2*pi*(0.125+f1).*n);

figure
stem(n,xn);
figure
subplot(2,1,1),plot(n,abs(fft(xn)));
title('f =1/16 幅频响应');
f2 = 1/64;

xn = sin(2*pi*0.125.*n)+ cos(2*pi*(0.125+f2).*n);
subplot(2,1,2),plot(n,abs(fft(xn)));
title('f =1/64 幅频响应')；
（1)当N=16时的结果图形显示

（2）当N=128时的结果图形显示

（3） 思考题回答
上述结果说明频谱混叠与连续时间信号的时域采样频率fs有关，频谱泄露与时域加窗截断的长度N有关。
当N变化而Δf不变时，频谱图像变得密集，而且幅频响应变得集中而且尖锐了，宽度也变小了。
5，用FFT分别实现xa(n),和 xb(n)的16点圆周卷积和线性卷积。
（1） 试验程序；
%线性卷积和圆周卷积相等的条件：N>=n1+n2

N=16;
n=0:1:15;
p=8;
q=2;
a=0.1;
f=0.0625;
xa=exp(-((n-p).^2)./q);
xb=exp(-a*n).*sin(2*pi*f*n);

%1 线性卷积
x=conv(xa,xb);
XDft= fft(x, 32);
XDftR = abs(XDft);
XDftPhase = angle(XDft);
XDftPhase = unwrap(XDftPhase);

figure(1)
%stem(n, x,'.');
stem(x,'.')
title('x(n)序列')
xlabel('n')
ylabel('x(n)')
grid on
figure(2)
subplot(2, 1, 1)
stem(XDftR, '.');
title('X(k)的幅度')
xlabel('k')
ylabel('|X(k)|')
grid on
subplot(2, 1, 2)
stem(XDftPhase, '.');
title('X(k)的相角')
xlabel('k')
ylabel('fai((X(k)))')
grid on

%2 圆周卷积？

XDft161 = fft(xa, N);
XDft16R1 = abs(XDft161);
XDft16Phase1 = angle(XDft161);
XDft16Phase1 = unwrap(XDft16Phase1);

XDft162 = fft(xb, N);
XDft16R2 = abs(XDft162);
XDft16Phase2 = angle(XDft162);
XDft16Phase2 = unwrap(XDft16Phase2);

XDft16=XDft161.*XDft162;
XDft16R=XDft16R1.*XDft16R2;
XDft16Phase=XDft16Phase2 +XDft16Phase1 ;

x = ifft(XDft16, N);

figure(3)
%stem(n, x,'.');

stem(x,'.')
title('x(n)序列')
xlabel('n')
ylabel('x(n)')
grid on

figure(4)
subplot(2, 1, 1)
t = 0 : 1 : N - 1;
stem(t, XDft16R, '.');

title('X(k)的幅度')
xlabel('k')
ylabel('|X(k)|')
grid on

subplot(2, 1, 2)
stem(t,XDft16Phase, '.');
title('X(k)的相角')
xlabel('k')
ylabel('fai((X(k)))')
grid on
（2）实验结果图像显示。

(3) 结果分析与问题回答
N=16时，xa(n)和 xb(n)的线性卷积与圆卷积结果不相同；当N>=32时，xa(n)和 xb(n)的线性卷积与圆卷积的结果相同
(6)、产生一512点的随机序列xe(n)，并用xc(n)和xe(n)作线性卷积，观察卷积前后xe(n)频谱的变化。要求将xe(n)分成8段，分别采用重叠相加法和重叠保留法。
6、分别采用重叠相加法和重叠保留法对xc(n)和xe(n)作线性卷积。

（1） 试验程序
clc
clear all;
close all;

xe=rand(1,512);
n1=0:1:3;
xc1=n1+1;
n2=4:7;
xc2=8-n2;
xc=[xc1,xc2];
%重叠相加法
yn=zeros(1,519);
for j=0:7
xj=xe(64*j+1:64*(j+1));
xak=fft(xj,71);
xck=fft(xc,71);
yn1=ifft(xak.*xck); %每段的卷积结果
temp=zeros(1,519);
temp(64*j+1:64*j+71)=yn1;
yn=yn+temp;%将每段的卷积结果加到yn中
end;
n=0:518;
figure(1)
subplot(2,1,1);
plot(n,yn);
xlabel('n');ylabel('y(n)');
title('xc(n)与xe(n)的线性卷积的时域波形—重叠相加法');
subplot(2,1,2);
plot(n,abs(fft(yn)));
xlabel('k');ylabel('Y(k)');
axis([0,600,0,300]);
title('xc(n)与xe(n)的线性卷积的幅频特性—重叠相加法');
%重叠保留法¨
k=1:7;
xe1=k-k;
xe_1=[xe1,xe];
yn_1=zeros(1,519);
for j=0:7
xj_1=xe_1(64*j+1:64*j+71);
xak_1=fft(xj_1);
xck_1=fft(xc,71);
yn1_1=ifft(xak_1.*xck_1);% 每段的卷积结果
temp_1=zeros(1,519);
temp_1(64*j+1:64*j+64)=yn1_1(8:71);
yn_1=yn_1+temp_1;% 将每段的卷积结果加到yn中
end;

n=0:518;
figure(2)
subplot(2,1,1);
plot(n,yn_1);
xlabel('n');ylabel('y(n)');
title(' xc(n)与xe(n)的线性卷积的时域波形—重叠保留法');
subplot(2,1,2);
plot(n,abs(fft(yn_1)));
xlabel('k');ylabel('Y(k)');
axis([0,600,0,300]);
title('xc(n)与xe(n)的线性卷积的幅频特性—重叠保留法');

(2)实验结果图形显示。

(3) 结果分析
重叠相加法和重叠保留法都是将长序列分段成短序列。不同的是重叠相加法是将分成的短序列分别与短序列计算线性卷积，最后将这些分段卷积结果进行重叠相加；而重叠保留法是将分成的短序列与短序列计算循环卷积，再从个循环卷积的结果中提取相当于线性卷积的部分。

实验三 用双线性变换法设计
IIR数字滤波器
一、 实验目的
(1) 熟悉用双线性变换法设计IIR数字滤波器的原理与方法。
(2) 掌握数字滤波器的计算机仿真方法。
(3) 通过观察对实际心电图信号的滤波作用， 获得数字滤波的感性知识。
二、 实验内容
(1) 用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。
(2) 以 0.02π为采样间隔， 打印出数字滤波器在频率区间［0, π/2］上的幅频响应特性曲线。
(3) 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图， 观察总结滤波作用与效果。
三、 实验步骤
(1) 复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容， 按照
例 6.4.2， 用双线性变换法设计数字滤波器系统函数H(z)。
例 6.4.2 中已求出满足本实验要求的数字滤波器系统函数：

A=0.090 36
B1=1.2686,C1=-0.7051
B2=1.0106,C2=-0.3583
B3=0.9044,C3=-0.2155

由(10.5.1)式和(10.5.2)式可见，滤波器H(z)由三个二阶滤波器H1(z), H2(z)和H3(z)级联组成， 如图10.5.1 所示。

图10.5.1 滤波器H(z)的组成

(2) 编写滤波器仿真程序, 计算H(z)对心电图信号采样序列x(n)的响应序列y(n)。设yk(n)为第k级二阶滤波器Hk(z)的输出序列， y k-1 (n)为输入序列， 如图10.5.1 所示。 由(10.5.2)式可得到差分方程：
yk(n)=Ay k-1 (n)+2Ay k-1 (n-1)+Ay k-1
(n-2)+Bkyk(n-1)+Ckyk(n-2) (10.5.3)
(3) 在通用计算机上运行仿真滤波程序， 并调用通用绘图子程序， 完成实验内容(2)和(3)。
四、思考题
用双线性变换法设计数字滤波器过程中， 变换公式

中T的取值， 对设计结果有无影响? 为什么?
五、实验报告要求
(1) 简述实验目的及原理。
(2) 由所打印的|H(ejω)|特性曲线及设计过程简述双线性变换法的特点。
(3) 对比滤波前后的心电图信号波形， 说明数字滤波器的滤波过程与滤波作用。
(4) 简要回答思考题。
六、心电图信号采样序列x(n)
人体心电图信号在测量过程中往往受到工业高频干扰， 所以必须经过低通滤波处理后， 才能作为判断心脏功能的有用信息。 下面给出一实际心电图信号采样序列样本x(n)， 其中存在高频干扰。 在实验中， 以x(n)作为输入序列， 滤除其中的干扰成分。
七，实验程序与结果图像
1、试验程序
x=[-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];
k=1;
close all;
figure(1)
subplot(2,2,1)
n=0:55;
stem(n,x,'.')
axis([0 55 -100 50]);
hold on;
n=0:60;
m=zeros(61);
plot(n,m);
xlabel('n');
ylabel('x(n)');
title('心电图信号采样序列x(n)');
B=[0.09036 2*0.09036 0.09036];
A=[1.2686 -0.7051];
A1=[1.0106 -0.3583];
A2=[0.9044 -0.2155];
while(k<=3)
y=filter(B,A,x);
x=y;
if k==2
A=A1;
end
if k==3
A=A2;
end
k=k+1;
end
subplot(2,2,3)
n=0:55;
stem(n,y,'.');
axis([0 55 -15 5]);
hold on;
n=0:60;
m=zeros(61);
plot(n,m);
xlabel('n');
ylabel('y(n)');
title('三级滤波会的心电图信号');
A=[0.09036,0.1872,0.09036];
B1=[1,-1.2686,0.7051];
B2=[1,-1.0106,0.3583];
B3=[1,-0.9044,0.2155];
[H1,w]=freqz(A,B1,100);
[H2,w]=freqz(A,B2,100);
[H3,w]=freqz(A,B3,100);
H4=H1.*(H2);
H=H4.*(H3);
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
subplot(2,2,2)
plot(w/pi,db);
axis([0,0.5,-50,10]);
grid on;
title('滤波器的幅频响应曲线');
实验绘图如下：

实验总结：
1、双线性变换法的特点是：
双线性变换比脉冲响应法的设计计算直接和简单。由于s与z之间的简单代数关系，所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。双线性换法的主要优点是S平面与Z平面一一单值对应，双线性变换法的缺点：Ω与ω的非线性关系，导致数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变。
2、 数字滤波器的滤波过程和滤波作用简述如下：
数字滤波器的作用是利用离散时间系统的特性对输入信号波形（或频谱）进行加工处理，或者说利用数字方法按预定的要求对信号进行变换。把输入序列x(n)变换成一定的输出序列y(n)从而达到改变信号频谱的目的。
3、 用双线性变换法设计数字滤波器中，变换公式中T的取值与最终的设计结果无关，

故参数T的取值和最终的设计结果无关。