一、概念

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 之间：

• 1 表示完美的正相关（完全同步）。
• -1 表示完全的负相关（反向同步）。
• 0 表示没有相关性。

三、示例 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 autocorrelationx = 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_datadef 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_lagmax_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 = 1000min_lag = 200plot_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 2return 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 2return 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;
}