درک ریاضی کانولوشن و روشهای محاسبه کانولوشن خطی

عمل کانولوشن (convolution)، در بسیاری از کاربردهای پردازش سیگنال دیده میشود. ریاضیات کانولوشن ریشه در عملیات روی چندجمله‌ای‌ها دارد. هدف از این مقاله، ارتقا سطح دانش و درک جزییات ریاضی کانولوشن و همچنین آشنایی با روشهای متنوع محاسبه کانولوشن خطی است.

توابع چند‌جمله‌ای

توابع چند‌جمله‌ای عبارتهایی شامل جمع بخشهای مختلف است که هر بخش شامل یک یا چند متغیر به توان یک مقدار نامنفی به همراه یک ضریب تغییر مقیاس است. جمع، تفریق و ضرب چندجمله‌ای‌ها میسر و ممکن است.

توابع چندجمله‌ای میتوانند دارای یک یا چند متغیر باشند. برای مثال، عبارت چندجمله‌ای که در ادامه می‌آید، تابعی از متغیر x است. این عبارت شامل 3 بخش است که هر بخش با یک ضریب، تغییر مقیاس داده شده است.

$f(x)=x^2+2x+1$

عبارت چندجمله‌ای بعدی دارای دو متغیر x و y است.

$f(x,y)=2x^4-4x^2y+3xy^2+8y^2+5$

نمایش توابع چندجمله‌ای تک-متغیره

توابع چندجمله‌ای که یک متغیر دارند در اینجا مدنظر است. به طور کلی چندجمله‌ای تک-متغیره (به نام x) به صورت جمع بخشهای مختلف که شامل ضرایب $a_0,a_1,a_2,...,a_{n-1}$ هستند، مشخص میشود:

$f(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_{n-1}x^{n-1}$

درجه یا مرتبه تابع چندجمله‌ای، بالاترین توان متغیر x که ضریب غیرصفر دارد، است. معادله بالا میتواند به شکل زیر نیز نوشته شود:

$f(x)=\sum_{i=0}^{n-1}a_ix^i$

این نمایش میتواند حتی به فرم بردار ضرایب $a=[a_0, a_1,a_2,...,a_{n-1}]$ نیز باشد.

چندجمله‌ای‌ها میتوانند با استفاده از ریشه‌هایشان که به فرم ضرب بخشهای خطی هستند، نیز نمایش داده شوند که این موضوع در ادامه بررسی میشود.

ضرب چندجمله‌ای‌ها و کانولوشن خطی

همانطور که قبلا نشان داده شد، عملیات ریاضی همچون جمع، تفریق و ضرب بر روی توابع چندجمله‌ای قابل انجام است. جمع یا تفریق چندجمله‌ای‌ها آسان و سرراست است. اما ضرب آنها با توجه به موضوع مورد بحث در اینجا، مورد توجه است.

محاسبه ضرب دو چندجمله‌ای با بردار ضرایب $a=[a_0,a_1,a_2,...,a_n]$ و $b=[b_0,b_1,b_2,...,b_n]$ مد نظر است. نمایش رایج این چندجمله‌ایها به شکل زیر است:

$p(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_nx^n\\ \\ {\ \ \ \ }q(x)=b_0+b_1x+b_2x^2+b_3x^3+...+b_nx^n$

ضرب آنها عبارت است از:

$p(x)q(x)=a_0b_0+(a_1b_0+b_1a_0)x+...+a_nb_mx^{n+m}$

یا به طور معادل:

$p(x)q(x)=(p*q)(x)=a.b=[x_0,c_1,c_2,...,c_{n+m}] \\ \\ {\ \ \ \ } c_k=\sum_{i,j:i+j=k}a_ib_j \ \ \ \ \ \ k=0,1,...,n+m$

از آنجاییکه اندیسها از رابطه $i+j=k$ تبعیت میکنند، تغییر اندیس $j$ به $k-i$ منجر به رابطه زیر میشود:

$c_k=\sum_{i=-\infty}^{\infty}a_ib_{k-i} \ \ \ \ \ k=0,1,...,n+m$

که وقتی بر حسب اندیسها نوشته میشود، رایج‌ترین فرم ظاهری دیده شده در متون پردازش سیگنال بدست می‌آید:

$c[k]=\sum_{i=-\infty}^{\infty}a[i]b[k-i] \ \ \ \ \ k=0,1,...,n+m$

این عمل تحت عنوان کانولوشن (کانولوشن خطی به طور دقیق) معرفی میشود که علامت آن * است. این عمل، خیلی به سایر عملیات روی بردارها نزدیک است – همچون همبستگی متقابل (cross-correlation)، خودهمبستگی (auto-correlation) و محاسبات میانگین متحرک (moving average).

بنابراین زمانیکه کانولوشن را محاسبه میکنیم به نوعی در حال ضرب دو چندجمله‌ای هستیم. دقت کنید که اگر چندجمله‌ایها شامل N و M بخش باشند، ضرب آنها M+N-1 بخش تولید میکند.

الگوریتم سرراست (در پایین) محاسبه عبارت ضرب به لحاظ پیچیدگی محاسباتی و زمانی از مرتبه $O(n^2)$ است. یک الگوریتم بهتر مورد نیاز است. به نظر میرسد که اگر به چندجمله‌ایها به فرم ضرب فاکتورهای خطی با ریشه‌های مختلط نگاه کنیم، بار محاسباتی کاهش خواهد یافت. این فرم، پایه روش مبتنی بر تبدیل فوریه سریع یا FFT است که ریشه nام مقدار واحد به طور بازگشتی برای محاسبه در حوزه فرکانس مورد استفاده قرار میگیرد.

for i = 0:n+m,

ci = 0;

for i = 0 :n,

for k = 0 to m,

c[i+k] = c[i+k] + a[i] · b[k];

ماتریس توپلیتز و کانولوشن

عمل کانولوشن دو دنباله میتواند به صورت ضرب دو ماتریس به شکل زیر در نظر گرفته میشود. با فرض داشن یک سیستم خطی تغییرناپذیر با زمان یا LTI با پاسخ ضربه h[n] و دنباله ورودی x[n]، خروجی سیستم y[n] از طریق کانولوشن دنباله ورودی و پاسخ ضربه حاصل میشود:

$y[k]=h[n]*x[n]=\sum_{i=-\infty}^{\infty}x[i]h[k-i]$

که دنباله x[n] طول N و h[n] طول M دارد.

فرض کنید که دنباله h[n] به طول 4 به شکل $h[n]=[h_0,h_1,h_2,h_3]$ است و دنباله x[n] با طول 3 به شکل $x[n]=[x_0,x_1,x_2]$ کانولوشن h[n]*x[n] به شکل زیر است:

$y[k]=h[n]*x[n] = \sum_{i=-\infty}^{\infty}x[i]h[k-i] \ \ \ \ k=0,1,2,3,4,5$

محاسبه هر نمونه از کانولوشن به شرح زیر است:

$\\ y[0]=\sum_{i=-\infty}^{\infty}x[i]h[-i]=x[0]h[0]+0+0 \\ \\ y[1]=\sum_{i=-\infty}^{\infty}x[i]h[1-i]=x[0]h[1]+x[1]h[0]+0 \\ \\ y[2]=\sum_{i=-\infty}^{\infty}x[i]h[2-i]=x[0]h[2]+x[1]h[1]+x[2]h[0] \\ \\ y[3]=\sum_{i=-\infty}^{\infty}x[i]h[3-i]=x[0]h[3]+x[1]h[2]+x[2]h[1] \\ \\ y[4]=\sum_{i=-\infty}^{\infty}x[i]h[4-i]=x[1]h[3]+x[2]h[1]+0 \\ \\ y[5]=\sum_{i=-\infty}^{\infty}x[i]h[5-i]=x[2]h[3]+0+0$

به نتیجه بالا دقت کنید. میبینید که هر سری x[n] و h[n] چگونه از جهت مخالف با یکدیگر ضرب میشوند. این مساله، علت معکوس کردن یک دنباله و جابه‌جایی زمانی آن (که معمولا به صورت گرافیکی در متون پردازش سیگنال نمایش داده میشود) را برای محاسبه کانولوشن، نشان میدهد.

بنابراین، به صورت گرافیکی، در کانولوشن، یکی از دنباله‌ها معکوس میشود. به چپ‌ترین مکان ممکن جابه‌جا میشود به نحوی که همپوشانی بین دو دنباله وجود نداشته باشد. حال در هر گام یک نمونه به راست حرکت میکنیم. در هر گام، بخشهای همپوشان h و x در هم ضرب شده و نمونه‌های بدست آمده با هم جمع میشوند. این فرآیند تا زمانیکه دنباله h به راست‌ترین مکان ممکن برسد و دیگر همپوشانی با x نداشته باشد، ادامه پیدا میکند.

نوشتن نتیجه نهایی معادله بالا به فرم ماتریسی به شکل زیر است:

$\begin{bmatrix} y[0]\\ y[1]\\ y[2]\\ y[3]\\ y[4]\\ y[5] \end{bmatrix} = \begin{bmatrix} h[0]&0 &0 \\ h[1]&h[0] &0 \\ h[2] &h[1] &h[0] \\ h[3]&h[2] &h[1] \\ 0&h[3] &h[2] \\ 0 & 0 & h[3] \end{bmatrix}\begin{bmatrix} x[0]\\ x[1]\\ x[2] \end{bmatrix}$

زمانیکه دنباله‌های h,x به شکل ماتریس نمایش داده شوند، عمل کانولوشن به طور معادل به شکل زیر قابل نمایش است:

$y=h*x=x*h = \begin{bmatrix} h_0&0 &0 \\ h_1&h_0 &0 \\ h_2 &h_1 &h_0 \\ h_3&h_2 &h_1 \\ 0&h_3 &h_2 \\ 0 & 0 & h_3 \end{bmatrix}\begin{bmatrix} x_0\\ x_1\\ x_2 \end{bmatrix}$

ماتریسی که تاخیرهای افزایشی h[n] (مورد استفاده در محاسبه کانولوشن) را نشان میدهد، فرم ماتریس توپلیتز نام دارد. ماتریس توپلیتز حاوی مقادیر ثابت در راستای قطر اصلی خود است. ماتریسهای توپلیتز در مدل کردن سیستمهایی که خواص تغییرناپذیری با زمان دارند، استفاده میشوند. خاصیت تغییرناپذیری با زمان از خود ساختار ماتریس مشهود است. از آنجاییکه یک سیستم خطی تغییرناپذیر با زمان در استفاده از کانولوشن فرض میشود، ماتریسهای توپلیتز انتخاب طبیعی ما خواهد بود. به عنوان یک نکته فرعی، توجه داشته باشید که یک فرم خاص ماتریس توپلیتز تحت عنوان ماتریس دایروی (circular matrix) وجود دارد که در کاربردهای مربوط به کانولوشن دایروی (circular) و تبدیل فوریه گسسته (DFT) استفاده میشود.

نمایش نحوه ساختن ماتریس توپلیتز T از دنباله h به شکل زیر است:

$T(h) = \begin{bmatrix} h0&0 &0 \\ h_1&h_0 &0 \\ h_2 &h_1 &h_0 \\ h_3&h_2 &h_1 \\ 0&h_3 &h_2 \\ 0 & 0 & h_3 \end{bmatrix}$

حال، کانولوشن x و h به سادگی برابر است با ضرب ماتریس توپلیتز T(h) در ماتریس (بردار) نمایش دهنده X:

$y=h*x=x*h=\begin{bmatrix} h0&0 &0 \\ h_1&h_0 &0 \\ h_2 &h_1 &h_0 \\ h_3&h_2 &h_1 \\ 0&h_3 &h_2 \\ 0 & 0 & h_3 \end{bmatrix}\begin{bmatrix} x_0\\ x_1\\ x_2 \end{bmatrix}=T(h).X$

عمل کانولوشن از این طریق را به سرعت میتوان در متلب به شکل برداری انجام داد. برای این کار از دستور toeplitz استفاده میکنیم:

y=toeplitz([h0 h1 h2 h3 0 0],[h0 0 0])*x.';

تعریف

با فرض سیستم LTI با پاسخ ضربه h[n] و دنباله ورودی x[n]، خروجی سیستم برابر با y[n] است که از طریق کانولوشن دنباله ورودی و پاسخ ضربه حاصل میشود:

$y[k]=h[n]*x[n]=\sum_{i=-\infty}^{\infty}x[i]h[k-i]$

که دنباله x[n] طول N و دنباله h[n] طول M دارد.

روشهای مختلف محاسبه کانولوشن به شرح زیر است:

روش 1- نیروی سهمگین (brute-force)

معادله بالا به شکل اصلی خود قابل پیاده‌سازی است که این کار از طریق حلقه‌های for تودرتو قابل انجام است. این روش بیشترین زمان محاسباتی را در بین تمامی روشها مصرف میکند. معمولا، پیچیدگی محاسباتی این روش $O(n^2)$ است.

کد پایتون:

def conv_brute_force(x,h):

“””

Brute force method to compute convolution

Parameters:

x, h : numpy vectors

Returns:

y : convolution of x and h

“””

N=len(x)

M=len(h)

y = np.zeros(N+M-1) #array filled with zeros

for i in np.arange(0,N):

for j in np.arange(0,M):

y[i+j] = y[i+j] + x[i] * h[j]

return y

کد متلب:

for i = 0:n+m,

yi = 0;

for i = 0 :n,

for k = 0 to m,

y[i+k] = y[i+k] + x[i] · h[k];

end

روش 2- استفاده از ماتریس توپلیتز

زمانیکه دنباله h[n] و x[n] به شکل ماتریس نمایش داده شوند، عمل کانولوشن خطی به طور معادل به شکل زیر بازنویسی میشود:

$y=h*x=x*h=Toeplitz(h)X$

با فرض طول 4 برای h و 3 برای x، داریم:

$y[k]=h[n]*x[n] = \sum_{i=-\infty}^{\infty}x[i]h[k-i] \ \ \ \ k=0,1,2,3,4,5$

نمایش معادل کانولوشن بالا به فرم ماتریسی به شکل زیر است:

ماتریس نشان‌دهنده تاخیرهای افزایشی h[n] که در معادله بالا استفاده شده است، همان فرم خاصی ماتریس توپلیتز است.

همانطور که قبلا اشاره شد، در متلب تابع درونی برای محاسبه فرم توپلیتز از یک بردار مشخص وجود دارد. ماتریس توپلیتز دنباله h به شکل زیر محاسبه میشود:

$Toeplitz(h)=\begin{bmatrix} h[0] &0 &0 \\ h[1]&h[0] &0 \\ h[2]&h[1] &h[0] \\ h[3]&h[2] &h[1] \\ 0& h[3] &h[2] \\ 0& 0 & h[3] \end{bmatrix}$

در تابع toeplitz، ورودی اول ستون اول ماتریس بالا است و ورودی دوم سطر اول ماتریس بالا.

y=toeplitz([h0 h1 h2 h3 0 0],[h0 0 0])*x.’;

در حالت کلی، که x و h دارای طولهای دلخواه M و N هستند، کد بالا به شکل زیر قابل اصلاح است (با فرض length(x)>=length(h)):

k=toeplitz([h zeros(1,length(x)-1)],[h(1) zeros(1,length(x)-1])*x.’

روش 3- استفاده از تبدیل فوریه سریع یا همان FFT

محاسبه کانولوشن با استفاده از FFT، مزیت کاهش بار محاسباتی برای دنباله‌های طولانی را دارد. برای محاسبه کانولوشن، ابتدا تبدیل FFT دنباله‌های x و h برای طول FFT برابر با length(x)+length(h)-1، محاسبه میشود. سپس این دو تبدیل در هم ضرب شده و تبدیل معکوس IFFT برای بازگشت به حوزه زمان بر روی دنباله حاصلضرب اعمال میشود. دقت کنید که FFT، پیاده‌سازی مستقیم کانولوشن دایروی در حوزه زمان است. در اینجا به دنبال محاسبه کانولوشن خطی با استفاده از کانولوشن دایروی (یا FFT) به همراه اضافه کردن صفر (zero-padding) به دنباله کوتاهتر هستیم. البته این روش در مقایسه با کانولوشن دایروی مقداری کاهش راندمان دارد ولی به هر حال دارای پیچیدگی محاسباتی $O(\frac{N}{log_2N})$ است که نسبت به روش brute-force راندمان بهتری دارد:

$y(n)=IFFT[FFT_L(X)*FFT_L(H)],\ \ \ \ \ L=2^{nextpower2(N+M-1)}$

معمولا، الگوریتم زیر که صفرهای اضافی در خروجی را حذف میکند، برای محاسبه کانولوشن کفایت میکند:

$y(n)=IFFT[FFT_L(X)*FFT_L(H)], \ \ \ \ L=N+M-1$

کد پایتون:

from scipy.fftpack import fft,ifft

y=ifft(fft(x,L)*(fft(h,L))) #Convolution using FFT and IFFT

کد متلب:

y=ifft(fft(x,L).*(fft(h,L)))

روشهای دیگر

اگر طول دنباله ورودی بینهایت یا خیلی طولانی باشد، همچنان که در بسیاری کاربردهای واقعی چنین است، روشهای پردازش بلوکی (block processing) همچون Overlap-Add و Overlap-Save برای محاسبه کانولوشن استفاده میشود که سریعتر و کارآمدتر است.

الگوریتمهای استاندارد برای کانولوشن سریع

الگوریتمهای استاندارد برای محاسبه کانولوشن که پیچیدگی محاسباتی را به شدت کاهش میدهند عبارت‌اند از:

– Cook-Toom Algorithm

– Modified Cook-Toom Algorithm

– Winograd Algorithm

– Modified Winograd Algorithm

– Iterated Convolution

کد متلب مقایسه روشها

روش 2 (استفاده از ماتریس توپلیتز) و روش 3 (استفاده از FFT) با روش استاندارد متلب از طریق تابع conv، در کد زیر مقایسه شده‌اند:

%Create random vectors for test

x=rand(1,7)+1i*rand(1,7);

h=rand(1,3)+1i*rand(1,3);

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

%Convolution Using Toeplitz matrix

y1=toeplitz([h zeros(1,length(x)-1)],[h(1) zeros(1,length(x)-1)])*x.’

%Convolution FFT

y2=ifft(fft(x,L).*(fft(h,L))).’

%Matlab’s standard function

y3=conv(h,x).’

نتیجه شبیه‌سازی

نتایج سه روش در زیر آمده است که همگی یکسان هستند:

منبع: www.gaussianwaves.com

درک ریاضی کانولوشن و روشهای محاسبه کانولوشن خطی

دیدگاه ها (0)

لغو پاسخ