井敏英， 龍姝明

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

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

井敏英， 龍姝明

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

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

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

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

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

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

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

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

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

(1)

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

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

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

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

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

(2)

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

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

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

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

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

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

(3)計(jì)算頻域序列X和Y的數(shù)組乘法Z=X·Y(即對(duì)應(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ù)乘法和加法，耗時(shí)代價(jià)近似為8N(1+3log2N)次實(shí)數(shù)乘法的耗時(shí)，計(jì)算需要的內(nèi)存大約為64N字節(jié)。

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

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

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

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

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

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

(3)設(shè)j=1，2，…，n+m-1，從yoo的第j個(gè)位置開(kāi)始順序取出n個(gè)數(shù)據(jù)記為yoo(j：j+n-1)，與xf做矢量標(biāo)積給出卷積和結(jié)果的第j個(gè)數(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ì)算步驟如下(其中含有程序語(yǔ)句)：

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

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

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

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

(3)從yoo的第j個(gè)位置開(kāi)始取n個(gè)數(shù)據(jù)與xf做矢量標(biāo)(或點(diǎn))積運(yùn)算，作為卷積計(jì)算結(jié)果的第j個(gè)數(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)為例，演示無(wú)數(shù)據(jù)移動(dòng)的矢量標(biāo)積算法計(jì)算卷積和的手動(dòng)步驟：

(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個(gè)數(shù)據(jù)與xf做矢量標(biāo)量積得到兩序列的卷積和，即

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

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

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

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

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

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

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

為了定量比較卷積和的4種算法的時(shí)間復(fù)雜度，我們?cè)贑PU為i5-4590，頻率3.3 GHz，內(nèi)存16 GB的計(jì)算機(jī)系統(tǒng)上用MATLAB編寫(xiě)如下程序來(lái)體驗(yàn)長(zhǎng)序列卷積和4種計(jì)算方法耗時(shí)情況(也即時(shí)間復(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)])

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

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

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

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

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

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

4 結(jié) 語(yǔ)

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

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

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

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

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

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

[5] 陳穎頻，程祝媛，王靈芝，等.基于動(dòng)態(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] 柳向，王杰貴，方建華.對(duì)數(shù)據(jù)后處理雷達(dá)的分段卷積轉(zhuǎn)發(fā)干擾研究[J].火力與指揮控制，2016，41(4)：15-24.

[10] 劉鳳，伍星，潘楠，等.改進(jìn)時(shí)域盲解卷積算法在軸承故障診斷中的應(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)樾盘?hào)處理；龍姝明(1955—)，男，陜西省城固縣人，陜西理工大學(xué)教授，主要研究方向?yàn)榱孔游锢?、?jì)算物理、數(shù)據(jù)處理。