محاسبه خودهمبستگی و ماتریس آن در متلب و پایتون

 

 

 

خودهمبستگی (Auto-correlation) که همبستگی سری هم نامیده میشود، همبستگی یک دنباله با خودش به عنوان تابعی از تاخیر زمانی است. همبستگی متقابل (Cross-correlation) عبارت جامع‌تری است که همبستگی بین دو دنباله مختلف را به عنوان تابعی از تاخیر زمانی میدهد.

با داشتن دو دنباله x و y، همبستگی متقابل در زمانهای جدا از هم به اندازه تاخیر i، توسط معادله زیر بدست می‌آید:

 

    \[ R_{xy_i} = Corr(x,y)_i = \sum_{j=-\infty}^{\infty}x_j y^{*}_{j+i} \]

 

خودهمبستگی، یک حالت خاص از همبستگی متقابل است که در آن x=y است. میتوان از روش نیروی سهمگین (Brute force) (که در آن از حلقه‌های for برای پیاده‌سازی معادله بالا استفاده میشود) برای محاسبه دنباله خودهمبستگی استفاده کرد. هر چند که، روشهای جایگزین دیگری نیز در دسترس هست.

 

 

روش اول: خودهمبستگی با استفاده از توابع آماده در متلب و پایتون

 

 

متلب (MATLAB)

 

برای یک بردار N-بعدی داده شده x، دستور xcorr(x,x) در متلب یا به شکل ساده‌تر، xcorr(x)، دنباله خودهمبستگی را ایجاد میکند. مثال زیر یک نمونه از این محاسبه است:

x = [1 2 3 4]

acf = xcorr(x)

acf = 4 11 20 30 20 11 4

پایتون (Python)

 

در پایتون، خودهمبستگی دنباله یک-بعدی را میتوان بوسیله تابع numpy.correlate محاسبه کرد. در این تابع، پارامتر mode باید در حالت ‘full’ قرار بگیرد که در شرایطی که میخواهیم خودهمبستگی را به صورت تابعی از تاخیر محاسبه کنیم، لازم است. یک نمونه از این محاسبه در پایتون در ادامه آمده است:

 

import numpy as np

x = np.asarray([1,2,3,4])

np.correlate(x, x, mode=’full’)

# output: array([4, 11, 20, 30, 20, 11, 4])

 

 

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

 

دنباله خودهمبستگی را میتوان به صورت کانولوشن (Convolution) بین دو دنباله اصلی و نسخه معکوس (flipped) شده مزدوج همان دنباله، محاسبه کرد. عمل مزدوج کردن در صورتیکه دنباله اصلی، حقیقی باشد، نیازی نیست.

 

 

متلب

 

x = [1 2 3 4]

acf = conv(x,fliplr(conj(x)))

acf = 4 11 20 30 20 11 4

 

 

پایتون

 

import numpy as np

x = np.asarray([1,2,3,4])

np.convolve(x, np.conj(x)[::-1])

# output: array([4, 11, 20, 30, 20, 11, 4])

 

 

 

پیشنهاد: فیلم آموزشی “کانولوشن و همبستگی در دنباله‌های زمان-گسسته”‌ را مشاهده کنید. 

 

 

روش سوم: خودهمبستگی با استفاده از ماتریس توپلیتز

 

دنباله خودهمبستگی را میتوان با استفاده از ماتریس توپلیتز (Toeplitz) بدست آورد. یک نمونه مثال از ماتریس توپلیتز برای محاسبه کانولوشن در اینجا قبلا ارائه شده است. همان روش در اینجا برای محاسبه خودهمبستگی توسعه می‌یابد که در آن یک سیگنال به عنوان دنباله ورودی در نظر گرفته میشود و دیگری همان نسخه معکوس شده مزدوج اولی است. مجدد اگر دنباله اصلی حقیقی است، نیازی به عمل مزدوج کردن نیست. 

 

 

متلب

 

x = [1 2 3 4]

acf = toeplitz([x zeros(1,length(x)-1)], [x(1) zeros(1,length(conj(x))-1)])*fliplr(x).’

acf = 4 11 20 30 20 11 4

 

 

پایتون

 

import numpy as np

from scipy.linalg import toeplitz

x = np.asarray([1,2,3,4])

toeplitz(np.pad(x, (0,len(x)-1), mode=’constant’), np.pad([x[0], (0,len(x)-1), mode=’constant’))@x[::-1]

# output: array([4, 11, 20, 30, 20, 11, 4])

 

 

روش چهارم: خودهمبستگی با استفاده از تبدیل فوریه سریع

 

دنباله خودهمبستگی را میتوان با استفاده از زوج تبدیل فوریه سریع و معکوس آن، یعنی FFT/IFFT بدست آورد. در اینجا یک مثال برای محاسبه کانولوشن از این طریق آمده است. همان روش برای خودهمبستگی توسعه یافته که در آن یک سیگنال به عنوان دنباله ورودی است و دیگری نسخه معکوس شده مزدوج سیگنال اولی است. مجدد اگر سیگنال حقیقی است، نیازی به عمل مزدوج نیست.

 

 

متلب

 

x = [1 2 3 4]

L=2*length(x)-1

acf = ifft(fft(x,L).*fft(fliplr(conj(x)),L))

acf = 4 11 20 30 20 11 4

 

 

پایتون

 

import numpy as np

from scipy.fftpack import fft, ifft

x = np.asarray([1,2,3,4])

L = 2*len(x)-1

ifft(fft(x,L)*fft(np.conj(x)[::-1],L))

# output: array([4, 11, 20, 30, 20, 11, 4])

 

 

دقت کنید که در تمام موارد بالا، به دلیل خاصیت متقارن (Symmetric) بودن تابع خودهمبستگی، مقدار واقع شده در مرکز نمودار نماینده R_{xx}(0) است.

 

 

ساخت ماتریس خودهمبستگی

 

ماتریس خودهمبستگی، شکل خاصی از ماتریس است که از دنباله خودهمبستگی ساخته میشود. این ماتریس به شکل زیر است:

 

 

 

در صورت داشتن دنباله خودهمبستگی، R_{xx}(0), R_{xx}(1), ..., R_{xx}(N-1)، ماتریس خودهمبستگی به سادگی ساخته میشود. این ماتریس، R_x یک ماتریس با تقارن هرمیتی (Hermitian) و در عین حال یک ماتریس توپلیتز است. از این خواص در کدهای زیر برای ساخت ماتریس خودهمبستگی استفاده شده است:

 

 

متلب

 

   x = [1+j 2+j 3-j] % x is complex

acf = conv(x,fliplr(conj(x)))

Rxx = acf(3:end); % R_xx(0) is the center element

Rx = toeplitz(Rxx,[Rxx(1) conj(Rxx(2:end))])

 

 

پایتون

 

import numpy as np

from scipy.fftpack import fft, ifft

x = np.asarray([1+1j,2+1j,3-1j]) # x is complex

acf = np.convolve(x,np.conj(x)[::-1])

Rxx = acf[2:] # R_xx[0] is the center element

Rx = toeplitz(Rxx, np.hstack((Rxx[0], np.conj(Rxx[1:]))))

 

 

نتیجه:

    \[ R_x = \begin{bmatrix} 17 & 8+6i & 2+4i \\ 8-6i & 17 & 8+6i \\ 2-4i & 8-6i & 17 \end{bmatrix} \]

 

 

 

 

 

 

منبع: www.gaussianwaves.com