基于滑动窗DFT的实时THD算法解析

在电力电子控制与电能质量监测中，总谐波失真（THD）是衡量信号纯度的核心指标。传统的THD计算通常依赖于对整个数据块进行快速傅里叶变换（FFT），但这往往需要等待一个完整的周期，导致计算延迟。

本文将解析一种基于滑动窗（Sliding Window）的高效算法，通过递归更新累加和，实现对THD的实时计算。

核心原理：从FFT到滑动DFT

传统的FFT算法虽然计算效率高，但属于“批处理”模式，必须等待 $N$ 个点采集完毕才能进行计算。对于实时控制系统而言，这种等待引入了至少一个周期的滞后。

相比之下，滑动离散傅里叶变换（Sliding DFT）利用了数据的重叠特性。当一个新的采样点进入窗口时，只有一个最旧的点移出窗口。因此，我们不需要重新计算整个频谱，只需在原有的累加和基础上，减去旧数据的贡献，加上新数据的贡献即可。

本算法的核心在于维护三个关键的累加变量，从而将计算复杂度从 $O(N)$ 降低到 $O(1)$。

算法实现解析

该算法主要包含三个步骤：状态保持、滑动更新与数值计算。

1. 状态保持与初始化
为了在函数调用之间保存状态，代码使用了 persistent 变量。

• buffer: 长度为 $N$ 的循环缓冲区，存储当前窗口内的时域信号。
• sum_sq, sum_sin, sum_cos: 分别存储 $u^2$、$u\cdot \sin (\theta )$ 和 $u\cdot \cos (\theta )$ 的累加和。

2. 滑动更新逻辑
这是算法的灵魂所在。当一个新采样值 uk 到来时：

• 定位旧值：读取缓冲区当前位置 idx 的旧值 old_val（即即将被覆盖的值）。
• 更新缓冲区：将 uk 写入 buffer(idx)。
• 递归更新累加和：
  $$Su{m}_{new}=Su{m}_{old}+Contribution(uk)-Contribution(old\_val)$$
  代码中分别更新了能量累加和 sum_sq 以及基波分量的正余弦累加和 sum_sin、sum_cos。
• 索引循环：更新 idx，若到达窗口末尾则归零。

3. 数值计算与THD推导
利用更新后的累加和，我们可以直接计算均方根值和基波幅值：

• 总均方值 ($V_{rms}^2$)：
  $$V_{rms}^2=\frac{1}{N}\sum u^2=\frac{\text{sum\_sq}}{N}$$

• 基波幅值 ($V_{fund}$)：
  利用DFT定义，基波幅值为：
  $$V_{fund\_amp}=\frac{2}{N}\sqrt{{\text{sum\_sin}}^2+{\text{sum\_cos}}^2}$$

• 基波均方值 ($V_{fund\_sq}$)：
  $$V_{fund\_sq}=\frac{V_{fund\_amp}^2}{2}$$

最终，根据THD的定义（谐波含量有效值与基波有效值之比），计算公式如下：
$$THD=\sqrt{\frac{V_{rms}^2-V_{fund\_sq}}{V_{fund\_sq}}}\times 100\%$$

关键代码片段

% 核心更新逻辑：加新减旧
old_val = buffer(idx);
buffer(idx) = uk;

% 计算当前相位角
angle = 2.0 * pi * (idx - 1) / N;
sin_val = sin(angle);
cos_val = cos(angle);

% 更新积分项
sum_sq = sum_sq + (uk * uk) - (old_val * old_val);
sum_sin = sum_sin + (uk * sin_val) - (old_val * sin_val);
sum_cos = sum_cos + (uk * cos_val) - (old_val * cos_val);

总结

该算法通过巧妙的数学变换，避免了繁重的矩阵运算和内存搬移。其优势在于：

• 低延迟：每来一个采样点即可更新一次THD值，无等待时间。
• 计算量小：仅需几次乘加运算，非常适合在DSP或MCU的中断服务程序中运行。
• 内存占用低：无需存储庞大的FFT蝶形运算中间变量。

这种方法是实时嵌入式系统中进行电能质量分析的理想选择。

✅ 最终代码 (实时滑动窗版)

function thd_result = calculate_THD_Sliding(uk)
    % 输入: uk (实时标量信号)
    % 输出: thd_result (THD值)

    % --- 1. 参数配置 ---
    % 假设采样率 5kHz (Ts=0.0002), 信号 50Hz (周期 0.02s)
    % 窗口长度 N = 0.02 / 0.0002 = 100
    N = 100; 
    
    persistent buffer       % 信号缓冲区 [1 x N]
    persistent idx          % 当前写入位置 (1 到 N)
    persistent sum_sq       % 窗口内 u^2 的累加和
    persistent sum_sin      % 窗口内 u*sin 的累加和
    persistent sum_cos      % 窗口内 u*cos 的累加和
    persistent is_first     % 第一次运行标志

    % --- 2. 初始化 ---
    if isempty(is_first)
        buffer = zeros(1, N);
        idx = 1;
        sum_sq = 0;
        sum_sin = 0;
        sum_cos = 0;
        is_first = true;
    end

    % --- 3. 滑动窗更新逻辑 (核心算法) ---
    
    % A. 找出“最老”的那个点，准备把它从积分中减去
    % 因为是循环缓冲区，最老的点就在当前 idx 的位置（即将被覆盖）
    old_val = buffer(idx);
    
    % B. 更新缓冲区：写入新数据
    buffer(idx) = uk;
    
    % C. 计算当前点的角度 (用于 DFT)
    % 角度 = 2 * pi * (位置-1) / N
    angle = 2.0 * pi * (idx - 1) / N;
    sin_val = sin(angle);
    cos_val = cos(angle);
    
    % D. 更新累加和 (加新减旧)
    % 1. 更新总能量 (u^2)
    sum_sq = sum_sq + (uk * uk) - (old_val * old_val);
    
    % 2. 更新基波分量 (u * sin/cos)
    sum_sin = sum_sin + (uk * sin_val) - (old_val * sin_val);
    sum_cos = sum_cos + (uk * cos_val) - (old_val * cos_val);
    
    % E. 更新索引 (循环移位)
    if idx < N
        idx = idx + 1;
    else
        idx = 1;
    end

    % --- 4. 计算 THD ---
    
    % 计算均方根值的平方 (Vrms^2)
    % 公式: (1/N) * sum(u^2)
    V_rms_sq = sum_sq / N;
    
    % 计算基波幅值 (Vfund)
    % DFT 幅值公式: (2/N) * sqrt(sum_sin^2 + sum_cos^2)
    V_fund_amp = (2.0 / N) * sqrt(sum_sin^2 + sum_cos^2);
    
    % 基波均方值 (Vfund_sq) = (幅值 / sqrt(2))^2 = 幅值^2 / 2
    V_fund_sq = (V_fund_amp * V_fund_amp) / 2.0;
    
    % 防止除零和负数
    if V_fund_sq < 1e-10
        thd_result = 0;
    else
        numerator = V_rms_sq - V_fund_sq;
        if numerator < 0
            numerator = 0; % 消除数值误差导致的微小负数
        end
        thd_result = sqrt(numerator / V_fund_sq);
    end

end

测试50HZ 3次0.05 5次0.02 理论计算值5.29%