井敏英， 龍姝明

(陜西理工大學 物理與電信工程學院， 陜西 漢中 723000)

序列卷積和計算新方法

井敏英， 龍姝明

(陜西理工大學 物理與電信工程學院， 陜西 漢中 723000)

求卷積和是在時域中計算離散線性時不變系統(tǒng)響應(yīng)的重要方法，但長序列卷積和計算的時間復雜度和空間復雜度都很高，實現(xiàn)計算時需要極大內(nèi)存的計算機系統(tǒng)，而且計算耗時很長。針對這一問題，提出計算長序列卷積和的無數(shù)據(jù)移動的矢量標積算法，并且在MATLAB環(huán)境中編程實現(xiàn)了算法。研究發(fā)現(xiàn)，無數(shù)據(jù)移動的矢量標積算法明顯地降低了長序列卷積和計算的空間復雜度和時間復雜度。用無數(shù)據(jù)移動的矢量標積算法編程，在小內(nèi)存計算機系統(tǒng)上可以快速計算長序列卷積和。

線性時不變系統(tǒng)； 序列； 卷積和； 矢量標積； MATLAB程序

線性時不變(Linear Time-Invariant，LTI)系統(tǒng)理論應(yīng)用非常廣泛[1]。離散LTI系統(tǒng)對輸入信號(序列)的響應(yīng)，在時域中計算的步驟是，先計算系統(tǒng)的單位脈沖響應(yīng)(序列)，再計算輸入信號(序列)與單位脈沖響應(yīng)(序列)的卷積和[2]。文獻[3-8]給出了卷積和計算的解釋，文獻[9-10]討論了卷積的應(yīng)用。但是這些方法計算的空間復雜度高，長序列卷積和計算在小內(nèi)存計算機系統(tǒng)上很難實現(xiàn)。事實上，長序列卷積和計算是信號處理的難題。

本文對幾種流行的序列卷積和計算方法進行了研究，分析了各自的優(yōu)缺點，通過深入研究提出了序列卷積和計算的無數(shù)據(jù)移動的矢量標積新算法。此算法計算空間復雜度極低，時間復雜度相對現(xiàn)有部分算法明顯有所降低。無數(shù)據(jù)移動的矢量標積算法使得長序列卷積和的計算在一般PC機上很容易編程實現(xiàn)。

1 三種常用的序列卷積和計算方法特點

長序列卷積和計算的解析法與圖解法，本質(zhì)上是依據(jù)序列卷積和定義進行計算。

方法一：按卷積和定義計算兩序列卷積和。

設(shè)x(k)和y(k)長度分別為n和m，k=0，1，…，n+m-2，其卷積和定義為

(1)

顯見，可以用多次乘法和加法計算兩序列的卷積和。

卷積和計算耗時近似為4n·m次乘法耗時(包含n·m次乘法、n·m次加法、n·m次數(shù)值比較和n·m次邏輯判斷所花費的時間代價)，計算占用存儲空間16(n+m)字節(jié)。

方法特點：時間復雜度極高，空間復雜度很低，編程實現(xiàn)有點難度。

方法二：將兩序列轉(zhuǎn)換為矩陣，利用矩陣乘法完成序列卷積和計算。

我們以x(k)={a，b}和y(k)={u，v，w}為例說明用矩陣乘法計算x(k)與y(k)的卷積和?？紤]到兩序列卷積和長度為2+3-1=4，首先由較長序列{u，v，w}構(gòu)造4列2行矩陣，可通過給較長序列末尾補零得到矩陣第一行，補零個數(shù)等于較短序列的長度減1，在這里補零個數(shù)為2-1=1，然后循環(huán)右移一個數(shù)據(jù)存儲位置得到第二行，依次類推來構(gòu)造較長序列對應(yīng)的矩陣，矩陣的行數(shù)與較短序列長度相同。最后將較短序列看成單行矩陣左乘較長序列對應(yīng)矩陣，即得到兩序列卷積和的結(jié)果序列。計算過程如下：

(2)

如果兩序列長度分別為n和m(設(shè)n≤m)，則采用矩陣乘法算法計算序列卷積和要完成的乘法和加法次數(shù)均為2n(n+m)次(含構(gòu)造矩陣時在內(nèi)存中多次循環(huán)移動數(shù)據(jù)消耗的時間代價)。計算需要內(nèi)存8((n+1)·m+2n)字節(jié)。

方法特點：算法的空間復雜度非常高，時間復雜度較低，編程難度很低，易理解。

方法三：利用快速傅里葉變換計算卷積和。

設(shè)要用快速傅里葉變換(FFT)計算長度為n和長度為m的兩序列x(k)和y(k)的卷積和，計算步驟是：

(1)給序列x(k)末尾補m-1個0成為序列x0，給序列y(k)末尾補n-1個0成為序列y0；

(2)分別對x0和y0取FFT得到復數(shù)頻率序列X和Y；

(3)計算頻域序列X和Y的數(shù)組乘法Z=X·Y(即對應(yīng)位置元素相乘)；

(4)計算復數(shù)序列Z的逆傅里葉變換(IFFT)，得到x(k)與y(k)的卷積和。

可用(3)式描述計算步驟

X=FFT[x0]，Y=FFT[y0]，Z=X·Y，x(k)*y(k)=IFFT[Z]，

(3)

設(shè)N=n+m-1，采用FFT算法計算兩序列卷積和，需要完成4N(1+3log2N)次實數(shù)乘法和加法，耗時代價近似為8N(1+3log2N)次實數(shù)乘法的耗時，計算需要的內(nèi)存大約為64N字節(jié)。

方法特點：時間復雜度極低，空間復雜度高于定義法又遠低于矩陣乘法，編程簡單，但初學者理解卷積和計算過程比較困難。

上述計算卷積和的3種算法，方法一適合理論教學，無工程實踐價值；方法二用于長序列卷積和計算非常困難，不能用于實踐；方法三計算長序列卷積和最快，有實踐應(yīng)用價值，但即便是短序列卷積和，手工計算也極其困難。我們期盼尋求長或短序列卷積和計算都很快、很省空間、又易于理解的卷積和計算方法。研究發(fā)現(xiàn)，無數(shù)據(jù)移動的矢量標積算法就是我們期盼的卷積和計算新方法。

2 序列卷積和計算的新方法——無數(shù)據(jù)移動的矢量標積算法

設(shè)要計算長度為n和長度為m的兩序列x(k)和y(k)的卷積和(設(shè)n≤m)，深入研究發(fā)現(xiàn)，用無數(shù)據(jù)移動的矢量標積算法，可以快速計算序列卷積和，計算步驟是：

(1)將較短序列x的數(shù)據(jù)位置反序得到xf；

(2)將較長序列y首尾都放置n-1個0數(shù)據(jù)，創(chuàng)建一個長度為m+2n-2的序列yoo；

(3)設(shè)j=1，2，…，n+m-1，從yoo的第j個位置開始順序取出n個數(shù)據(jù)記為yoo(j：j+n-1)，與xf做矢量標積給出卷積和結(jié)果的第j個數(shù)據(jù)，這一步需要重復計算n+m-1次，即

z=x*y，z(j)=xf·yoo(j：j+n-1)，j=1，2，…，n+m-1。

(4)

在MATLAB 2012a軟件環(huán)境下，具體計算步驟如下(其中含有程序語句)：

(1)計算兩序列長度并將第一序列反序再轉(zhuǎn)置為列矢量

n=length(x)； m=length(y)； L=n+m-1； xf=x(n：-1：1)′ ；

(2)給y前后各補n-1個0存于yoo中

yoo=zeros(1，L+n-1)； yoo(n：L)=y；

(3)從yoo的第j個位置開始取n個數(shù)據(jù)與xf做矢量標(或點)積運算，作為卷積計算結(jié)果的第j個數(shù)據(jù)z(j)，j=1，…，L，重復做L次標(或點)積運算完成卷積和計算，

for j=1：L； z(j)=yoo(j：j+n-1)*xf； end %此處的*代表矢量標積運算

以序列x={a，b}(n=2)與序列y={u，v，w}(m=3)為例，演示無數(shù)據(jù)移動的矢量標積算法計算卷積和的手動步驟：

(3) {0,u,v,w,0}

{b,a}z(1)={b,a}·{0,u}=au

{b,a}z(2)={b,a}·{u,v}=bu+av

{b,a}z(3)={b,a}·{v,w}=bv+aw

{b,a}z(1)={b,a}·{w,0}=bw

順序從首尾都補0的序列中取n個數(shù)據(jù)與xf做矢量標量積得到兩序列的卷積和，即

z=x*y={au，bu+av，bv+aw，bw}。

特別注意，這里討論的矢量標積算法編程實現(xiàn)時，不需要在內(nèi)存中移動數(shù)據(jù)，因為移動數(shù)據(jù)是要耗時耗能的舉動，正因為這樣稱新方法為無數(shù)據(jù)移動的矢量標積算法。

采用“無數(shù)據(jù)移動的矢量標積”算法計算兩長序列的卷積和，需要完成n(n+m-1)次實數(shù)乘法和加法，耗時代價近似為2n(n+m)次實數(shù)乘法的耗時，計算需要的內(nèi)存為16(m+2n)字節(jié)。

例估計用無數(shù)據(jù)移動的矢量標積算法對30 min單聲道錄音數(shù)據(jù)進行濾波的時空復雜度。

如果用無數(shù)據(jù)移動的矢量標積算法編程計算卷積和來完成對聲音數(shù)據(jù)進行濾波(注意n=1024，m=22 050×1800)，要完成的乘法次數(shù)為2n(n+m)= 8.13×1010次，計算過程需要內(nèi)存16(m+2n)=0.64 GB。相對于卷積和定義法算法，無數(shù)據(jù)移動的矢量標積算法的時間復雜度至少降低50%，但空間復雜度沒有變化。在一般PC機上，這一聲音濾波的無數(shù)據(jù)移動的矢量標積算法非常容易實現(xiàn)，因此該方法既可以用于教學(易于理解)又可以用于工程實踐解決信號濾波問題。

無數(shù)據(jù)移動的矢量標積算法特點：時間復雜度和空間復雜度都相對較低，易于理解易于編程實現(xiàn)，既可用于短序列卷積和的手工計算，又可以用于長序列卷積和的編程上機計算。

3 長序列卷積計算方法時間復雜度實驗比較

為了定量比較卷積和的4種算法的時間復雜度，我們在CPU為i5-4590，頻率3.3 GHz，內(nèi)存16 GB的計算機系統(tǒng)上用MATLAB編寫如下程序來體驗長序列卷積和4種計算方法耗時情況(也即時間復雜度)：

clear； clc； f=22050； t=1800；

m=f*t； n=1024； L=n+m-1；

x=rand(n，1)； y=rand(m，1)；

t0=tic； z0=conv(x，y)； dt=toc(t0)

% ========= define method ==========

t1=tic； z1=zeros(L，1)；

for k=1：L

z1(k)=0；

for p=1：n

if (k+1-m<=p)&&(p<=k)

z1(k)=z1(k)+x(p).*y(k+1-p)；

end

dt=toc(t1)； error=max(abs(z0-z1))；

disp([′time(define)=′ num2str(dt) ′s， error=′ num2str(error)])

% ========= FFT method ==============

t1=tic； x0=zeros(L，1)； y0=x0；

x0(1：n)=x； y0(1：m)=y；

X=fft(x0，L)； Y=fft(y0，L)； Z=X.*Y； z2=ifft(Z，L)；

dt=toc(t1)； error=max(abs(z2-z0))；

disp([′time(FFT)=′ num2str(dt) ′s， error=′ num2str(error)])

% ==========Vector Scale Product==========

t1=tic； xf=x(n：-1：1).′ ； yoo=zeros(L+n-1，1)；

yoo(n：L)=y； z3=zeros(L，1)；

for k=1：L； z3(k)=xf*yoo(k：k+n-1)；

end

dt=toc(t1)； error=max(abs(z3-z0))；

disp([′time(VSP)= ′ num2str(dt) ′s， error=′ num2str(error)])

其中長度n=1024的隨機實數(shù)序列模擬低通濾波器的單位脈沖響應(yīng)序列，長度m=22 050×1800=3.969×107的隨機實數(shù)序列模擬30 min單聲道語音采樣序列，卷積結(jié)果序列模擬實驗數(shù)據(jù)濾波結(jié)果。程序中用到4種卷積和計算方法，分別是：(1)MATLAB內(nèi)部卷積計算函數(shù)conv；(2)卷積定義方法；(3)FFT算法；(4)無數(shù)據(jù)移動的矢量標積算法。程序運行結(jié)果見表1。

表1 4種卷積和計算方法耗時實驗數(shù)據(jù)表

由于16 GB計算機內(nèi)存限制和大數(shù)據(jù)量限制(4千萬數(shù)據(jù)量)，矩陣乘法算法未納入實驗。實際內(nèi)存占用，是用任務(wù)管理器讀出加載MATLAB(占用內(nèi)存0.52 GB)后再運行卷積和計算程序時，MATLAB占用內(nèi)存的峰值減去0.52 GB測算出來的。

表1實驗數(shù)據(jù)與理論估算，數(shù)量級是吻合的，實驗還表明，F(xiàn)FT算法的計算耗時沒有理論預期那么短，內(nèi)存占用比理論預期還要高。

由表1實驗數(shù)據(jù)表明，我們提出的卷積計算新方法(即無數(shù)據(jù)移動的矢量標量積算法)確實綜合性能優(yōu)秀，既省空間又省時間，而且算法容易理解。這意味著，對短序列可快速手工計算、對長序列又可以快速編程上機計算，教學和工程應(yīng)用兼顧。由表1的數(shù)據(jù)轉(zhuǎn)換成文字，我們得到表2。

表2 3種卷積和計算方法特點比較表

4 結(jié) 語

表2的比較，全面地刻畫了本文提出的新卷積和計算方法(即無數(shù)據(jù)移動的矢量標量積算法)的特點。新卷積計算方法，算法本身的特點可以用“反序補零標量積”7個字概括，算法應(yīng)用及其理論特點也可用“省時省地易實現(xiàn)”7個字來概括。

當然，無數(shù)據(jù)移動的矢量標量積算法與世界級大師們凝練在MATLAB中的內(nèi)部conv、Mathematica內(nèi)部函數(shù)ListConvolve中的優(yōu)秀卷積算法還有巨大的差距。從科學研究層面看，無數(shù)據(jù)移動的矢量標量積算法僅僅是個初級的中間產(chǎn)物，還需要我們?nèi)ゲ粩嗟嘏μ剿鳌?/p>

[1] 吳大正.信號與線性系統(tǒng)分析[M].北京：高等教育出版社，2009.

[2] 高西全，丁玉美.數(shù)字信號處理[M].西安：西安電子科技大學出版社，2013.

[3] 楊永生，趙梅.從時域和頻域兩種角度探討卷積積分[J].通信技術(shù)，2010，43(11)：165-168.

[4] 黎明.探討卷積和的求解方法[J].北京工商大學學報(自然科學版)，2005，23(2)：49-51.

[5] 陳穎頻，程祝媛，王靈芝，等.基于動態(tài)坐標和圖解法的卷積積分計算[J].閩南師范大學學報(自然科學版)，2016(1)：29-38.

[6] 李春然.基于Mathematica的卷積計算[J].現(xiàn)代電子技術(shù)，2010(19)：81-86.

[7] 王益艷.離散序列線性卷積和循環(huán)卷積的關(guān)系[J].四川文理學院學報，2015，25(5)：32-35.

[8] HITZER E.General Steerable Two-sided Clifford Fourier Transform，Convolution and Mustard Convolution[J].Advances in Applied Clifford Algebras，2016，7(9)：1-20.

[9] 柳向，王杰貴，方建華.對數(shù)據(jù)后處理雷達的分段卷積轉(zhuǎn)發(fā)干擾研究[J].火力與指揮控制，2016，41(4)：15-24.

[10] 劉鳳，伍星，潘楠，等.改進時域盲解卷積算法在軸承故障診斷中的應(yīng)用[J].機械強度，2016，38(2)：207-214.

[責任編輯：謝 平]

Method of computing convolution-summation of sequence

JING Min-ying, LONG Shu-ming

(School of Physics and Telecommunication Engineering, Shaanxi University of Technology, Hanzhong 723000, China)

Computing convolution-summation is a method frequently used in discrete linear time-invariant systems response in time domain computed, but computation of convolution-summation is very complicated both in time and in space, especially in computing longer sequences. This paper proposes a new method of computing finite sequence convolution-summation and actualizes the computing method in MATLAB. The result shows that the new method has significantly reduced the computing time and space complexity and achieves greater speed in computing convolution-summation of sequence in small memory computer system.

linear time-invariant systems; series; convolution-summation; vector scalar product; MATLAB procedure

O173

A

2096-3998(2017)05-0081-05

2017-03-10

2017-07-05

井敏英(1978—)，女，陜西省富平縣人，陜西理工大學講師，碩士，主要研究方向為信號處理；龍姝明(1955—)，男，陜西省城固縣人，陜西理工大學教授，主要研究方向為量子物理、計算物理、數(shù)據(jù)處理。