代码收藏家 常见算法——自相关的含义及Python、C实现
常见算法——自相关的含义及C实现
一、概念
1. 自相关概念
自相关（Autocorrelation）计算是一种用于分析时间序列数据中的模式和周期性的方法。它通过计算当前序列和其自身在不同滞后时间上的相关性，以找出数据是否在某个时间间隔内呈现重复的周期性行为。
如果自相关值较高，说明数据在这些滞后时间点上有相似的变化趋势，可能存在周期性。
如果自相关值较低，说明数据在这些滞后点之间没有显著的相关性。
简单例子：假设我们有一组表示水位的时间序列数据。如果水位是周期性变化的，那么当前时刻的水位值与前几个时刻的水位值会有某种关联，这种关联就可以通过自相关来捕捉。
2. 滞后期
滞后期（lag） 是指在计算自相关时，当前数据序列与它自身的一个滞后版本进行比较的时间差。滞后期代表两个序列之间的位移。具体来说，滞后期为 lag 的自相关计算的是数据序列中每个点与其前 lag 个数据点的相似性。
滞后期的选择通常根据数据的周期性特征来决定。例如，若我们怀疑数据有周期性波动（如正弦波），我们会尝试多个滞后期，看看在哪个滞后期的自相关值最大，进而推测数据的周期。
在常见的分析中，滞后期可能从 1 开始增加，直到找到显著的自相关值或达到数据长度的一半。
示例说明：
假设有一组有规律变化的水位数据：
waterLevel = {10, 12, 14, 16, 18, 16, 14, 12, 10, 8, 6, 8, 10, 12, 14, 16, 18, 16, 14, 12}
这种数据呈现出周期为 10 的波动。当滞后期为 10 时，自相关值会非常高，因为第 i 个数据和第 i-10 个数据会很相似。
二、自相关的计算步骤：
1. 确定滞后期 (Lag)：
自相关的滞后期是指数据点之间的间隔。举例来说，如果滞后期为 3，则我们将当前水位数据与 3 个时间步长前的水位数据进行比较。
2. 计算平均值：
在计算自相关时，通常需要减去数据的均值，来去除掉整体趋势的影响，关注数据的波动。
3. 计算自相关：
在不同的滞后期，计算当前数据和滞后数据的相关性。自相关值的范围在 -1 到 1 之间：
三、示例 Python自相关计算
本示例将使用Python生成 1000*10 个数据点，数据点呈正弦分布，但加入了随机参数。
然后程序计算在lag=994的时候自相关度；
最后程序分别计算了200-1000的时候自相关度变化趋势，以方便找出自相关度最大的时候的lag值。
1. 代码
这些点略呈正弦规律，但添加大量噪声：
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# 模拟生成 1000*10 的水位数据 (略有规律 + 随机波动 + 周期噪声)
def generate_data(num_points=10000):
waterLevelData = []
base_frequency = 2 * np.pi / 1000 # 基础周期
base_water_level = 27.5 # 基础水位
amplitude = 12.5 # 振幅
for i in range(num_points):
# 每个点的周期有小幅度的随机波动
noisy_frequency = base_frequency * (1 + np.random.normal(0, 0.01)) # 1% 随机噪声
water_level = base_water_level + amplitude * np.sin(noisy_frequency * i) + np.random.normal(0, 1) # 随机波动的水位值
waterLevelData.append(water_level)
return waterLevelData
# 绘制水位波形图
def plot_waveform(waterLevelData):
time = np.arange(0, len(waterLevelData))
# 创建图形和轴
plt.figure(figsize=(15, 6))
# 绘制水位随时间变化的曲线
plt.plot(time, waterLevelData, label="Water Level", color='blue', marker=None, linewidth=0.5)
# 添加图形细节
plt.title("Water Level Over Time with Random Fluctuations and Noisy Periodicity")
plt.xlabel("Time")
plt.ylabel("Water Level")
plt.legend()
# 显示图形
plt.grid(True)
plt.show()
def autocorrelation(data, lag=1):
n = len(data)
mean = np.mean(data)
var = np.var(data)
data = np.array(data)
# Compute autocorrelation
x = data[:-lag]
y = data[lag:]
autocorr = np.correlate(x - mean, y - mean, mode='valid') / ((n - lag) * var)
return autocorr[0]
def plot_autocorrelation(data, min_lag, max_lag):
autocorrs = [autocorrelation(data, lag) for lag in range(min_lag, max_lag+1)]
plt.figure(figsize=(10, 5))
plt.stem(range(min_lag, max_lag+1), autocorrs, use_line_collection=True)
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title('Autocorrelation function (ACF)')
plt.grid(True)
plt.show()
def resample_data(data, num_points):
resampled_data = signal.resample(data, num_points)
return resampled_data
def find_max_autocorrelation_lag_and_value(data, min_lag, max_lag):
autocorrs = [autocorrelation(data, lag) for lag in range(min_lag, max_lag+1)]
max_lag = np.argmax(autocorrs) + min_lag
max_value = np.max(autocorrs)
return max_lag, max_value
# 主函数
if __name__ == "__main__":
waterLevelData = generate_data()
plot_waveform(waterLevelData)
# 计算指定lag的自相关性
autocorr = autocorrelation(waterLevelData, 994)
print("Autocorrelation: ", autocorr)
# 绘制自相关性图
max_lag = 1000
min_lag = 200
plot_autocorrelation(waterLevelData, min_lag, max_lag)
# 找到最大自相关性的lag和值
max_autocorr_lag, max_autocorr_value = find_max_autocorrelation_lag_and_value(waterLevelData, min_lag, max_lag)
print("Lag with maximum autocorrelation: ", max_autocorr_lag)
print("Maximum autocorrelation: ", max_autocorr_value)
2. 运行结果
原始数据图形：
运行结果：
Autocorrelation: 0.8718133244362108
Lag with maximum autocorrelation: 1000
Maximum autocorrelation: 0.8790579058186793
相关度趋势：
当数据周期比较大时，使用较小的lag将可能计算出错误的自相关度。本示例中数据周期近1000，查找最佳lag时不能从1开始找。。
四、C语言实现自相关算法
1. 代码
下面将实现一个功能，用于检测传入的水位数据是否呈现规律变化。每次传入一个新的水位数据，缓存在一个固定大小的环形缓冲区。当缓冲区满了以后，开始计算自相关。
每次计算自相关可以检测自相关值是否显示出明显的周期性。如果自相关值较高，则输出对应位置索引和水位值，及自相关值。
输入数据模拟的是一段无规律值、一段有规律值、再来一段无规律值、一段固定值，波形如下所示：
代码：
#include <stdio.h>
#include <math.h>
#define bufferSize 2000 // 定义水位记录缓冲区大小
#define autocorrThreshold 0.9 // 自相关阈值，用于判断是否有规律
int waterLevelBuffer[bufferSize]; // 水位数据缓冲区
int currentIndex = 0; // 当前缓冲区索引
int dataCount = 0; // 当前已传入的数据计数
// 计算滞后期为 lag 的自相关值
float computeAutocorrelation(int lag) {
float mean = 0;
float autocorrelation = 0;
float variance = 0;
// 计算水位数据的平均值
for (int i = 0; i < bufferSize; i++) {
mean += waterLevelBuffer[i];
}
mean /= bufferSize;
// 计算方差（归一化因子）
for (int i = 0; i < bufferSize; i++) {
variance += (waterLevelBuffer[i] - mean) * (waterLevelBuffer[i] - mean);
}
// 计算滞后期为 lag 的自相关值
for (int i = 0; i < bufferSize; i++) {
int index1 = (currentIndex + i) % bufferSize;
int index2 = (currentIndex + i + lag) % bufferSize;
autocorrelation += (waterLevelBuffer[index1] - mean) * (waterLevelBuffer[index2] - mean);
}
// 返回归一化后的自相关值
return autocorrelation / variance;
}
// 检查水位数据是否有规律,返回自相关度
float checkForPattern() {
if (dataCount < bufferSize) {
// 缓冲区还未满，无法计算自相关
return 0;
}
// 检测滞后期为 900 到 1100 的自相关
for (int lag = 900; lag <= 1100; lag++) {
float autocorr = computeAutocorrelation(lag);
if (autocorr > autocorrThreshold) {
// printf("detected at lag = %d, Autocorrelation = %.4f\n", lag, autocorr);
return autocorr;
}
}
// 没有检测到规律，继续收集数据
return 0;
}
// 模拟传入一个水位数据,返回自相关度
float addWaterLevelData(int newWaterLevel) {
// 将新的水位数据存入缓冲区
waterLevelBuffer[currentIndex] = newWaterLevel;
currentIndex = (currentIndex + 1) % bufferSize;
dataCount++;
// 检查当前缓冲区数据是否有规律
return checkForPattern();
}
// 模拟生成有规律变化的水位数据
int generateRegularWaterLevel(int step) {
float noise = ((float)rand() / RAND_MAX) * 8 - 4; // Generate a random float between -2 and 2
return 25 + 12.5 * sin(2 * M_PI * step / 1000)+noise; // 周期为 1000 的正弦波
}
// 模拟生成无规律变化的水位数据
int generateIrregularWaterLevel() {
return 25 + rand() % 5; // 随机生成0到20的水位
}
int main() {
// 模拟生成1000*12个水位数据，逐个传入检测
printf("Simulating 12000 water level data points...\n");
for (int i = 0; i < 1000*12; i++) {
// 模拟产生数据，可以根据需要生成有规律或无规律的数据
int waterLevel;
if (i > 1000*3 && i<1000*6) {
// 当中3000有规律值
waterLevel = generateRegularWaterLevel(i);
} else if(i>8000){
// 最后常数值
waterLevel=25;
}else{
// 其它是无规律值
waterLevel = generateIrregularWaterLevel();
}
// 传入数据
float val = addWaterLevelData(waterLevel);
if(val>0) {
// 打印当前传入的水位数据
printf("Water level at step %d: %d, autoCorrelation=%f\n", i + 1, waterLevel, val);
}
// printf(",%d",waterLevel);
}
return 0;
}
2. 运行结果：
Simulating 10000 water level data points...
Water level at step 4880: 14, autoCorrelation=0.900089
Water level at step 4881: 17, autoCorrelation=0.901121
Water level at step 4882: 14, autoCorrelation=0.901935
Water level at step 4883: 19, autoCorrelation=0.900232
Water level at step 4884: 14, autoCorrelation=0.900028
...
Water level at step 6193: 26, autoCorrelation=0.900189
Water level at step 6194: 26, autoCorrelation=0.901144
Water level at step 6195: 28, autoCorrelation=0.900636
Water level at step 6196: 28, autoCorrelation=0.900139
3. 优化
在新的数据进来时，不需要对所有数据计算均值和方差，只需要对新数据进行计算，下面是改进的算法：
#include <stdio.h>
#include <math.h>
#define bufferSize 2000 // 定义水位记录缓冲区大小
#define autocorrThreshold 0.9 // 自相关阈值，用于判断是否有规律
int waterLevelBuffer[bufferSize]; // 水位数据缓冲区
int currentIndex = 0; // 当前缓冲区索引
int dataCount = 0; // 当前已传入的数据计数
float mean = 0;
float variance = 0;
// 计算滞后期为 lag 的自相关值
float computeAutocorrelation(int lag) {
float autocorrelation = 0;
if(lag>bufferSize){
return 0;
}
// 计算滞后期为 lag 的自相关值
for (int i = 0; i < bufferSize; i++) {
int index1 = (currentIndex + i) % bufferSize;
int index2 = (currentIndex + i + lag) % bufferSize;
autocorrelation += (waterLevelBuffer[index1] - mean) * (waterLevelBuffer[index2] - mean);
}
// 返回归一化后的自相关值
return autocorrelation / variance;
}
// 检查水位数据是否有规律,返回自相关度
float checkForPattern() {
if (dataCount < bufferSize) {
// 缓冲区还未满，无法计算自相关
return 0;
}
// 检测滞后期为 900 到 1100 的自相关
for (int lag = 900; lag <= 1100; lag++) {
float autocorr = computeAutocorrelation(lag);
if (autocorr > autocorrThreshold) {
return autocorr;
}
}
// 没有检测到规律，继续收集数据
return 0;
}
// 模拟传入一个水位数据,返回自相关度
float addWaterLevelData(int newWaterLevel) {
int oldWaterLevel = waterLevelBuffer[currentIndex];
// 更新均值和方差
if(dataCount ==0){
mean = newWaterLevel;
variance = 0;
}
else if (dataCount < bufferSize) {
float oldMean = mean;
mean = (mean * dataCount + newWaterLevel) / (dataCount+1);
variance = variance + (newWaterLevel - oldMean) * (newWaterLevel - mean);
} else {
mean = mean + (newWaterLevel - oldWaterLevel) / bufferSize;
variance = variance - (oldWaterLevel - mean) * (oldWaterLevel - mean)
+ (newWaterLevel - mean) * (newWaterLevel - mean);
}
// 将新的水位数据存入缓冲区
waterLevelBuffer[currentIndex] = newWaterLevel;
if (dataCount < bufferSize-1) {
currentIndex++;
} else {
currentIndex = (currentIndex + 1) % bufferSize;
}
dataCount++;
// 检查当前缓冲区数据是否有规律
return checkForPattern();
}
// 模拟生成有规律变化的水位数据
int generateRegularWaterLevel(int step) {
float noise = ((float)rand() / RAND_MAX) * 8 - 4; // Generate a random float between -2 and 2
return 25 + 12.5 * sin(2 * M_PI * step / 1000)+noise; // 周期为 1000 的正弦波
}
// 模拟生成无规律变化的水位数据
int generateIrregularWaterLevel() {
return 25 + rand() % 5; // 随机生成0到20的水位
}
int main() {
// 模拟生成1000*12个水位数据，逐个传入检测
printf("Simulating 12000 water level data points...\n");
for (int i = 0; i < 1000*12; i++) {
// 模拟产生数据，可以根据需要生成有规律或无规律的数据
int waterLevel;
if (i > 1000*3 && i<1000*6) {
// 3000个数据是有规律的
waterLevel = generateRegularWaterLevel(i);
} else if(i>8000){
// 4000个常数值
waterLevel=25;
}else{
// 其它数据是无规律的
waterLevel = generateIrregularWaterLevel();
}
// 传入数据
float val = addWaterLevelData(waterLevel);
if(val>0) {
// 打印当前传入的水位数据
printf("Water level at step %d: %d, autoCorrelation=%f\n", i + 1, waterLevel, val);
}
// printf(",%d",waterLevel);
}
return 0;
}
4. 检测规律波动
上面的算法里，在最后一段数据是常数不变时，算法仍会返回高相关值。有时候我们希望检测到水位发生了规则变化，而不是仅高相关值。
下面采用局部标准差的滑动窗口检测，即通过检测一段数据中的局部波动来判断数据是否有足够的变化。方法是计算水位数据的局部标准差，如果局部标准差过小，说明数据可能是常值或变化不大，从而跳过自相关计算。
另外下面代码实现连续100次检测到自相关值达到目标值，则输出检测成功。
标准差计算的是方差的平方根
改进后的代码：
#include <stdio.h>
#include <math.h>
#define bufferSize 2000 // 定义水位记录缓冲区大小
#define autocorrThreshold 0.9 // 自相关阈值，用于判断是否有规律
#define alertThreshold 100 // 连续达到自相关阈值的次数，用于报警
#define windowSize 100 // 用于计算局部标准差的滑动窗口大小
int waterLevelBuffer[bufferSize]; // 水位数据缓冲区
int currentIndex = 0; // 当前缓冲区索引
int dataCount = 0; // 当前已传入的数据计数
int alertCounter = 0; // 记录连续达到自相关阈值的次数
float mean = 0;
float variance = 0;
int patternDetected = 0; // 标志，是否检测到规律变化
// 计算滞后期为 lag 的自相关值
float computeAutocorrelation(int lag) {
// 检查方差是否接近0。如果是，那么就直接返回0，表示数据没有变化，自相关值没有意义。
if (lag >= bufferSize || variance < 1e-6) { // 检查滞后期和方差
return 0;
}
float autocorrelation = 0;
// 计算滞后期为 lag 的自相关值
for (int i = 0; i < bufferSize; i++) {
int index1 = (currentIndex + i) % bufferSize;
int index2 = (currentIndex + i + lag) % bufferSize;
autocorrelation += ((float)waterLevelBuffer[index1] - mean) * ((float)waterLevelBuffer[index2] - mean);
}
// 返回归一化后的自相关值
return autocorrelation / variance;
}
// 计算滑动窗口内的局部标准差
float computeLocalStandardDeviation() {
if (dataCount < windowSize) {
return 0; // 数据不足时无法计算局部标准差
}
// 计算窗口内的均值
float localMean = 0;
for (int i = 0; i < windowSize; i++) {
int index = (currentIndex + i) % bufferSize;
localMean += (float)waterLevelBuffer[index];
}
localMean /= windowSize;
// 计算窗口内的标准差
float localVariance = 0;
for (int i = 0; i < windowSize; i++) {
int index = (currentIndex + i) % bufferSize;
localVariance += (float)pow((float)waterLevelBuffer[index] - localMean, 2);
}
localVariance /= windowSize;
return sqrt(localVariance); // 返回局部标准差
}
// 检查水位数据是否有规律，返回自相关度
float checkForPattern() {
if (dataCount < bufferSize) {
// 缓冲区还未满，无法计算自相关
return 0;
}
// 计算局部标准差，确保数据有足够的波动
float localStdDev = computeLocalStandardDeviation();
if (localStdDev < 1.0) {
return 0; // 数据变化幅度过小，跳过自相关计算
}
// 检测滞后期为 900 到 1100 的自相关
for (int lag = 900; lag <= 1100; lag++) {
float autocorr = computeAutocorrelation(lag);
if (autocorr > autocorrThreshold) {
return autocorr;
}
}
// 没有检测到规律，继续收集数据
return 0;
}
// 模拟传入一个水位数据，返回自相关度
float addWaterLevelData(int newWaterLevel) {
int oldWaterLevel = waterLevelBuffer[currentIndex];
// 更新均值和方差
if (dataCount == 0) {
mean = (float)newWaterLevel;
variance = 0;
} else if (dataCount < bufferSize) {
float oldMean = mean;
mean = (mean * (float)dataCount + (float)newWaterLevel) / ((float)dataCount + 1);
variance = variance + ((float)newWaterLevel - oldMean) * ((float)newWaterLevel - mean);
} else {
mean = mean + (float)(newWaterLevel - oldWaterLevel) / bufferSize;
variance = variance - ((float)oldWaterLevel - mean) * ((float)oldWaterLevel - mean)
+ ((float)newWaterLevel - mean) * ((float)newWaterLevel - mean);
}
// 将新的水位数据存入缓冲区
waterLevelBuffer[currentIndex] = newWaterLevel;
currentIndex = (currentIndex + 1) % bufferSize;
dataCount++;
// 检查当前缓冲区数据是否有规律
return checkForPattern();
}
// 模拟生成有规律变化的水位数据
int generateRegularWaterLevel(int step) {
float noise = (float)(rand() / (float)RAND_MAX) * 8 - 4; // 生成随机噪声
return (int)(25 + 12.5 * sin(2 * M_PI * step / 1000) + noise); // 周期为 1000 的正弦波
}
// 模拟生成无规律变化的水位数据
int generateIrregularWaterLevel() {
return (int)(25 + rand() % 5); // 随机生成 0 到 5 的水位
}
int main() {
// 模拟生成 12000 个水位数据，逐个传入检测
printf("模拟生成 12000 水位数据...\n");
for (int i = 0; i < 12000; i++) {
// 模拟产生数据，可以根据需要生成有规律或无规律的数据
int waterLevel;
if (i > 3000 && i < 6000) {
// 中间这段数据是有规律的
waterLevel = generateRegularWaterLevel(i);
} else if (i > 8000) {
waterLevel = 25; // 这一段数据是常值
} else {
// 其余部分数据是无规律的
waterLevel = generateIrregularWaterLevel();
}
// 传入数据
float val = addWaterLevelData(waterLevel);
// 检查自相关是否超过阈值，并进行计数
if (val > autocorrThreshold) {
alertCounter++;
if (alertCounter >= alertThreshold && !patternDetected) {
printf("第 %d 步连续 %d 次检测到自相关值超阈值!\n", i + 1, alertThreshold);
patternDetected = 1; // 标志已检测到规律
}
} else {
alertCounter = 0; // 自相关值未达到阈值时，重置计数器
}
}
return 0;
}
作者：编程圈子
