在离散时间信号与系统的分析与处理中,经常需要对信号进行操作,也就是对序列进行运算,常用的序列运算有序列的加与乘、移位与翻转、尺度变换、累加、差分、卷积和以及相关等运算,所有的这些序列运算都可以利用MATLAB加以实现。
序列的加法与乘法是指两个或两个以上的序列对应序号的样点值相加或相乘的运算。 y 1 ( n )= x 1 ( n )+ x 2 ( n ), y 1 ( n )就是两个序列 x 1 ( n )、 x 2 ( n )对应项相加后得到的和序列。 y 2 ( n )= x 1 ( n )· x 2 ( n ), y 2 ( n )就是两个序列 x 1 ( n )、 x 2 ( n )对应项相乘后得到的积序列。
例1—1 利用MATLAB实现下面两个序列的相加运算。
x 1 ( n ) = [1, 0.7, 0.4, 0.1, 0, 0.1], n= [1, 2, 3, 4, 5, 6],
x 2 ( n ) = [0.2, 0.1, 0.3, 0.5, 0.7, 0.9, 1], n= [2, 3, 4, 5, 6, 7, 8],
求 x ( n )= x 1 ( n )+ x 2 ( n )。
解: MATLAB实现程序如下:
n1=1:6; x10=[1 0.7 0.4 0.1 0 0.1]; n2=2:8; x20=[0.2 0.1 0.3 0.5 0.7 0.9 1]; n=1:8; x1=[x10 zeros(1,8-length(n1))];%对x10进行右侧补零 x2=[zeros(1,8-length(n2))x20];%对x20进行左侧补零 x=x1+x2 subplot(3,1,1);stem(n,x1,'.k');subplot(3,1,2);stem(n,x2,'.k');subplot(3,1,3);stem(n,x, '.k');
计算结果为
x= 1.0000 0.9000 0.5000 0.4000 0.5000 0.8000 0.9000 1.0000
结果如图1—11所示。
图1—11 序列的加法
例1—2 利用MATLAB实现下面两个正弦序列的相加运算。
x 1 ( n ) = sin(0.06π n ), x 2 ( n ) = sin(0.24π n )
求 y ( n )= x 1 ( n )+ x 2 ( n )。
解: MATLAB实现程序如下:
t=0:0.001:0.6; x1= sin(2 * pi * 30 * t); figure; subplot(3,1,1); plot(x1); x2=sin(2 * pi * 120 * t); subplot(3,1,2); plot(x2); y=x1+x2; subplot(3,1,3); plot(y);
结果如图1—12所示。
图1—12 正弦序列的加法
例1—3 利用MATLAB实现下面两个序列函数的相乘运算。
x ( n ) = e n , y ( n ) = sin(10π n ),
求 z ( n )= x ( n )· y ( n )。
解: MATLAB实现程序如下:
clear n=[-10:1:10]; x=exp(n); y=sin(2 * pi * 5 * n); z=x. * y; subplot(3,1,1);stem(n,x,'.k');subplot(3,1,2);stem(n,y,'.k');subplot(3,1,3);stem(n,z,'.k')
结果如图1—13所示。
图1—13 序列的乘法
例1—4 利用MATLAB实现数字0和数字2的语音信号的相加与相乘运算。
解: 首先给出数字0和数字2的语音信号(见图1—4和图1—15)。
[au0,fs,bits]=wavread('au0.wav'); [au2,fs,bits]=wavread('au2.wav'); l1=size(au0); l2=size(au2); l=max(l1,l2); au0=[au0'zeros(1,l-l1)]'; figure plot(au0); xlabel('n') ylabel('au0') figure plot(au2); xlabel('n') ylabel('au2');
图1—14 数字0的语音信号序列图
图1—15 数字2的语音信号序列图
语音信号的相加运算(见图1—16):
au02=au0+au2; figure plot(au02); xlabel('n'); ylabel('au02'); soundsc(au02,fs);
图1—16 数字0和数字2的语音信号相加后的信号序列图
语音信号的相乘运算(见图1—17):
au_02=au0. * au2; figure plot(au_02); xlabel('n'); ylabel('au_02'); soundsc(au_02,fs);
图1—17 数字0和数字2的语音信号相乘后的信号序列图
序列的移位也称作序列的延时,在数学上表示为
当 m >0时, x ( n - m )是序列 x ( n )右移 m 位后的序列; x ( n + m )是序列 x ( n )左移 m 位后的序列。
例1—5 设序列 x ( n )=[1, 1, 2, 2, 4, 4, 5, 4, 2, 2, 1], n =[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]。利用MATLAB编程实现 x ( n )的移位 x ( n +2)和 x ( n -2)。
解: MATLAB实现程序如下:
n=[0:10]; x=[1,1,2,2,4,4,5,4,2,2,1]; subplot(3,1,1);stem(n,x,'.k');axis([-3,13,0,5]) n1=n-2; subplot(3,1,2);stem(n1,x,'.k');axis([-3,13,0,5]) n2=n+2; subplot(3,1,3);stem(n2,x,'.k');axis([-3,13,0,5])
结果如图1—18所示。
图1—18 序列的移位
例1—6 实现“0”~“9”数字语音信号的移位相加运算。
解: MATLAB实现程序如下:
[x,fs,bits]=wavread('1to9.wav'); z1=[zeros(1000,1);x;zeros(4000,1)]; z2=[zeros(2000,1);x;zeros(3000,1)]; z3=[zeros(3000,1);x;zeros(2000,1)]; z4=[zeros(4000,1);x;zeros(1000,1)]; z5=[zeros(5000,1);x]; x=[x;zeros(5000,1)]; y=z1+z2+z3+z4+z5+x; soundsc(y,fs);
结果如图1—19所示。
图1—19 原始语音信号与移位相加后的语音信号
序列的翻转是以 n =0的纵轴为对称轴,将序列左右两边加以对调,即
序列 y ( n )就是 x ( n )的翻转序列。
例1—7 已知函数 x ( n )=e n -2 ,求 x ( n )的翻转序列。
解: MATLAB实现程序如下:
n=[-4:1:10]; x=exp(n-2); subplot(2,1,1);stem(n,x,'.k');axis([-10,10,0,10]) n=fliplr(-n); x=fliplr(x); subplot(2,1,2);stem(n,x,'.k');axis([-10,10,0,10])
结果如图1—20所示。
图1—20 序列的翻转
序列的尺度变换包括幅度尺度变换和时间尺度变换。
1)幅度尺度变换
表示将序列 x ( n )的序列值按比例放大 m 倍。体现为序列按一定比例伸长或缩短。
例1—8 设 x ( n )=sin(2π n ),求利用MATLAB实现 x ( n )的幅度尺度变换。
解: MATLAB实现程序如下:
n=[-1:0.02:1]; x=sin(2 * pi * n); subplot(2,1,1);stem(n,x,'.k'); y=2 * x; %对信号进行幅度加倍 subplot(2,1,2);stem(n,y,'.k');
结果如图1—21所示。
图1—21 序列的幅度尺度加倍
2)时间尺度变换
时间尺度变换表示为序列按一定比例伸长或缩短。
表示序列每隔 m 点取一点形成的抽样序列,相当于将时间轴压缩为原来的 倍,即序列的抽取。
表示序列每两点之间插入 m -1点的序列值形成的插值序列,相当于将时间轴扩展了 m 倍,即序列的插值。
例1—9 实现信号的抽取运算,其中 x ( n )=sin(2π50 n ), y ( n )=sin(4π50 n )。
解: MATLAB实现程序如下:
n=0:0.0007:2; x=sin(2 * pi * 50 * n); y=decimate(x,2);%对信号进行2倍抽取 subplot(2,1,1);stem(x(1:250),'.k');subplot(2,1,2);stem(y(1:140),'.k ');
运行结果如图1—22所示。
例1—10 实现信号的插值运算,其中 x ( n )=sin(2π50 n ), y ( n )=sin(π50 n )。
解: MATLAB实现程序如下:
n=0:0.0016:2; x=sin(2 * pi * 50 * n); y=interp(x,2);%对信号进行2倍插值 subplot(2,1,1);stem(x(1:140),'.k');subplot(2,1,2);stem(y(1:250),'.k ');
图1—22 序列的抽取
运行结果如图1—23所示。
图1—23 序列的插值
例1—11 对数字0的语音信号进行幅度尺度与时间尺度的变换。
解: 对语音信号进行时间尺度变换就是相当于改变数字语音信号的采样频率,其MATLAB实现程序如下:
[au0,fs,bits]=wavread('au0.wav'); y1=5 * au0; y2=decimate(au0,5); soundsc(y1,fs); soundsc(y2,fs); figure subplot(1,2,1) plot(y1);xlabel('n');ylabel('y1');title('幅度尺度变换'); subplot(1,2,2) plot(y2);xlabel('n');ylabel('y2');title('时间尺度变换');
结果如图1—24所示。
图1—24 数字0的语音信号的幅度尺度和时间尺度的变换结果
一个序列 x ( n )的累加序列定义为
表示在某一 n 0 时刻上的序列值 y ( n 0 )等于这一时刻上 x ( n 0 )的值以及 n 0 以前时刻所有序列值的和。
例1—12 实现信号的累加运算,其中 x ( n )=[0.1, 0.2, 0.3, 0.5, 0.5, 0.1, 0.5, 0.6]; n =[1, 2, 3, 4, 5, 6, 7, 8],求 x ( n )的累加和序列。
解: MATLAB实现程序如下:
n=[1:8]; x=[0.1 0.2 0.3 0.5 0.5 0.1 0.5 0.6]; subplot(2,1,1);stem(n,x,'.k'); y=zeros(1,8); n1=1; for n2=1:8 y(n2)=sum(x(n1:n2)) end subplot(2,1,2);stem(n,y,'.k');
运行结果如图1—25所示。
图1—25 序列的累加
例1—13 计算数字0的语音信号的累加和信号。
解: MATLAB实现程序如下:
[au0,fs,bits]=wavread('au0.wav'); N=length(au0); n=[1:N]; for n=1:N y(n)=sum(au0(1:n)); end figure plot(y);xlabel('n') soundsc(y,fs)
运行结果如图1—26所示。
图1—26 数字0的语音信号的累加和信号
差分运算有两种,即向前差分和向后差分。
一阶向前差分:
一阶向后差分:
可见,
▽ x ( n )=Δ x ( n -1)
二阶向前差分:
在连续系统的分析与处理中,系统零状态响应常常利用卷积积分来求解。对于离散系统来说,卷积和是求解离散时间线性移不变系统零状态响应的有效方法。卷积和定义为
其中“∗”为卷积和运算符号,通常这种离散卷积也叫做离散线性卷积。如果序列 x ( n )为有限长序列,序列长度为 L 1 ,序列的取值区间为[ N 1 , N 2 ]。序列 h ( n )为有限长序列,序列长度为 L 2 ,序列的取值区间为[ N 3 , N 4 ],即在卷积和的运算中,
x ( m ): N 1 ≤ m ≤ N 2
h ( n - m ): N 3 ≤ n - m ≤ N 4
则
可见两个序列的卷积和序列 y ( n )也是有限长序列,其取值区间为[ N 1 + N 3 , N 2 + N 4 ],序列 y ( n )的长度为 L 1 + L 2 -1。
下面具体讨论几种离散卷积的运算方法。
1)公式法
公式法就是直接利用卷积和的定义来计算离散卷积的方法。
例1—14 已知两个序列
x ( n )=2 n , 0≤ n ≤2
h ( n )=1, 0≤ n ≤2
试计算这两个序列的卷积和序列 y ( n )。
解: 由卷积和计算公式
2)图解法
图解法就是在序列图上实现离散卷积的计算。通过序列的翻转、移位、相乘以及相加来完成卷积和的计算。
例1—15 将例1—14用图解法实现。
解: 利用MATLAB实现基于图解法的卷积和计算程序如下:
n=-5:5; x=zeros(1,length(n)); x(6)=1;x(7)=2;x(8)=4; h=zeros(1,length(n)); h([find((n > =0)&(n < =2))])=1; subplot(4,2,1);stem(n,x,' * k'); subplot(4,2,3);stem(n,h,'.k'); n1=fliplr(-n);h1=fliplr(h); subplot(4,2,5);stem(n,x,' * k');hold on;stem(n1,h1,'.k'); h2=[0,h1];h2(length(h2))=[];n2=n1; subplot(4,2,7);stem(n,x,' * k');hold on;stem(n2,h2,'.k'); h3=[0,h2];h3(length(h3))=[];n3=n2; subplot(4,2,2);stem(n,x,' * k');hold on;stem(n3,h3,'.k'); h4=[0,h3];h4(length(h4))=[];n4=n3; subplot(4,2,4);stem(n,x,' * k');hold on;stem(n4,h4,'.k'); h5=[0,h4];h5(length(h5))=[];n5=n4; subplot(4,2,6);stem(n,x,' * k');hold on;stem(n5,h5,'.k'); n6=-n;nmin=min(n1)-max(n6); nmax=max(n1)-min(n6); n=nmin:nmax; y=conv(x,h);%M ATLAB下卷积和计算函数 subplot(4,2,8); stem(n,y,'.k');
运算过程如图1—27所示。
图1—27 图解法实现卷积和的运算
3)列表法
列表法是计算两个有限长序列离散卷积和的方法,该方法具有简单、有效的特点。
列表法的实现:将参加卷积运算的两个序列作为表中的行和列,然后将行与列的元素对应相乘作成乘积表,最后将乘积表中对角线的元素相加得到卷积和的结果。
例1—16 将例1—14用列表法实现卷积和的计算。
解: 根据列表法规则,实现如图1—28所示。
图1—28 列表法实现序列卷积和的计算
因此,卷积和 y ( n )=[1, 3, 7, 6, 4], 0≤ n ≤4。
4)相乘对位相加法
相乘对位相加法适用于有限长序列的卷积和计算。
将参加卷积和运算两个序列的序列值按乘法运算的竖式排列,然后按照从左到右的顺序逐项将竖式相乘的乘积对位相加,所得结果就是离散卷积和。
例1—17 将例1—14用相乘对位相加法实现。
解: 根据相乘对位相加法运算规则,得如图1—29所示的乘法运算。
图1—29 相乘对位相加法实现序列卷积和的计算
因此,利用相乘对位相加法得卷积和运算结果为
y ( n ) = [1, 3, 7, 6, 4], 0≤ n ≤4
在信号的分析与处理中,经常需要研究两个信号的相似性,或者一个信号与其自身的时延信号的相似性,这就引出了相关函数。相关函数是研究随机信号的一种重要运算。
相关函数有两种形式:自相关函数(auto-correlation function, ACF)和互相关函数(cross-correlation function, CCF)。
自相关函数只涉及一个信号,提供时域中信号结构或其行为的有关信息。自相关函数在检测隐藏的周期信号时具有特殊的用途。
自相关函数的定义为
互相关函数是两个信号之间的相似性或共享性的量度。互相关可以应用于互谱分析、噪声信号的检测与恢复、模式匹配以及延迟测量等。
互相关函数的定义为
通过简单的变量代换,就可以得到互相关函数的另一种定义
上述的两种定义称为序列信号 x ( n )和 y ( n )的互相关函数,式(1.32)和式(1.33)中, r xy ( m )不能写成 r yx ( m ),因为
r yx ( m )称为信号 y ( n )和 x ( n )的互相关函数。
例1—18 计算序列 x ( n )和 y ( n )的自相关和互相关函数。
其中 x ( n )=sin(π n /8+π/3)+2·cos(π n /6), y ( n )= x ( n )+ w ( n ), w ( n )为均值为0方差为1的白噪声。
解: MATLAB实现程序如下:
n=[1:50]; x=sin(pi/8 * n+pi/3)+2 * cos(pi/6 * n); w=randn(1,length(n)); y=x+w; rxx=xcorr(x); rxy=xcorr(x,y); ryy=xcorr(y); figure(1),plot(rxx),grid; figure(2),plot(rxy),grid; figure(3),plot(ryy),grid; figure(4),plot(y);
运行结果如图1—30所示。
图1—30 相关计算
例1—19 计算数字0的语音信号的自相关函数以及其与数字2的语音信号的互相关函数。
解: MATLAB实现程序如下:
[au0,fs,bits]=wavread('au0.wav'); [au2,fs,bits]=wavread('au2.wav'); au00=xcorr(au0); au02=xcorr(au0,au2); subplot(2,2,1) plot(au0);ylabel('au0'); subplot(2,2,2) plot(au2);ylabel('au2'); subplot(2,2,3) plot(au00);ylabel('ACF 0f au0'); subplot(2,2,4) plot(au02);ylabel('CCF 0f au0 and au2');
运行结果如图1—31所示。
图1—31 数字0的语音信号的自相关函数以及其与数字2的语音信号的互相关函数