基于Ｍａｔｌａｂ的離散卷積

摘 要：卷積運(yùn)算廣泛用于通訊、電子、自動(dòng)化等領(lǐng)域的線性系統(tǒng)的仿真、分析及數(shù)字信號(hào)處理等方面。在Matlab中可以使用線性卷積、圓周卷積和快速傅里葉運(yùn)算實(shí)現(xiàn)離散卷積。線性卷積是工程應(yīng)用的基礎(chǔ)，但圓周卷積和快速傅里葉運(yùn)算實(shí)現(xiàn)線性離散卷積具有速度快等優(yōu)勢(shì)，圓周卷積采用循環(huán)移位，在Matlab中沒有專用函數(shù)，需要根據(jù)圓周卷積的運(yùn)算過程編制程序代碼；快速傅里葉運(yùn)算(FFT)是DSP的核心算法，在序列比較長(zhǎng)時(shí)FFT是一種最合適的方法，運(yùn)算速度快、程序簡(jiǎn)單，序列越長(zhǎng)其優(yōu)勢(shì)越明顯。以同一個(gè)例子介紹了進(jìn)行離散卷積仿真運(yùn)算的兩種方法與特點(diǎn)。

關(guān)鍵詞：Matlab；數(shù)字信號(hào)處理；離散卷積運(yùn)算；快速傅里葉運(yùn)算

中圖分類號(hào)：TP311文獻(xiàn)標(biāo)識(shí)碼：B

文章編號(hào)：1004-373X(2009)05-125-02

Discreted Convolution Based on Matlab

LIU Guoliang

(Luoyang Institute of Science and Technology,Luoyang,471000,China)

Abstract：The convolution operation is widely applied in the field of linear system simulations,analysis and DSP design of communication,electronics,automatic system and so on,linear convolution,circular convolution and FFT for discreted convolution emulation operation can be adopted in Matlab.Linear convolution is the foundation of project application,but circular convolution and FFT in linear convolution have advantages of fast speed.Circular convolution adopts circulating shift,there is not privately function in Matlab,it needs to draw up program code according to the operational course of circular convolution.FFT is the major operation of DSP,when sequence grows fairly,FFT is a kind of most suitable method,it is more rapidly in operational speed and program simple,the longer of the sequence the more of its advantage.With a same example,to carry out characteristic and two kinds of method of emulation operation for discreted convolution are introduced.

Keywords：Matlab;DSP;discreted convolution operation;FFT

0 引 言

在數(shù)字信號(hào)處理(DSP)系統(tǒng)中，無(wú)論在時(shí)域或頻域都離不開卷積運(yùn)算和快速傅里葉(FFT)運(yùn)算。Matlab具有強(qiáng)大的矩陣運(yùn)算能力，方便實(shí)用的繪圖功能和語(yǔ)言的高度集成性，在DSP開發(fā)中，使用Matlab可以快速對(duì)系統(tǒng)進(jìn)行仿真運(yùn)算。

設(shè)輸入信號(hào)為x(t)，其頻譜函數(shù)為X(jΩ)，該信號(hào)通過濾波器h(t)后，其輸出信號(hào)y(t)的頻譜函數(shù)Y(jΩ)是頻譜函數(shù)X(jΩ)與濾波器的頻譜函數(shù)H(jΩ)的乘積，即:

Y(jΩ)=X(jΩ)H(jΩ)(1)

該信號(hào)通過低通濾波器的示意圖如圖1所示。

而在時(shí)域，輸出信號(hào)y(t)實(shí)際是輸入信號(hào)x(t)與濾波器h(t)的卷積。就是說頻譜函數(shù)的乘積相當(dāng)于時(shí)間函數(shù)的卷積，反之亦然。即：

y(t)=x(t)*h(t)=∫∞-∞x(τ)h(t-τ)dτ(2)

圖1 信號(hào)通過低通濾波器的示意圖

由于DSP主要依靠計(jì)算機(jī)完成，而計(jì)算機(jī)無(wú)論在時(shí)域或頻域只能處理有限長(zhǎng)的離散信號(hào)，因此上述關(guān)系可表示為：

Y(ω)=X(ω)H(ω)(3)

y(n)=∑∞m=-∞x(m)h(n-m)=x(n)*h(n)(4)

式中：n表示離散時(shí)間；數(shù)字頻率ω=ΩT璼=2πf/f璼；f為輸入信號(hào)頻率；f璼,T璼為抽樣頻率和周期。

1 實(shí)現(xiàn)離散卷積的方法

可以使用線性卷積、圓周卷積實(shí)現(xiàn)離散卷積，也可以根據(jù)時(shí)域信號(hào)函數(shù)與頻譜函數(shù)的關(guān)系，使用快速傅里葉變換(FFT，IFFT)實(shí)現(xiàn)離散卷積。

1.1 圓周卷積

線性卷積是求離散系統(tǒng)響應(yīng)的主要方法之一，許多重要應(yīng)用都建立在這一理論基礎(chǔ)上,如卷積濾波等，專用函數(shù)conv(x,h)可完成線性卷積過程。

圓周卷積的計(jì)算速度要遠(yuǎn)遠(yuǎn)快于線性卷積，如果選擇圓周卷積的長(zhǎng)度L≥N1+N2-1,則可以用圓周卷積取代線性卷積。方法如下:

(1) 定義圓周卷積的長(zhǎng)度：選擇L=N1+N2-1。

(2) 將兩個(gè)序列的長(zhǎng)度都補(bǔ)足為L(zhǎng):

將長(zhǎng)為N1的序列x(n)延長(zhǎng)到L,補(bǔ)L-N1個(gè)零，

將長(zhǎng)為N2的序列h(n)延長(zhǎng)到L,補(bǔ)L-N2個(gè)零。

(3) 反轉(zhuǎn)x(n)，周期延拓為序列 (n)，取主周期。

(4) 循環(huán)移位：與線性卷積不同，圓周卷積運(yùn)算中采用的是循環(huán)移位，有限長(zhǎng)序列x(n)的循環(huán)移位定義為：

f(n)=x((n+m))璑R璑(n)

其含義如下：

① x((n+m))璑 表示 x(n) 的周期延拓序列 (n)的移位：

x((n+m))璑=(n+m)

② x((n+m))璑R璑(n)表示對(duì)移位的周期序列x((n+m))璑取主值序列，所以f(n)仍然是一個(gè)長(zhǎng)度為N的有限長(zhǎng)序列。

為此定義一個(gè)位移函數(shù)cirshiftd.m如下：

%保存函數(shù)為cirshiftd.m

function fm= cirshiftd(x,m,N);

n=［0:1:N-1］;

fm=x(mod(n-m,N)+1);

這里的x即表示對(duì)向右移位N位的周期序列x((n+m))璑取主值序列的運(yùn)算。

(5) 將序列x與h進(jìn)行循環(huán)移位生成的矩陣H的轉(zhuǎn)置矩陣相乘，即完成將x(m)和h(n-m)所有對(duì)應(yīng)項(xiàng)相乘之后相加，完成圓周卷積。

例：已知x=［1 2 3 2 ］和h=［1 1 1］，求兩序列的卷積。

用圓周卷積的計(jì)算程序如下：

%′圓周卷積′主程序

x=［1 2 3 2 ］;

h=［1 1 1］;

N=length(x)+length(h)-1;

x=［x,zeros(1,N- length(x))］;

h=［h,zeros(1,N- length(h))］;

n=［0:1:N-1］;

h=h(mod(-n,N)+1);

H=zeros(N,N);

for n=1:1:N

H(n,:)=cirshiftd(h,n-1,N);%調(diào)用位移函數(shù)

end

y=x*H′;

n0=［0:1:N-1］;

stem(n0,y);

title(′圓周卷積′);

xlabel(′( n )′);ylabel(′ y(n)′);

運(yùn)行結(jié)果：y=［1 3 6 7 5 2 ］

1.2 線性卷積的FFT算法

FFT運(yùn)算具有速度快、計(jì)算量小等優(yōu)越性，是DSP的核心算法。根據(jù)信號(hào)在時(shí)域與頻域的關(guān)系，可以用FFT計(jì)算線性卷積。方法如下:

(1) 定義長(zhǎng)度：選擇N=N1+N2-1。

(2) 求兩個(gè)要卷積序列的N點(diǎn)頻譜函數(shù)：

求X(k)=FFT［x(n)］

求H(k)=FFT［h(n)］

(3) 將兩個(gè)頻譜函數(shù)相乘：

求Y(k)=H(k)X(k)

(4) 對(duì)乘積進(jìn)行N點(diǎn)快速傅里葉反變換(IFFT)：

求y(n)=IFFT［Y(k)］

可見，只要進(jìn)行二次FFT,一次IFFT就可完成線性卷積計(jì)算。在x(n)，h(n) 兩序列長(zhǎng)度比較接近或相等的情況，當(dāng)運(yùn)算的序列長(zhǎng)度 N＞32時(shí),上述計(jì)算線性卷積的方法比直接計(jì)算線性卷積有明顯的優(yōu)越性，長(zhǎng)度值越大優(yōu)勢(shì)越明顯。

上例用FFT運(yùn)算的程序如下：

x=［1 2 3 2 ］;

h=［1 1 1］;

N= length(x) + length(h)-1;

X=fft(x,N);

H=fft(h,N);

Y=X.*H;

y=ifft(Y,N);

n=0:1:N-1;

stem(n,y);

title(′FFT 線性卷積′);

xlabel(′( n )′);ylabel(′ y(n)′);

2 結(jié) 語(yǔ)

用圓周卷積和FFT運(yùn)算離散卷積具有很大優(yōu)勢(shì)，在Matlab中圓周卷積沒有專用函數(shù)，需要自己編制，在序列比較長(zhǎng)時(shí)FFT是一種最合適的方法，可以快速得出仿真運(yùn)算的結(jié)果。選擇圓周卷積或FFT運(yùn)算線性離散卷積可提高DSP設(shè)計(jì)的效率。

參考文獻(xiàn)

［1］徐科軍.信號(hào)分析與處理［M］.北京:清華大學(xué)出版社,2006.

［2］劉順蘭,吳杰.數(shù)字信號(hào)處理［M］.西安:西安電子科技大學(xué)出版社,2003.

［3］朱仁峰.精通Matlab 7［M］.北京:清華大學(xué)出版社,2006.

［4］燕慶明.信號(hào)與系統(tǒng)［M］.北京:高等教育出版社,2004.

［5］俞卞章.數(shù)字信號(hào)處理［M］.2版.西安:西北工業(yè)大學(xué)出版社,2006.

［6］徐莉,羅新民,徐燕紅.卷積碼的Matlab仿真及其性能研究.現(xiàn)代電子技術(shù),2006,29(11):64-66.

［7］黃瓊玲,劉振興,尉宇.一種分?jǐn)?shù)階傅里葉變換快速算法的研究.現(xiàn)代電子技術(shù)，2008，31(9)：156-157.

作者簡(jiǎn)介

劉國(guó)良 男，1958年出生，計(jì)算機(jī)高級(jí)工程師。研究方向?yàn)橛?jì)算機(jī)CAD設(shè)計(jì)、電子技術(shù)應(yīng)用和無(wú)線電等。