|
Tôi nghĩ rằng có 2 điều làm tăng thêm sự nhầm lẫn cho chủ đề này:
- định nghĩa thống kê so với xử lý tín hiệu: như những người khác đã chỉ ra, trong thống kê, chúng tôi chuẩn hóa tự động tương quan thành [-1,1].
- phương sai / trung bình một phần so với không một phần: khi thời gian thay đổi ở độ trễ> 0, kích thước chồng chéo của chúng sẽ luôn <�độ dài ban đầu. Chúng ta có sử dụng giá trị trung bình và điểm chuẩn của giá trị gốc (không một phần) hay luôn tính toán giá trị trung bình và điểm trung bình mới bằng cách sử dụng chồng chéo luôn thay đổi (một phần) tạo ra sự khác biệt. (Có lẽ có một thuật ngữ chính thức cho điều này, nhưng tôi sẽ sử dụng "một phần" bây giờ).
Tôi đã tạo 5 hàm tính toán tự động tương quan của một mảng 1d, với sự khác biệt một phần và không một phần. Một số sử dụng công thức từ thống kê, một số sử dụng tương quan theo nghĩa xử lý tín hiệu, cũng có thể được thực hiện thông qua FFT. Nhưng tất cả các kết quả đều là tương quan tự động trong định nghĩa thống kê , vì vậy chúng minh họa cách chúng được liên kết với nhau. Mã bên dưới:
import numpy import matplotlib.pyplot as plt def autocorr1(x,lags): '''numpy.corrcoef, partial''' corr=[1. if l==0 else numpy.corrcoef(x[l:],x[:-l])[0][1] for l in lags] return numpy.array(corr) def autocorr2(x,lags): '''manualy compute, non partial''' mean=numpy.mean(x) var=numpy.var(x) xp=x-mean corr=[1. if l==0 else numpy.sum(xp[l:]*xp[:-l])/len(x)/var for l in lags] return numpy.array(corr) def autocorr3(x,lags): '''fft, pad 0s, non partial''' n=len(x) # pad 0s to 2n-1 ext_size=2*n-1 # nearest power of 2 fsize=2**numpy.ceil(numpy.log2(ext_size)).astype('int') xp=x-numpy.mean(x) var=numpy.var(x) # do fft and ifft cf=numpy.fft.fft(xp,fsize) sf=cf.conjugate()*cf corr=numpy.fft.ifft(sf).real corr=corr/var/n return corr[:len(lags)] def autocorr4(x,lags): '''fft, don't pad 0s, non partial''' mean=x.mean() var=numpy.var(x) xp=x-mean cf=numpy.fft.fft(xp) sf=cf.conjugate()*cf corr=numpy.fft.ifft(sf).real/var/len(x) return corr[:len(lags)] def autocorr5(x,lags): '''numpy.correlate, non partial''' mean=x.mean() var=numpy.var(x) xp=x-mean corr=numpy.correlate(xp,xp,'full')[len(x)-1:]/var/len(x) return corr[:len(lags)] if __name__=='__main__': y=[28,28,26,19,16,24,26,24,24,29,29,27,31,26,38,23,13,14,28,19,19,\ 17,22,2,4,5,7,8,14,14,23] y=numpy.array(y).astype('float') lags=range(15) fig,ax=plt.subplots() for funcii, labelii in zip([autocorr1, autocorr2, autocorr3, autocorr4, autocorr5], ['np.corrcoef, partial', 'manual, non-partial', 'fft, pad 0s, non-partial', 'fft, no padding, non-partial', 'np.correlate, non-partial']): cii=funcii(y,lags) print(labelii) print(cii) ax.plot(lags,cii,label=labelii) ax.set_xlabel('lag') ax.set_ylabel('correlation coefficient') ax.legend() plt.show()
Đây là con số đầu ra:
Chúng ta không nhìn thấy tất cả 5 dòng vì 3 trong số chúng chồng lên nhau (ở màu tím). Các phần trùng lặp là tất cả các tương quan tự động không từng phần. Điều này là do các tính toán từ các phương pháp xử lý tín hiệu ( np.correlate, FFT) không tính toán giá trị trung bình / std khác nhau cho mỗi chồng chéo.
Cũng lưu ý rằng fft, no padding, non-partialkết quả (đường màu đỏ) là khác nhau, bởi vì nó không đệm thời gian bằng 0 trước khi thực hiện FFT, vì vậy nó là FFT tròn. Tôi không thể giải thích chi tiết tại sao, đó là những gì tôi học được từ những nơi khác.
21 hữu ích 0 bình luận chia sẻ
numpy.correlate(a, v, mode='valid')[source]#
Cross-correlation of two 1-dimensional sequences.
This function computes the correlation as generally defined in signal processing texts:
\[c_k = \sum_n a_{n+k} \cdot \overline{v_n}\]
with a and v sequences being zero-padded where necessary and \(\overline x\) denoting complex conjugation. Parametersa, varray_like
Input sequences. mode{‘valid’, ‘same’, ‘full’}, optional
Refer to the convolve docstring. Note that the default is ‘valid’, unlike convolve, which uses ‘full’. old_behaviorbool
old_behavior was removed in NumPy 1.10. If you need the old behavior, use multiarray.correlate.
Returnsoutndarray
Discrete cross-correlation of a and v.
See also convolve
Discrete, linear convolution of two one-dimensional sequences. multiarray.correlate
Old, no conjugate, version of correlate. scipy.signal.correlate
uses FFT which has superior performance on large arrays.
Notes
The definition of correlation above is not unique and sometimes correlation may be defined differently. Another common definition is:
\[c'_k = \sum_n a_{n} \cdot \overline{v_{n+k}}\]
which is related to \(c_k\) by \(c'_k = c_{-k}\).
numpy.correlate may perform slowly in large arrays (i.e. n = 1e5) because it does not use the FFT to compute the convolution; in that case, scipy.signal.correlate might be preferable.
Examples
>>> np.correlate([1, 2, 3], [0, 1, 0.5])
array([3.5])
>>> np.correlate([1, 2, 3], [0, 1, 0.5], "same")
array([2. , 3.5, 3. ])
>>> np.correlate([1, 2, 3], [0, 1, 0.5], "full")
array([0.5, 2. , 3.5, 3. , 0. ])
Using complex sequences:
>>> np.correlate([1+1j, 2, 3-1j], [0, 1, 0.5j], 'full')
array([ 0.5-0.5j, 1.0+0.j , 1.5-1.5j, 3.0-1.j , 0.0+0.j ])
Note that you get the time reversed, complex conjugated result (\(\overline{c_{-k}}\)) when the two input sequences a and v change places:
>>> np.correlate([0, 1, 0.5j], [1+1j, 2, 3-1j], 'full')
array([ 0.0+0.j , 3.0+1.j , 1.5+1.5j, 1.0+0.j , 0.5+0.5j])
|
