常见算法——自相关的含义及Python、C实现

常见算法——自相关的含义及C实现

• 一、概念

• 1. 自相关概念
• 2. 滞后期

• 示例说明：

• 二、自相关的计算步骤：

• 1. 确定滞后期 (Lag)：
• 2. 计算平均值：
• 3. 计算自相关：

• 三、示例 Python自相关计算

• 1. 代码
• 2. 运行结果

• 四、C语言实现自相关算法

• 1. 代码
• 2. 运行结果：
• 3. 优化
• 4. 检测规律波动

一、概念

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 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;
}