Given two 2d arrays I want to calculate the cross-correlation similar as done in scipy.signal.correlate2d but instead of using a fill_value=0 I want to calculate the Pearson R from the data overlap of the two input arrays. This will lead to higher values further away from the center of the data. If the data shape is rows*cols and the correlation is calculated for a row lag of row_lag and a col lag of col_lag, then the data overlap is (rows-|row_lag|) * (cols-|col_lag|). Of course, a simple approach would be to write two nested for loops to calculate the lagged data a_lagged = a[row_lag:, col_lag:], b_lagged = b [:-row_lag, :-col_lag] and then calculate the Pearson R of a_lagged, b_lagged. However this implementation is very slow. So my approach was to first use the scipy implementation and then adjust the result by a factor that accounts for the different data size caused by lagging the inputs. def zscore(x): return (x-np.mean(x))/np.std(x) def corr2d_adjusted(a,b): """ calculates 2d cross-correlation for a,b adjusting the value without padding but normalization to data overlap size """ assert a.shape==b.shape assert a.ndim==2 # calculate data overlap for each lag used in correlation to normalize rows,cols = a.shape #2*rows-1, 2*cols-1 # total size of correlation result # correlation is performed by these row & col lags (minimum 1 row and col) row_lags = np.arange(-rows+1,rows) col_lags = np.arange(-cols+1,cols) #len(row_lags), len(col_lags) # the data overlap size is calculated by substracting the lag from the number of rows/cols rows_lagged = rows - np.abs(row_lags) cols_lagged = cols - np.abs(col_lags) data_size = rows_lagged.reshape((-1,1)) * cols_lagged.reshape((1,-1)) #data_size.shape # same as corr # calculate correlation using scipy a_ = zscore(a) # zscore both inputs b_ = zscore(b) corr = correlate2d(a_,b_) corr_adj = corr/data_size # adjust result: divide by overlap data size used for correlation calculation return corr_adj However this does not perform a zscore calculation on each lag but only globally. Note that if data_size was a constant of a.size (which is the case only in the center), then this would be equivalent to pearson correlation. Is there an idea for a fast implementation? Continue reading...
