井敏英， 龍姝明

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

序列卷積和計(jì)算新方法

井敏英， 龍姝明

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

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

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

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

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

1 三種常用的序列卷積和計(jì)算方法特點(diǎn)

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

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

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

(1)

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

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

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

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

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

(2)

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

方法特點(diǎn)：算法的空間復(fù)雜度非常高，時間復(fù)雜度較低，編程難度很低，易理解。

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

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

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

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

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

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

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

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

(3)

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

方法特點(diǎn)：時間復(fù)雜度極低，空間復(fù)雜度高于定義法又遠(yuǎn)低于矩陣乘法，編程簡單，但初學(xué)者理解卷積和計(jì)算過程比較困難。

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

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

設(shè)要計(jì)算長度為n和長度為m的兩序列x(k)和y(k)的卷積和(設(shè)n≤m)，深入研究發(fā)現(xiàn)，用無數(shù)據(jù)移動的矢量標(biāo)積算法，可以快速計(jì)算序列卷積和，計(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做矢量標(biāo)積給出卷積和結(jié)果的第j個數(shù)據(jù)，這一步需要重復(fù)計(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)境下，具體計(jì)算步驟如下(其中含有程序語句)：

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

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

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

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

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

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

以序列x={a，b}(n=2)與序列y={u，v，w}(m=3)為例，演示無數(shù)據(jù)移動的矢量標(biāo)積算法計(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

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

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

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

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

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

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

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

3 長序列卷積計(jì)算方法時間復(fù)雜度實(shí)驗(yàn)比較

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

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

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

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

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

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

表2 3種卷積和計(jì)算方法特點(diǎn)比較表

4 結(jié) 語

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

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

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

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

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

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

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

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

[7] 王益艷.離散序列線性卷積和循環(huán)卷積的關(guān)系[J].四川文理學(xué)院學(xué)報(bào)，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ù)后處理雷達(dá)的分段卷積轉(zhuǎn)發(fā)干擾研究[J].火力與指揮控制，2016，41(4)：15-24.

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

[責(zé)任編輯：謝 平]

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—)，女，陜西省富平縣人，陜西理工大學(xué)講師，碩士，主要研究方向?yàn)樾盘柼幚恚积堟?1955—)，男，陜西省城固縣人，陜西理工大學(xué)教授，主要研究方向?yàn)榱孔游锢怼⒂?jì)算物理、數(shù)據(jù)處理。