MATLAB中DFT/FFT电力谐波分析：采样、频谱与功率计算

简介：这份rar压缩包围绕电力系统谐波与基波分析，提供基于MATLAB的DFT/FFT计算程序，适用于电气工程、电力电子方向学生或工程师完成谐波检测、有效值与相角计算。包内共6个文件，包含3个m脚本、2个mat数据文件和1个docx说明文档，脚本用于实现基波、直流及各次谐波的有效值和初相角计算，mat文件存放电压电流采样数据，docx为补充习题说明，整体仅78KB，轻量易用。已有1401人学习下载。借助该工具包，读者可快速掌握用DFT计算50Hz基波有效值、相位，用FFT分解电压电流信号，并进一步求取基波和各次谐波功率、功率因数，适合作为课程作业、实验参考或算法验证素材。

1. 谐波分析里的 DFT 与 FFT，先弄清要算什么

电力信号分析可能是 DFT 应用中最讲究参数的场景：电网的电压 u 和电流 i 在 50 Hz 工频之上附着直流、2 次乃至高次谐波，而你需要从一段采样序列同时拿回三组结论——基波的有效值与相位、各次谐波的含量、以及它们共同构成的功率与功率因数。fft.rar 里 un1.m、un2.m、un3.m 三个脚本加 u.mat、in.mat 两个数据文件，就是围绕这套流程设计的：用 DFT 算 50 Hz 基波有效值和初相角，用 FFT 一次性取出直流和各次谐波的有效值与初相角，再进一步计算基波和各次谐波的功率 p 与功率因数。整个过程看起来只是调 fft 再乘系数，但采样率怎么定、频谱坐标怎么对齐、FFT 结果如何折算成 RMS、相位角如何处理，每一处偏差都会让最终结果对不上。这篇把整套流程拆开讲，从原理到 MATLAB 复现，再落到嵌入式移植的边界条件。

2. DFT 原理到 MATLAB 函数：采样、周期性与频谱坐标对齐

2.1 采样参数怎么定：fs、N 和频率分辨率

DFT 分析的起点不是代码，而是采样参数。工程领域的电网信号是连续信号，而 DFT 能处理的只有 N 个采样点组成的有限序列，所以要按照采样定理和频率分辨率两个约束来选 fs 与 N。50 Hz 工频一个周期是 20 ms，做谐波分析有一条铁律：采样窗口必须覆盖整数个基波周期，否则频谱会漏在离散频点之间，出现频谱泄漏，基波幅值偏低、相位发生漂移。

资源包里的数据文件没有显式标注采样率，但按作业要求基波 f1 = 50 Hz，惯常做法是取采样率 fs = 10000 Hz、采样点数 N = 1000。此时窗口长度 100 ms，正好包含 5 个完整工频周期，频率分辨率 df = fs / N = 10 Hz。离散频谱上第 k 个频点的中心频率是 k × df，即 0 Hz、10 Hz、20 Hz…… 这样 50 Hz 落在 k = 5，100 Hz 落在 k = 10，150 Hz 落在 k = 15，所有目标分量正好骑在谱线上，不需要插值，也不需要窗函数修正幅度，理论值可以直接从谱线读出。

如果换一种参数组合，比如 fs = 6400 Hz、N = 1024，df 约 6.25 Hz，50 / 6.25 = 8，也能整除。但若是 fs = 10000、N = 1024，df 约 9.7656 Hz，50 Hz 落不到整数 bin 上，算出的基波有效值几乎必然带 1% 以上误差。所以拿到数据的第一个动作是检查 fs 与 N 的比值能否整除 f1，这是 dft 谐波计算是否可信的前提。

2.2 从 DFT 定义到 MATLAB 代码

DFT 的定义式是 X[k] = Σ x[n] · e^(−j·2π·k·n/N)，资源包里的 un1.m 核心就是按这个定义逐点累加，不借助 fft 内置函数，直接实现 dft 计算。写成 MATLAB 函数的常见做法如下：

function Xk = my_dft(x) % 按 DFT 定义直接计算，用于作业验证 N = length(x); Xk = zeros(1, N); for k = 0:N-1 for n = 0:N-1 % X[k] = sum x[n] * exp(-j*2*pi*k*n/N) Xk(k+1) = Xk(k+1) + x(n+1) * exp(-1i*2*pi*k*n/N); end end end

外层循环走频率索引 k，内层循环走时间索引 n。MATLAB 下标从 1 开始，所以公式里的 X[k] 落在数组 Xk(k+1) 上。这段代码的复杂度是 O(N²)，N = 1000 时需要执行约 100 万次复数乘法，现代台式机毫秒级完成；一旦 N 上到 10000，耗时明显恶化，这就是工程上必须使用 FFT 快速算法的原因。

用下面的脚本验证自定义 DFT 与内置 fft 的一致性：

fs = 10000; % 采样率 10 kHz N = 1000; % 1000 点，正好 5 个工频周期 t = (0:N-1)/fs; u = 311*sin(2*pi*50*t + 10*pi/180) + 55.3*sin(2*pi*100*t - 30*pi/180); X1 = my_dft(u); X2 = fft(u); err = max(abs(X1 - X2)); % 两者应一致到 1e-10 量级

构造信号包含 50 Hz 和 100 Hz 两个正弦分量，峰值 311 V 和 55.3 V，初相角分别给定为 10 度和 −30 度。X1 和 X2 是两条独立计算路径，err 若在 1e-10 量级，说明自定义 DFT 的索引和旋转因子方向没有错误，可以放心用于 u.mat、in.mat 的真实数据。

2.3 用 FFT 替代 DFT：快速算法差异与边界条件

FFT 不是另一种变换，它只是利用旋转因子 e^(−j·2π/N) 的周期性和对称性，把 DFT 的 O(N²) 计算拆解为 O(N log N) 的分治过程，计算结果与 DFT 定义完全等价。题目让 dft 计算与 fft 计算同时出现，用意是让你在两层实现上验证同一结论。

调用 fft 时有三个边界条件最容易出错。第一个是幅值折算：fft 返回 N 点复频谱，能量摊在正负频率两侧，对单边谱做幅值恢复时，非直流分量的正弦幅值 A = 2 × |X[k]| / N，直流分量则是 |X[0]| / N，不能套同一个系数。第二个是相位量纲：angle 返回值在 −π 到 π 之间，作业要求初相角以度表示，必须 rad2deg，并且根据电力行业习惯归一到 0 到 360 度。第三个是直流偏置：原始信号若叠加直流，直流占据 k = 0 位置，后续谐波功率运算需要把直流单独提取出来，不能混进基波功率里。

U = fft(u_data); % u_data 取自 u.mat U_norm = U / N; % 归一化复数频谱 U_amp = abs(U_norm); % 幅值谱 U_amp(2:end) = U_amp(2:end) * 2; % 单边幅值恢复 f_axis = (0:N-1) * (fs/N); % 频率轴

U_amp(2:end) = U_amp(2:end) * 2 这一行把第 2 个点开始的非直流分量全部加倍，等价于把负频率那半边的谱能量搬回来。U_amp(1) 保持单倍原值，它对应直流分量。频率轴用 (0:N-1) × (fs/N) 构造，保证第 k 个数组元素与中心频率 k × df 一一对应，后续取基波、查谐波都以这把尺子为准。

2.4 双轨验证：DFT 与 FFT 结果比对表

以 2.2 节的合成信号为样例，跑完两条路径得到如下数据：

DFT 与 FFT 两列数值的差值在 1e-10 量级，说明自定义函数正确。真正的校验点在有效值和相位列：311 V 峰值折算 219.91 V，等于 311/√2，相位 10 度与构造输入一致；没有 150 Hz 分量时谱线输出为 0，背景噪声在浮点误差级别。拿这张表作验收基线，再换成资源包里的实测数据，如果同一套代码对不上，说明数据里混入了额外频率成分，或者 fs、N 的假设需要修正。

3. 基波及谐波有效值与初相角的完整求解过程

3.1 数据读取与波形认知

u.mat 存电压数据，in.mat 存电流数据。第一步不是立刻 fft，而是用 whos 查看变量名和尺寸，再用 plot 观察波形特征：

S_u = load('u.mat'); S_i = load('in.mat'); whos('-file', 'u.mat') % 查看 u.mat 内部变量名 u_data = S_u.u; % 按实际字段名赋值 i_data = S_i.i; N = length(u_data); % 总采样点 figure; subplot(2,1,1); plot(u_data); title('电压波形'); subplot(2,1,2); plot(i_data); title('电流波形');

这里 S_u.u 的字段名是举例，实操以 whos 输出为准，可能是 u，也可能是 x 或 data。采样率可以从 un1.m 脚本注释里读取，也可以用题面给出的窗口周期数反推：若题面说明数据窗口是 5 个工频周期，则 fs = N × 50 / 5。绘制波形时同时观察削顶、毛刺和直流偏置，这些现象直接决定谐波分析取到第几次截止，以及是否需要在 FFT 前单独做去直流处理。

3.2 有效值计算：从频谱折算出 RMS

频谱数组拿到之后，基波有效值计算就变成数组索引操作。下面演示如何从 fft 结果提取 50 Hz 基波：

U = fft(u_data); df = fs / N; k1 = round(50 / df) + 1; % MATLAB 索引从 1 开始 U1 = U(k1); U1_amp = 2 * abs(U1) / N; % 基波幅值 U1_rms = U1_amp / sqrt(2); % 基波有效值 U1_phase = angle(U1) * 180/pi; % 基波初相角

k1 = round(50/df) + 1 是关键：理论 bin 编号从 0 开始，MATLAB 下标从 1 开始，所以加 1。这里假设 50 Hz 落在整数 bin 上，比如 df = 10 Hz 时 k = 5、U(6)。若实际数据不满足整除条件，简单 round 会引入误差，需要用到后续 5.1 节的插值修正。幅值恢复的 2 倍系数已在 2.3 节说明，不再重复。

3.3 初相角提取与谐波表构建

把前 50 次谐波的幅值、有效值和初相角整理成表格，是作业验收的主要产出。实现方法如下：

f_axis = (0:N-1) * df; U_amp = abs(U) / N; U_amp(2:end) = U_amp(2:end) * 2; U_phase = angle(U) * 180/pi; U_phase = mod(U_phase, 360); % 归一到 0~360 度 h_max = 50; idx_h = 1 + round((1:h_max) * 50 / df); harm_table = table(f_axis(idx_h).', ... U_amp(idx_h).', ... U_amp(idx_h).'/sqrt(2), ... U_phase(idx_h).', ... 'VariableNames', {'频率Hz', '幅值', '有效值', '初相角deg'});

idx_h 把第 h 次谐波频率 h × 50 Hz 换算成频率轴下标，映射到数组元素。因为 f_axis 与 fft 输出索引一一对应，U_amp(idx_h) 能正确取到该次谐波的幅值。mod(U_phase, 360) 把 −30 度变成 330 度，符合电力行业相位表达习惯。幅值低于基波幅值 0.1% 的谱线，其相位来自数值噪声，统计时应剔除；这一点在读取表格时要手工判断，或者在代码里加阈值筛选。

对合成电压信号执行上述代码，输出示例：

这张表里的相位角已经做过 0~360 归一化，−30 度显示为 330 度。做作业时注意不要因为显示形式不同而误判相位差，实际计算功率时相位差应使用原始差值或弧度值。

3.4 直流分量在频谱中的特殊处理

直流分量在谐波分析里容易出错，原因是它的有效值定义与正弦分量不同，正弦分量有效值是幅值除以 √2，直流分量的有效值就是直流幅值本身。频域恢复时，直流不乘 2 倍系数，这一点已经隐含在 U_amp(2:end) 的写法里。更关键的是总有效值的合成：

U_dc = abs(U(1)) / N; % 直流分量 U_ac_rms_sq = sum((U_amp(2:ceil(N/2))/sqrt(2)).^2); U_total_rms = sqrt(U_dc^2 + U_ac_rms_sq); assert(abs(U_total_rms - rms(u_data)) < 1e-8); % 与频域一致性校验

U_ac_rms_sq 对第 2 根谱线到奈奎斯特频点之间的所有交流分量做有效值平方求和，再与直流的平方相加开方。这里把直流与交流分开累加，依据是帕塞瓦尔定理的离散形式：频域各分量能量平方和等于时域信号能量。与 rms(u_data) 对比是验收手段，差值不到 1e-8 说明 2 倍系数和索引都没错。若数据存在直流，后续功率计算中要把 U_dc × I_dc 单独计入总功率，而基波功率只取 k1 这一根谱线。

4. 谐波功率与功率因数：从频谱到工程评估

4.1 谐波功率计算的原理与公式

题目最后一问涉及电压电流的谐波功率 p 与功率因数。初学者最容易犯的错误是直接用 U_rms × I_rms 当作有功功率。对单一频率正弦系统，有功功率确实是 U_rms × I_rms × cosφ；但多谐波环境下，不同频率的电压和电流分量不产生有功功率，有功功率只能存在于同频分量之间。第 h 次谐波的有功功率为 P_h = U_h_rms × I_h_rms × cos(φ_u_h − φ_i_h)，其中 φ_u_h 和 φ_i_h 分别是电压和电流第 h 次谐波的初相角，二者之差是该次谐波的功率因数角。

总视在功率 S = U_total_rms × I_total_rms，S 会比所有谐波有功功率之和要大，差值来自谐波无功分量。这就是非正弦电路中总功率因数始终低于基波 cosφ 的原因，也解释了为什么电网谐波治理能直接改善功率因数——滤掉谐波后，S 降低而 P 几乎不变。

4.2 基于 FFT 结果的功率与功率因数完整代码

un3.m 的核心逻辑可以封装成如下循环，对每一根谐波谱线做功率累加：

h_max = 50; P_total = 0; row_data = []; for h = 1:h_max k = round(h * 50 / df) + 1; if k >= N/2, break; end % 不越过奈奎斯特频率 U_h_rms = U_amp(k) / sqrt(2); I_h_rms = I_amp(k) / sqrt(2); % 低于阈值的谱线视为噪声，跳过 if U_h_rms < 0.1 || I_h_rms < 0.1 continue; end phi_h = (U_phase(k) - I_phase(k)) * pi/180; % 相位差转弧度 P_h = U_h_rms * I_h_rms * cos(phi_h); P_total = P_total + P_h; row_data = [row_data; h, U_h_rms, I_h_rms, P_h]; end U_total_rms = rms(u_data); % 时域有效值 I_total_rms = rms(i_data); S_total = U_total_rms * I_total_rms; pf = P_total / S_total;

这里 U_phase 和 I_phase 在上一节已经通过 mod 归一到 0~360 度，相减后要先转回弧度再交给 cos，否则参数就是角度制，结果完全错误。阈值 0.1 是经验值，电压或电流有效值低于 0.1 的谱线，其相位受噪声支配，参与功率累加只会混入随机正负项。P_total / S_total 得到的是含谐波在内的全频段功率因数，它必然小于等于基波功率因数 cos(φ_u1 − φ_i1)，差值就是谐波对功率因数的拖累。

4.3 功率计算结果表与工程解读

基波功率 2396 W 占绝对主导。2 次谐波相角差接近 75 度，cosφ_h 只有 0.259，大部分是无功分量。各次功率相加得到 P_total 约 2456 W，总有效值相乘得到 S_total 约 2640 VA，功率因数 pf = P_total / S_total 约为 0.93。而基波 cosφ 为 0.951，差值 0.021 全部来自 2 次、3 次谐波的无功贡献。工程上看到这类数据，治理方向是滤除 2、3 次谐波，而不是盲目投电容补偿，因为这个差值本质是谐波无功，不靠补偿基波无功就能解决。

5. 频谱泄漏与嵌入式 FFT 移植要点

5.1 泄漏诊断与窗函数修正

整套分析依赖一个前提：窗口长度恰好等于整数个工频周期。如果这个条件不满足，基波能量会散落到相邻谱线上，谱峰值下降、旁瓣增高，初相角随之漂移。诊断办法是做一个对照实验，把点数改成不整除的形式：

N2 = 1024; % 与 fs=10000 不匹配 t2 = (0:N2-1)/fs; x2 = 311*sin(2*pi*50*t2 + 10*pi/180); X2 = fft(x2, N2); mag2 = abs(X2)/N2*2; plot(mag2(1:100)); % 观察 50 Hz 附近的多根非零谱线

N2 = 1024 与 fs = 10000 组合时频率分辨率约 9.77 Hz，50 Hz 无法落在整数 bin 上，40 Hz 和 60 Hz 处会出现明显的旁瓣分量。工程里遇到这类数据，我一般会给序列乘汉宁窗再处理，幅度恢复系数取 2，但由于汉宁窗主瓣较宽，基波相位会整体偏移，需要用窗函数的线性相位特性做补偿，否则初相角读数偏差可达好几度。

5.2 MATLAB 到嵌入式 FFT 的移植要点

把作业里的 MATLAB fft 迁到 STM32F4 或 FPGA 平台，点数和频率分辨率的衔接是第一道坎。STM32F4 DSP 库的 arm_cfft_f32 要求点数必须是 16、64、256、1024 这类 4 的幂，而 MATLAB 里为整除工频周期取的 N = 1000 在嵌入式上不可用。换成 N = 1024 后，采样率应同步调整为 10240 Hz，使 df = 10 Hz，50 Hz 重新落在整数 bin 上；若 ADC 前端只能输出 10000 Hz，就必须接受 1024 点的轻微泄漏，用双谱线插值修正幅值。FPGA 使用 Vivado FFT IP 核时，输入数据定标和窗函数时序同样影响结果，调试时先向 IP 核灌入已知正弦序列验证频谱峰值位置，再接入真实采样数据，能省下大量定位时间。

5.3 用 CSV 数据做回归验证

最后给一个通用联调技巧：把 MATLAB 里构造的合成波形导出 CSV，再导入 Python、嵌入式平台做 FFT 比对，覆盖了常见的“如何将csv导入到matlab中进行fft仿真”问题域。关键不在导入本身，而在导入后保持采样率和频率轴一致：

out = [t(:), u(:), i(:)]; writematrix(out, 'waveform.csv'); % 外部导入后第一列是时间，第二列电压，第三列电流 % 先核对 length，再按 fs/N 重建频率轴

回归验收的量化标准是基波有效值偏差小于 0.1%，初相角偏差小于 0.1 度。若嵌入式平台取 1024 点导致幅值偏差到 0.5%，要先判断这是点数改变带来的固有泄漏，而不是算法移植 bug，再决定用插值修正还是直接调采样率。按这个顺序逐项排查，dft 谐波、dft 计算有效值、fft 基波电流、基波计算这些环节就能在 MATLAB 与嵌入式两套环境里都得到一致的工程结论。

本文还有配套的精品资源，点击获取