信号运算与离散时间系统表示方法

实验 二 信号运算与离散时间系统表示方法
实验目的：
掌握基于MATLAB环境下信号的各种运算，特别是卷积运算；熟悉离散时间系统模型的各种表示方法及其相互转换；掌握离散时间傅里叶变换、快速傅氏变换和Z变换。
实验环境：
硬件环境:计算机,软件环境：MATLAB平台。
实验内容与步骤：
1 信号的运算
在数字信号处理领域，对信号所做的基本运算主要包括:信号加、信号乘、移位、采样和、采样积、翻转以及卷积等。下面，就对怎样通过MATLAB来实现这些具体的操作进行详细说明。
1.1 信号加
MATLAB实现 x= x1 + x2
注意:x1和x2应该具有相同的长度，位置对应，才能相加。否则，需要先通过zeros函数左右补零后再相加。例如：
n1=1:1:5;
x10=[1 0.7 0.4 0.1 0];
n2=3:8; % n2=3:1:8;
x20=[0.1 0.3 0.5 0.7 0.9 1];
n=1:8;
x1=[x10 zeros(1,8-length(n1))];
x2=[zeros(1,8-length(n2)) x20];
x=x1+x2;
subplot(3,1,1); stem(n,x1);
subplot(3,1,2); stem(n,x2);
subplot(3,1,3); stem(n,x);
1.2 信号延迟
MATLAB实现 y(n)=x(n-k)
1.3 信号乘
MATLAB实现 x=xl.*x2
这是信号的点乘运算，所以同样需要x1和x2的长度要相等这一前提条件。例如:
n1=1:5;
x10=[1 0.7 0.4 0.1 0];
n2=3:8;
x20=[0.1 0.3 0.5 0.7 0.9 1];
x1=[x10 zeros(1,8-length(n1))],
x2=[zeros(1,8-length(n2)) x20],
x=x1.x2
运行结果
x =
0 0 0.0400 0.0300 0 0 0 0
1.4 信号乘以常数，即信号的变化幅度
MATLAB实现 y=kx
1.5 信号的翻转
MATLAB实现 y=fliplr(x)
A = 1:10 %若A为行向量，y为A的翻褶。
A =

 1     2     3     4     5     6     7     8     9    10

y =

10     9     8     7     6     5     4     3     2     1

1.6 信号采样和
MATLAB实现 y=sum(x(ns:ne))
1.7 信号采样积
MATLAB实现 y = prod(x(ns:ne))
1.8 信号能量
MATLAB实现 Ex=sum(abs(x).^2)
1.9 信号功率
MATLAB实现 Px=sum((abs(x)).^2)/N
1.10 卷积
MATLAB实现 y=conv(x, h)
例如：
x = [1 3 5 7 6 4 2];
h = [ 3 5 9 2];
y=conv(x,h)
运行结果：
y =
3 14 39 75 104 115 94 58 26 4
上面后来介绍的这几种信号计算，应用起来相当简单直接。同学们可以直接根据MATLAB实现方法来应用。
2 离散时间系统
MATLAB内部，为一个线性离散系统模型提供了很多种表示方法，分别可以用在不同的场合。MATLAB 提供了6种表达线性离散系统模型的方法。
• 传递函数法
• 零极点增益法
• 状态空间法
• 部分分式法
• 二次分式法
• Lattice结构法

2.1 传输函数
传输函数是主要的离散时间系统模型，它是离散系统的 z 域表示形式，为两个多项式值。系统函数的表达式为

在这里，bi 和 ai 是离散系统的系数，分别存储在两个向量中b、a。
N = max(n, m)−1为滤波器阶数。
2.2 零极点增益法
零极点增益法，其实也是系统函数的一种方法，它不过是把分子和分母上的多项式表示为一系列一阶因式相乘的形式，其表达式如下

其中k 为系统的增益，zi 和 pi 是离散系统的零极点，分别存储在两个向量中z、p。

2.3 状态空间法 略
2.4 部分分式法
任何一个系统函数表示的系统，都存在一个带余数多项式的部分分式与之对应。带余数多项式的部分分式表达式如下所示

此式成立的条件是H(z)没有重复的极点，n 是系统的阶数。当存在 sr 阶重复的极点时，H(z)则表示为

所以，使用一个极点列向量 p =[p(1) , p(2) , …, p(n)]，一个与极点对应的分子列向量 r =[r(1) , r(2) , …, r(n)]和余数多项式系数行向量k =[k(1) , k(2) , …, k(m―n+1)]，就可以完整地标识整个系统了。
2.5 二次分式法
任何系统函数H(z)都可以用如下形式的二次分式来表示
→
这种表达式实际上是系统的二阶子系统级联实现，其中L是描述系统的二次分式的数目，Hk(z)为各阶子系统。
MATLAB用一个L×6的矩阵sos来表示离散系统的二次分式形式，sos矩阵的每一行代表一个二次分式，前半部分为分子系数，后半部分为分母系数。矩阵的形势如下所示

2.6 Lattice结构法

3 离散时间系统模型变换
MATLAB的工具箱，为同一离散系统的不同系统模型间的多种表达方式的变换，提供了需多功能丰富而快捷的MATLAB函数，如表2-1所示。
表 2-1 离散系统模型变换MATLAB函数
系统原模型 系统转换模型 函 数
传输函数 状态空间 tf2ss
传输函数 零极点增益 tf2zp
零极点增益 状态空间 zp2ss
零极点增益 传输函数 zp2tf
零极点增益 二次分式 zp2sos
状态空间 传输函数 ss2tf
状态空间 零极点增益 ss2zp
状态空间 二次分式 ss2sos
二次分式 传输函数 sos2tf
二次分式 零极点增益 sos2zp
二次分式 状态空间 sos2ss
多项式 Lattice结构 poly2rc
Lattice结构 多项式 rc2poly

3.1 tf2ss：系统函数模型转换为状态空间模型

3.2 tf2zp：传输函数模型转换为零极点增益模型
调用方式
[z，p，k]=tf2zp(num, den)：输入向量num和向量den分别为系统传输函数模型的分子和分母所组成的向量。需要注意的是：系统传输函数模型的分子多项式和分母多项式的长度必须相等，否则需要进行补零。
输出向量z，p和k为系统的零极点增益模型的零点向量z、极点向量p和系统增益k。
应用说明
例：将系统

变换成零极点增益模型。
num=[1 3 5 0]; %进行补零;
den=[1 8 1 3];
[z,p,k]=tf2zp(num, den)
运行结果
z =

0.0000 + 0.0000i
-1.5000 + 1.6583i
-1.5000 - 1.6583i

p =

-7.9216 + 0.0000i
-0.0392 + 0.6141i
-0.0392 - 0.6141i

k =
1

；

3.3 zp2ss：零极点增益模型转换为状态空间模型

3.4 zp2tf零极点增益模型转换为传输函数模型
调用方式
[num，den] = zp2tf(z, p, k)：将系统零极点增益模型的零点向量z、极点向量p和系统增益k，转换为传输函数模型[num, den]。
应用说明

z =[-1.5000+1.6583i; -1.5000-1.6583i];
p =[-7.4641; -0.5359];
k = 1;
[num,den]=zp2tf(z,p,k)
运行结果
num =
1.0000 3.0000 5.0000
den =
1.0000 8.0000 4.0000

3.5 zp2sos：零极点增益模型转换为二次分式模型

调用方式
(1) [sos, g] = zp2sos(z，p，k)：将系统零极点增益模型的零点向量z、极点向量p和系统增益k，转换为二次分式模型。g为整个系统的增益，即H(z)=g*H1(z)H2(z)…*HL(z)。
注意:零、极点必须以共扼复数对的形式出现。
(2) [sos, g] = zp2sos(z，p，k，dir_flag)：参数dir_flag指定sos中行的顺序，即
• up，首行中所包含的极点离原点最近，离单位圆最远；
• down，首行中所包含的极点离原点最远，离单位圆最近。
参数dir_flag的默认值为up。
(3) [sos，g] = zp2sos(z，p，k，dir_flag，scale)：参数scale指定所需要的增益因子，其取值为none、Inf 或 2；默认值为none。
应用说明
[z, p, k] = butter(5, 0.4);
[sos1, g1] = zp2sos(z, p, k)
[sos2, g2] = zp2sos(z, p, k, ‘down’)
[sos3, g3] = zp2sos(z, p, k, ‘down’, 2)
运行结果
sos1 =
1.0000 1.0000 0 1.0000 -0.1584 0
1.0000 2.0000 1.0000 1.0000 -0.3493 0.1303
1.0000 2.0000 1.0000 1.0000 -0.4777 0.5457
g1 =
0.0219
sos2 =
1.0000 2.0000 1.0000 1.0000 -0.4777 0.5457
1.0000 2.0000 1.0000 1.0000 -0.3493 0.1303
1.0000 1.0000 0 1.0000 -0.1584 0
g2 =
0.0219
sos3 =
0.2887 0.5775 0.2887 1.0000 -0.4777 0.5457
0.2781 0.5562 0.2781 1.0000 -0.3493 0.1303
0.3429 0.3429 0 1.0000 -0.1584 0
g3 =
0.7970
3.6 ss2tf：状态空间模型转换为系统函数模型
3.7 ss2zp：状态空间模型转换为零极点增益模型
3.8 ss2sos：状态空间模型转换为二次分式模型

3.9 sos2tf：二次分式模型转换为传输函数模型

调用方式
[num, den] = sos2tf(sos)：将系统的二次分式模型sos转换为系统的传输函数模型[num，den]。
应用说明
sos =[1 3 2 -2 0 1; -2 1 1 4 6 3; -3 2 1 0 5 4; 4 1 2 -2 -1 -3; 1 -1 1 -1 -1 2];
[num, den] = sos2tf(sos)
运行结果
num =
24 26 -51 5 17 -42 -3 -5 9 16 4
den =
0 -80 -304 -392 -170 292 635 406 -93 -222 -72

3.10 sos2zp：二次分式模型转换为零极点增益模型
调用方式
[z, p, k] = sos2zp(sos, g)：将系统的二次分式模型sos转换为系统的零极点增益模型[z, p, k]，输入参数g为二次分式模型系统的增益因子，默认值为1。
应用说明
sos=[1 3 2 -2 0.5 1; -2 1 1 4 6 3; -3 -2 1 0.4 5 4; 4 1 2 -2 -1 -3];
[z, p, k] = sos2zp(sos,2)
运行结果
z =
-2.0000
-1.0000
1.0000
-0.5000
-1.0000
0.3333
-0.1250 + 0.6960i
-0.1250 - 0.6960i
p =
0.8431
-0.5931
-0.7500 + 0.4330i
-0.7500 - 0.4330i
-11.6410
-0.8590
-0.2500 + 1.1990i
-0.2500 - 1.1990i
k =
7.5000
3.11 sos2ss：二次分式模型转换为状态空间模型

3. 12 poly2rc:由多项式计算Lattice网状结构的映射系数

3.13 poly2rc：由Lattice网状结构的映射系数计算多项式

4 DFT变换、FFT 变换与Z变换
4.1 离散傅里叶变换
若将DFT变换的定义写成矩阵形式，则得到
X =A•x
x为信号。
Dftmtx函数：用来计算DFT变换矩阵A的函数。
调用方式
(1) A = dftmtx (n)：返回n×n的DFT变换矩阵A。若x为给定长度的行向量，则y=x*A，返回x的DFT变换y。
(2) iA =conj(dftmtx(n))/n：返回n×n的IDFT变换矩阵iA。
应用说明
例1：
A = dftmtx(4);
iA = conj(dftmtx (4))/4;
运行结果
A =
1.0000 1.0000 1.0000 1.0000
1.0000 0 - 1.0000i -1.0000 0 + 1.0000i
1.0000 -1.0000 1.0000 -1.0000
1.0000 0 + 1.0000i -1.0000 0 - 1.0000i

iA =
0.2500 0.2500 0.2500 0.2500
0.2500 0 + 0.2500i -0.2500 0 - 0.2500i
0.2500 -0.2500 0.2500 -0.2500
0.2500 0 - 0.2500i -0.2500 0 + 0.2500i
例2：如果x(n) = sin(nπ/8) + sin(nπ/4)是一个N=16的有限序列，用MATLAB求其DFT的结果，并画出其结果图，如图2-3所示。
程序
N=16;
n=0:1:N-1; %时域采样
xn=sin(npi/8)+sin(npi/4);
k=0:1:N-1; %频域采样
WN=exp(-j2pi/N);
nk=n’k;
WNnk=WN.^nk;
Xk=xnWNnk;
subplot(2, 1, 1)
stem(n,xn);
subplot(2, 1, 2)
stem(k, abs(Xk));
运算结果
Xk =
-0.0000 -0.0000 - 8.0000i -0.0000 - 8.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 8.0000i 0.0000 + 8.0000i
也可以编写出有限长序列的离散傅里叶逆变换MATLAB程序。这里给出计算离散傅里叶逆变换的MATLAB函数文件。

n=0:1:N-1;
k=0:1;N-1;
WN=exp(-j2pi/N);
nk=n’k;
WNnk=WN.^(-nk);
xn=(XkWNnk)./N;
subplot(2,1,1)
stem(k,abs(Xk));
subplot(2,1,2)
stem(n,real(xn));

例：假设现有含有3种频率成分，f1=20 Hz，f2=20.5Hz，f3=40Hz，采样频率为fs=100Hz，x(n)=sin(2nπ f1/ fs)+sin(2nπ f2/ fs)+sin(2nπ f3/ fs)
(1) 求在记录中最少的点数;
(2) 求x(n)在0～128之间的DFT的X(k)；
(3) 求把(2)中的x(n)以补零方式使其加长到0～512之间后的DFT的X(k)；
(4) 求x(n)在0～512之间的DFT的X(k)。
解:
(1) 要想分辨出信号的所有频谱成分，则信号采样的频谱分辨率应为
F =0.5 Hz
所以记录中的最少点数为
N = 2fh /F = 2.40/0.5=160
(2)、(3)和(4)的解决程序和运行结果如下所示。
程序
% 第2问题的程序
N=128;
fs=100;
n=0:1:N-1;
f1=20;
f2=20.5;
f3=40;
xn=sin(2pif1n/fs)+sin(2pif2n/fs)+sin(2pif3n/fs);
Xk=fft(xn);
AXk=abs(Xk(1:N/2));
figure(1)
subplot(2,1,1)
plot(n,xn);
k=(0:N/2-1)fs/N;
subplot(2,1,2)
plot(k,AXk);
% 第3问题的程序
M=512;
xn=[xn zeros(1,M-N)];
Xk=fft(xn);
AXk=abs(Xk(1:M/2));
m=0:1:M-1;
figure(2)
subplot(2,1,1)
plot(m,xn);
k=(0:M/2-1)fs/M;
subplot(2,1,2)
plot(k,AXk);
% 第4问题的程序
n=0:1:M-1;
xn=sin(2pif1n/fs)+sin(2pif2n/fs)+sin(2pif3n/fs);
Xk=fft(xn);
AXk=abs(Xk(1:M/2));
figure(3)
subplot(2,1,1)
plot(n,xn);
k=(0:M/2-1)*fs/M;
subplot(2,1,2)
plot(k,AXk);

上面程序执行后的结果分别如图2-4、图2-5和图2-6所示。
图2-4中N为不满足取样点数公式N≧2fh/F，中所有的频率成分。图2-5为以补零方式便其加长到0～512之间后序列的DFT结果图，从图中可以看出，补零对分辨率没有影响，但是对频谱起到了平滑的作用。图2-6为取样512点的序列的DFT结果图，由于N满足采样点数的公式，所以结果图中能区分序列中所有的频率成分。
4.2 快速傅里叶变换(FFT)
在信号处理中，DFT的计算具有举足轻重的地位，信号的相关、滤波、谱估计等都要通过DFT来实现。然而，当N很大的时候，求一个N点的DFT要完成N×N次复数乘法和N(N-l)次复数加法，其计算量相当大。1965年J. W. Cooley和J.W. Tukey巧妙地利用WN因子的周期性和对称性，构造了一个DFT快速算法，即快速傅里叶变换(FFT)。
MATLAB为计算数据的离散快速傅里叶变换，提供了一系列丰富的数学函数，主要有fft、ifft，fft2、ifft2，fftn、ifftn和fftshift、iffshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能便所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。
fft和ifft函数
调用方式
(l) Y = fft(X)
参数说明
· 如果X是向量，则采用傅里叶变换来求解X的离散傅里叶变换；
· 如果X是矩阵，则计算该矩阵每一列的离散傅里叶变换；
· 如果X是(ND)维数组，则是对第一个非单元素的维进行离散傅里叶变换。
(2) Y = fft(X, N)
参数说明
N是进行离散傅里叶变换的x的数据长度，可以通过对x进行补零或截取来实现。
(3) Y = fft(X, [], dim) 或Y = fft(X, N, dim)
参数说明
· 在参数dim指定的维上进行离散傅里叶变换；
· 当X为矩阵时，dim用来指定变换的实施方向，dim=l，表明变换按列进行；dim=2表明变换按行进行。
函数ifft的参数应用与函数fft完全相同。
应用说明
例1：fft的应用
X = [2 1 2 8];
Y = fft(X);
运行结果
Y =
13.0000 0 + 7.0000i -5.0000 0 - 7.0000i
例2；fft(x, N, dim)的应用
A = [2 5 7 8;
1 4 0 5;
3 8 5 1;
9 1 2 7];
Y = fft(A, [], 1), %对矩阵A按列进行变换，变换的数据长度为默认值，即其长度。
运行结果
Y =
15.0000 18.0000 14.0000 21.0000
-1.0000 + 8.0000i -3.0000 - 3.0000i 2.0000 + 2.0000i 7.0000 + 2.0000i
-5.0000 8.0000 10.0000 -3.0000
-1.0000 - 8.0000i -3.0000 + 3.0000i 2.0000 - 2.0000i 7.0000 - 2.0000i
Z = fft (A, 3, 2) %对矩阵A按行进行变换，变换的数据长度N=3。
Z =
14.0000 -4.0000 + 1.7321i -4.0000 - 1.7321i
5.0000 -1.0000 - 3.4641i -1.0000 + 3.4641i
16.0000 -3.5000 - 2.5981i -3.5000 + 2.5981i
12.0000 7.5000 + 0.8660i 7.5000 - 0.8660i
例3：fft在信号分析中的应用
使用频谱分析方法从受噪声污染的信号x(t)中鉴别出有用的信号。
程序
t = 0:0.001:1; %采样周期为0.001s，即采样频率为1000Hz;
%产生受噪声污染的正弦波信号;
x=sin(2pi100t)+sin(2pi200t)+rand(size(t));
subplot(2,1,1)
plot(x(1:50)); %画出时域内的信号;
Y=fft(x, 512); %对x进行512点的傅里叶变换;
f=1000(0:256)/512; %设置频率轴(横轴)坐标，1000为采样频率
subplot(2,1,2)
plot(f,Y(1:257)); %画出频域内的信号;

由上面的图2-7可以看出，从受噪声污染信号的时域形式中，很难看出正弦波的成分。但是通过对x(t)作傅里叶变换，把时域信号变换到频域进行分析，可以明显看出信号中100Hz和200Hz的两个频率分量。
fft2和iff2函数
fftshift和ifftshift函数
调用方式
Z = fftshift(Y)
此函数可用于将傅里叶变换结果Y(频域数据)中的直流成分(即频率为0处得值)移到频谱的中间位置。
4.3 Z变换
在MATLAB函数中，可直接利用Residuez函数进行Z变换。
例：求解

程序
b =1;
a = poly([-0.2 0.5 0.5 -0.6 -0.6]);
[z0, p, k]=residuez(b, a)
运行结果
z0 =
0.1826
0.4463
0.2032
0.1476
0.0204
p =
-0.6000
-0.6000
0.5000
0.5000
-0.2000
k =
[]
所以，X(z)的部分分式表示为

则相对应的Z反变换为

思考题：
1．分别产生一个周期为1的标准方波、锯齿波以及标准三角波？
2．分别画出典型sinc函数和Dirichlet函数的波形图？
3．首先产生两个信号A与信号B，然后求A的信号翻转、信号延迟、信号能量和信号功率，以及求解信号A与信号B的信号加、信号乘和信号卷积等运算。
4．通过MATLAB的help文档掌握impz求Z反变换的方法。