- 數字信號處理原理與實踐(修訂版)
- 劉紀紅等編著
- 3469字
- 2023-07-17 20:43:05
1.1.3 序列的運算
在離散時間信號與系統的分析與處理中,經常需要對信號進行操作,也就是對序列進行運算,常用的序列運算有序列的加與乘、移位與翻轉、尺度變換、累加、差分、卷積和以及相關等運算,所有的這些序列運算都可以利用MATLAB加以實現。
1.序列的加法與乘法
序列的加法與乘法是指兩個或兩個以上的序列對應序號的樣點值相加或相乘的運算。y1(n)=x1(n)+x2(n),y1(n)就是兩個序列x1(n)、x2(n)對應項相加后得到的和序列。y2(n)=x1(n)·x2(n),y2(n)就是兩個序列x1(n)、x2(n)對應項相乘后得到的積序列。
例1-1 利用MATLAB實現下面兩個序列的相加運算。
x1(n)=[1, 0.7, 0.4, 0.1, 0, 0.1], n=[1, 2, 3, 4, 5, 6],
x2(n)=[0.2, 0.1, 0.3, 0.5, 0.7, 0.9, 1], n=[2, 3, 4, 5, 6, 7, 8],
求x(n)=x1(n)+x2(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)
例1-2 利用MATLAB實現下面兩個正弦序列的相加運算。
x1(n)=sin(0.06πn), x2(n)=sin(0.24πn)
求y(n)=x1(n)+x2(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)=en, 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的語音信號相乘后的信號序列圖
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 原始語音信號與移位相加后的語音信號
3.序列的翻轉
序列的翻轉是以n=0的縱軸為對稱軸,將序列左右兩邊加以對調,即

序列y(n)就是x(n)的翻轉序列。
例1-7 已知函數x(n)=en-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 序列的翻轉
4.序列的尺度變換
序列的尺度變換包括幅度尺度變換和時間尺度變換。
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π50n),y(n)=sin(4π50n)。
解: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π50n),y(n)=sin(π50n)。
解: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的語音信號的幅度尺度和時間尺度的變換結果
5.序列的累加
一個序列x(n)的累加序列定義為

表示在某一n0時刻上的序列值y(n0)等于這一時刻上x(n0)的值以及n0以前時刻所有序列值的和。
例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的語音信號的累加和信號
6.差分運算
差分運算有兩種,即向前差分和向后差分。
一階向前差分:

一階向后差分:

可見,
▽x(n)=Δx(n-1)
二階向前差分:

7.序列的卷積和
在連續系統的分析與處理中,系統零狀態響應常常利用卷積積分來求解。對于離散系統來說,卷積和是求解離散時間線性移不變系統零狀態響應的有效方法。卷積和定義為

其中“?”為卷積和運算符號,通常這種離散卷積也叫做離散線性卷積。如果序列x(n)為有限長序列,序列長度為L1,序列的取值區間為[N1,N2]。序列h(n)為有限長序列,序列長度為L2,序列的取值區間為[N3,N4],即在卷積和的運算中,
x(m): N1≤m≤N2
h(n-m): N3≤n-m≤N4
則

可見兩個序列的卷積和序列y(n)也是有限長序列,其取值區間為[N1+N3,N2+N4],序列y(n)的長度為L1+L2-1。
下面具體討論幾種離散卷積的運算方法。
1)公式法
公式法就是直接利用卷積和的定義來計算離散卷積的方法。
例1-14 已知兩個序列
x(n)=2n, 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
8.序列信號的相關
在信號的分析與處理中,經常需要研究兩個信號的相似性,或者一個信號與其自身的時延信號的相似性,這就引出了相關函數。相關函數是研究隨機信號的一種重要運算。
相關函數有兩種形式:自相關函數(auto-correlation function, ACF)和互相關函數(cross-correlation function, CCF)。
自相關函數只涉及一個信號,提供時域中信號結構或其行為的有關信息。自相關函數在檢測隱藏的周期信號時具有特殊的用途。
自相關函數的定義為

互相關函數是兩個信號之間的相似性或共享性的量度。互相關可以應用于互譜分析、噪聲信號的檢測與恢復、模式匹配以及延遲測量等。
互相關函數的定義為

通過簡單的變量代換,就可以得到互相關函數的另一種定義

上述的兩種定義稱為序列信號x(n)和y(n)的互相關函數,式(1.32)和式(1.33)中,rxy(m)不能寫成ryx(m),因為

ryx(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的語音信號的互相關函數