网址:https://ww2.mathworks.cn/help/signal/ug/cross-correlation-of-two-moving-average-processes.html
描述:本案例由1个示例构成
-
针对以上案例,采用Python语言实现。
-
该示例将样本互相关与理论互相关进行比较。使用两个不同的移动平均滤波器过滤N(0,1)个白噪声输入。绘制样本和理论互相关序列。 创建一个$N(0,1)$白噪声序列。将随机数生成器设置为可再现结果的默认设置。创建两个移动平均值过滤器。一个滤波器具有脉冲响应$δ(n)+δ(n−1)$. 另一个滤波器具有脉冲响应$δ(n)−δ(n−1)$.
python
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
初始化数据
python
rng = np.random.default_rng()
xn = np.random.randn(100,1)
x = signal.lfilter([1,1], 1, xn)
x = np.array(x).flatten()
y = signal.lfilter([1,-1], 1, xn)
y = np.array(y).flatten()
定义函数detrend_none
返回x:不删除趋势
参数
----------
x : 任意对象
axis : 整数
另请参见
--------
detrend_mean : 其他detrend算法
detrend_linear : 其他detrend算法
detrend :所有detrend算法的包装
python
def detrend_none(x, axis=None):
return x
定义函数M_xcorr
绘制x和y之间的互相关
参数
----------
x, y : 长度为n的类数组
maxlags : 整数, 默认值: 10
要显示的滞后数。如果没有,将返回所有 ``2 * len(x) - 1``
滞后
Returns
-------
lags : 矩阵 (length 2*maxlags+1)
滞后向量
c : 矩阵 (length 2*maxlags+1)
滞后向量
python
def M_xcorr( x, y, normed=True, detrend=detrend_none,
maxlags=10, **kwargs):
Nx = len(x)
if Nx != len(y):
raise ValueError('x and y must be equal length')
x = detrend(np.asarray(x))
y = detrend(np.asarray(y))
correls = np.correlate(x, y, mode="full")
if normed:
correls /= np.sqrt(np.dot(x, x) * np.dot(y, y))
if maxlags is None:
maxlags = Nx - 1
if maxlags >= Nx or maxlags < 1:
raise ValueError('maxlags must be None or strictly '
'positive < %d' % Nx)
lags = np.arange(-maxlags, maxlags + 1)
correls = correls[Nx - 1 - maxlags:Nx + maxlags]
return correls, lags
python
[xc,lags] = M_xcorr(x, y, maxlags=20)
Xc = np.zeros(np.size(xc))
Xc[19] = -1
Xc[21] = 1
plt.stem(lags,xc, label = '$Sample cross-correlation$')
markerline, stemlines, baseline = plt.stem(lags,Xc, linefmt = 'r-', label = '$Theoretical autocorrelation$')
plt.setp(stemlines, 'linewidth', 2)
plt.grid(True)
plt.legend(loc = "upper left")
plt.axhline(0, color='blue', lw=2)
plt.show()
python