简介这份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 5100 Hz 落在 k 10150 Hz 落在 k 15所有目标分量正好骑在谱线上不需要插值也不需要窗函数修正幅度理论值可以直接从谱线读出。如果换一种参数组合比如 fs 6400 Hz、N 1024df 约 6.25 Hz50 / 6.25 8也能整除。但若是 fs 10000、N 1024df 约 9.7656 Hz50 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(k1) Xk(k1) x(n1) * exp(-1i*2*pi*k*n/N); end end end外层循环走频率索引 k内层循环走时间索引 n。MATLAB 下标从 1 开始所以公式里的 X[k] 落在数组 Xk(k1) 上。这段代码的复杂度是 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 幅值 (V)FFT 幅值 (V)有效值 (V)初相角 (deg)基波 50 Hz311.000311.000219.9110.002 次 100 Hz55.30055.30039.10−30.003 次 150 Hz0.0000.0000—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% 的谱线其相位来自数值噪声统计时应剔除这一点在读取表格时要手工判断或者在代码里加阈值筛选。对合成电压信号执行上述代码输出示例谐波次数频率 (Hz)幅值 (V)有效值 (V)初相角 (deg)基波 150311.00219.9110.0210055.3039.10330.0315022.1015.6360.042000.000.00—这张表里的相位角已经做过 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_rmsS 会比所有谐波有功功率之和要大差值来自谐波无功分量。这就是非正弦电路中总功率因数始终低于基波 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 功率计算结果表与工程解读谐波次数 hU_h_rms (V)I_h_rms (A)P_h (W)φ_u−φ_i (deg)cosφ_h1219.9111.032396.418.00.951239.102.2754.1−75.00.259315.631.274.9−30.00.866直流2.500.250.63——基波功率 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; % 与 fs10000 不匹配 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 Hz50 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 Hz50 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 与嵌入式两套环境里都得到一致的工程结论。本文还有配套的精品资源点击获取
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 5100 Hz 落在 k 10150 Hz 落在 k 15所有目标分量正好骑在谱线上不需要插值也不需要窗函数修正幅度理论值可以直接从谱线读出。如果换一种参数组合比如 fs 6400 Hz、N 1024df 约 6.25 Hz50 / 6.25 8也能整除。但若是 fs 10000、N 1024df 约 9.7656 Hz50 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(k1) Xk(k1) x(n1) * exp(-1i*2*pi*k*n/N); end end end外层循环走频率索引 k内层循环走时间索引 n。MATLAB 下标从 1 开始所以公式里的 X[k] 落在数组 Xk(k1) 上。这段代码的复杂度是 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 幅值 (V)FFT 幅值 (V)有效值 (V)初相角 (deg)基波 50 Hz311.000311.000219.9110.002 次 100 Hz55.30055.30039.10−30.003 次 150 Hz0.0000.0000—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% 的谱线其相位来自数值噪声统计时应剔除这一点在读取表格时要手工判断或者在代码里加阈值筛选。对合成电压信号执行上述代码输出示例谐波次数频率 (Hz)幅值 (V)有效值 (V)初相角 (deg)基波 150311.00219.9110.0210055.3039.10330.0315022.1015.6360.042000.000.00—这张表里的相位角已经做过 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_rmsS 会比所有谐波有功功率之和要大差值来自谐波无功分量。这就是非正弦电路中总功率因数始终低于基波 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 功率计算结果表与工程解读谐波次数 hU_h_rms (V)I_h_rms (A)P_h (W)φ_u−φ_i (deg)cosφ_h1219.9111.032396.418.00.951239.102.2754.1−75.00.259315.631.274.9−30.00.866直流2.500.250.63——基波功率 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; % 与 fs10000 不匹配 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 Hz50 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 Hz50 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 与嵌入式两套环境里都得到一致的工程结论。本文还有配套的精品资源点击获取
