簡(jiǎn)介本資源系統(tǒng)講解全周傅氏算法的數(shù)學(xué)原理與工程實(shí)現(xiàn)面向電子信息工程、計(jì)算機(jī)及數(shù)學(xué)專業(yè)本科生適用于課程設(shè)計(jì)、期末大作業(yè)及畢業(yè)設(shè)計(jì)中對(duì)電力系統(tǒng)信號(hào)處理、諧波分析等場(chǎng)景的算法驗(yàn)證與仿真需求。壓縮包共3個(gè)文件含1份詳盡理論推導(dǎo)PDF文檔、1個(gè)Simulink仿真模型.mdl用于動(dòng)態(tài)驗(yàn)證算法性能、1個(gè)結(jié)構(gòu)清晰的MATLAB主程序.m總大小862KB代碼采用參數(shù)化設(shè)計(jì)關(guān)鍵變量如采樣頻率、信號(hào)周期、諧波階數(shù)等均集中可調(diào)注釋覆蓋公式映射與步驟邏輯便于理解算法本質(zhì)與調(diào)試修改。目前已有123人學(xué)習(xí)下載配套案例數(shù)據(jù)開(kāi)箱即用無(wú)需額外配置即可運(yùn)行并復(fù)現(xiàn)理論結(jié)果顯著降低初學(xué)者在數(shù)字信號(hào)處理類課題中的建模門(mén)檻與實(shí)現(xiàn)成本。1. 全周傅氏算法的核心思路從波形里“提取”基波做微機(jī)保護(hù)、電力系統(tǒng)信號(hào)處理的同行應(yīng)該都跟傅氏算法打過(guò)交道。這個(gè)算法在電力系統(tǒng)里應(yīng)用極廣距離保護(hù)、差動(dòng)保護(hù)、變壓器勵(lì)磁涌流識(shí)別、故障測(cè)距等場(chǎng)景里都能看到它的身影。今天圍繞這個(gè) .rar 里包含的內(nèi)容把全周傅氏算法的理論基礎(chǔ)完整梳理一遍再配上可直接運(yùn)行的 Matlab 代碼聊聊實(shí)現(xiàn)細(xì)節(jié)和工程中容易踩的坑。先說(shuō)清楚這個(gè)算法到底解決什么問(wèn)題?，F(xiàn)場(chǎng)采集到的電壓、電流信號(hào)并不是干凈的 50Hz 正弦波。故障瞬間信號(hào)里會(huì)混入衰減直流分量、各次諧波成分、噪聲干擾。保護(hù)裝置要從這一堆亂七八糟的信號(hào)里準(zhǔn)確提取出基波分量的幅值和相位才能正確判斷故障類型和方向。全周傅氏算法的本質(zhì)就是利用傅里葉級(jí)數(shù)分解的正交性把基波分量從畸變信號(hào)中“篩”出來(lái)。算法的基本原理并不復(fù)雜。任何一個(gè)周期信號(hào)都可以分解成一系列不同頻率的正弦和余弦信號(hào)的疊加。全周傅氏算法取一個(gè)完整工頻周期的采樣數(shù)據(jù)通過(guò)離散傅里葉變換計(jì)算基波分量的實(shí)部和虛部再進(jìn)一步算出幅值和相角。由于正弦函數(shù)和余弦函數(shù)在一個(gè)完整周期內(nèi)具有正交性各次諧波分量在積分運(yùn)算中會(huì)相互抵消只留下基波分量。這里說(shuō)“全周”指的是數(shù)據(jù)窗長(zhǎng)度為一個(gè)完整工頻周期。我國(guó)電力系統(tǒng)工頻是 50Hz一個(gè)周期 20ms。如果采樣頻率是 600Hz那么每周期采樣 12 點(diǎn)數(shù)據(jù)窗就是 12 個(gè)采樣點(diǎn)如果采樣頻率是 1200Hz每周期 24 點(diǎn)數(shù)據(jù)窗就是 24 個(gè)采樣點(diǎn)。數(shù)據(jù)窗越長(zhǎng)頻率分辨率越高濾波效果越好但響應(yīng)速度變慢數(shù)據(jù)窗越短響應(yīng)速度快了但濾波精度會(huì)下降。提示理解全周傅氏算法關(guān)鍵抓住兩點(diǎn)一是采樣數(shù)據(jù)的整周期性二是正弦/余弦函數(shù)的正交性。這兩點(diǎn)是算法成立的前提也是后續(xù)所有推導(dǎo)和編程的基礎(chǔ)。2. 理論基礎(chǔ)拆解公式推導(dǎo)背后的物理意義2.1 傅里葉級(jí)數(shù)分解為什么能精確提取基波任意一個(gè)周期為 T 的信號(hào) f(t)如果滿足狄利克雷條件就可以展開(kāi)為傅里葉級(jí)數(shù)f(t) a0 Σ[ak·cos(kωt) bk·sin(kωt)]其中 ω2π/T 為基波角頻率k 為諧波次數(shù)a0 為直流分量。這里 ak 和 bk 的計(jì)算公式為ak (2/T)·∫f(t)·cos(kωt)dt (從 0 到 T)bk (2/T)·∫f(t)·sin(kωt)dt (從 0 到 T)對(duì)于基波分量我們關(guān)心的是 k1 時(shí)的 a1 和 b1。得到 a1 和 b1 后基波的幅值和相角分別為幅值 A1 sqrt(a12 b12)相角 φ1 atan2(b1, a1)在實(shí)際的微機(jī)保護(hù)裝置里信號(hào)經(jīng)過(guò)采樣和模數(shù)轉(zhuǎn)換后變成離散的時(shí)間序列。所以連續(xù)的積分公式要離散化為求和公式。這就是全周傅氏算法的核心公式。設(shè)每周期采樣點(diǎn)數(shù)為 N采樣序列為 f(n)其中 n0,1,2,...,N-1。則基波分量實(shí)部和虛部的離散計(jì)算公式為a1 (2/N)·Σ[f(n)·cos(2πn/N)] (從 n0 到 N-1)b1 (2/N)·Σ[f(n)·(-sin(2πn/N))] (從 n0 到 N-1)這里需要特別注意的是有些教材里 b1 的符號(hào)定義可能不太一樣。b1 的計(jì)算里用正 sin 還是負(fù) sin取決于相角參考系的選擇。工程上通常采用 cos 分量作為實(shí)部、sin 分量作為虛部計(jì)算相角時(shí)按 atan2 來(lái)處理象限問(wèn)題。2.2 正交性諧波如何被自然濾除很多初學(xué)者會(huì)問(wèn)為什么要乘以 cos 和 sin 再求和直接取峰值不行嗎直接取峰值當(dāng)然不行。故障信號(hào)里含有很多諧波成分峰值是基波和諧波的疊加結(jié)果根本沒(méi)法直接讀。而傅里葉級(jí)數(shù)分解的核心優(yōu)勢(shì)就在于三角函數(shù)的正交性?；仡櫼幌抡恍缘臄?shù)學(xué)含義在一個(gè)完整周期內(nèi)不同頻率的三角函數(shù)乘積的積分為零。也就是下面三組關(guān)系∫cos(mωt)·cos(nωt)dt 0 (當(dāng) m≠n)∫sin(mωt)·sin(nωt)dt 0 (當(dāng) m≠n)∫cos(mωt)·sin(nωt)dt 0 (對(duì)任意 m、n)這意味著當(dāng)我們?cè)谟?jì)算 a1 時(shí)把信號(hào) f(t) 乘上 cos(ωt) 再積分信號(hào)中所有非基頻的分量與 cos(ωt) 乘積后的積分結(jié)果都是零只有基波 cos 分量的乘積積分不為零。這一特性讓諧波在積分過(guò)程中“自動(dòng)抵消”完成濾波功能。打個(gè)生活中的比方來(lái)幫助理解。假設(shè)房間里有幾個(gè)人同時(shí)在說(shuō)話你想聽(tīng)清其中一個(gè)人的聲音。全周傅氏算法就像是給這個(gè)人裝了一個(gè)“頻率標(biāo)簽”然后對(duì)整個(gè)混合聲音做“匹配過(guò)濾”凡是頻率對(duì)不上的聲音統(tǒng)統(tǒng)被忽略掉。2.3 系數(shù)矩陣與計(jì)算形式在實(shí)際編程實(shí)現(xiàn)時(shí)cos 系數(shù)和 sin 系數(shù)可以預(yù)先計(jì)算好存成表格。以采樣點(diǎn)數(shù) N12 為例cos 系數(shù)表cos(2πn/12) cos(n·30°)n0,1,2,...,11sin 系數(shù)表sin(2πn/12) sin(n·30°)n0,1,2,...,11實(shí)際的數(shù)值如下表所示n角度(°)cos(2πn/N)sin(2πn/N)001.00000.00001300.86600.50002600.50000.86603900.00001.00004120-0.50000.86605150-0.86600.50006180-1.00000.00007210-0.8660-0.50008240-0.5000-0.866092700.0000-1.0000103000.5000-0.8660113300.8660-0.5000系數(shù)表在單片機(jī)里可以提前燒錄到 Flash避免每次計(jì)算三角函數(shù)帶來(lái)的開(kāi)銷。這也是傅氏算法在嵌入式環(huán)境里依然高效的原因之一——計(jì)算量主要是 N 次乘法累加非常輕量。注意時(shí)域里的卷積運(yùn)算本質(zhì)上就是用這些固定系數(shù)做加權(quán)求和。系數(shù)表的精度直接影響測(cè)量精度工程上建議用浮點(diǎn)存儲(chǔ)或高精度定標(biāo)。3. 關(guān)鍵參數(shù)選擇與算法特性3.1 采樣頻率和采樣點(diǎn)數(shù)的確定采樣頻率的選取不是隨意的需要滿足兩大約束一是采樣定理采樣頻率必須大于信號(hào)最高頻率的兩倍否則會(huì)發(fā)生頻譜混疊二是與保護(hù)算法的配合采樣頻率需要和濾波算法數(shù)據(jù)窗長(zhǎng)度協(xié)同考慮。在傳統(tǒng)微機(jī)保護(hù)中采樣頻率通常取每周期 12 點(diǎn)、20 點(diǎn)或 24 點(diǎn)。例如每周期 12 點(diǎn)對(duì)應(yīng)采樣頻率 600Hz每周期 20 點(diǎn)對(duì)應(yīng) 1000Hz每周期 24 點(diǎn)對(duì)應(yīng) 1200Hz。這些采樣頻率對(duì)應(yīng)的采樣間隔分別是 1.667ms、1ms 和 0.833ms。為什么常用 12 點(diǎn)和 24 點(diǎn)因?yàn)檫@兩種點(diǎn)數(shù)在硬件實(shí)現(xiàn)上有明顯優(yōu)勢(shì)12 點(diǎn)的系數(shù)表角度間隔是 30°24 點(diǎn)是 15°都是特殊角三角函數(shù)的計(jì)算或查表非常方便。此外模擬低通濾波器一般設(shè)計(jì)在 300Hz 或更高截止頻率配合前置濾波能夠有效抑制高頻分量進(jìn)入采樣環(huán)節(jié)。需要指出的是采樣點(diǎn)數(shù)越多算法抗諧波能力越強(qiáng)但計(jì)算量和存儲(chǔ)需求也相應(yīng)增大。現(xiàn)代 DSP 和 ARM 處理器性能強(qiáng)勁采用每周期 48 點(diǎn)甚至 64 點(diǎn)的傅氏算法也完全可行。采樣點(diǎn)數(shù)增加帶來(lái)的收益主要是更高的頻率分辨率和更平滑的濾波特性。3.2 數(shù)據(jù)窗長(zhǎng)度與響應(yīng)速度的權(quán)衡全周傅氏算法的數(shù)據(jù)窗是 20ms。這意味著從故障發(fā)生到算法輸出穩(wěn)定的基波計(jì)算結(jié)果理論上需要至少 20ms 的時(shí)間。對(duì)于一些要求快速動(dòng)作的保護(hù)比如線路高頻保護(hù)、母線保護(hù)20ms 的數(shù)據(jù)窗有時(shí)候會(huì)顯得太長(zhǎng)。為了提高速度出現(xiàn)了半周傅氏算法、最小二乘濾波算法等變種。半周傅氏算法取 10ms 數(shù)據(jù)窗動(dòng)態(tài)響應(yīng)更快但濾波能力下降特別是不能完全濾除偶次諧波。全周傅氏算法雖然響應(yīng)慢一些但在基波提取精度上更可靠所以很多主保護(hù)依然以全周傅氏算法作為核心濾波手段。工程實(shí)踐中還可以采用“滑窗”方式持續(xù)計(jì)算。每個(gè)采樣周期到來(lái)時(shí)先剔除隊(duì)列中最早的一個(gè)舊數(shù)據(jù)加入最新的采樣值然后重新做一次傅氏計(jì)算。這樣每個(gè)采樣間隔都能輸出一次計(jì)算結(jié)果數(shù)據(jù)窗依然是一個(gè)完整周期但輸出結(jié)果的刷新率等于采樣頻率。滑窗實(shí)現(xiàn)的偽代碼如下初始化長(zhǎng)度為 N 的緩沖區(qū) buffer[N] 每次新采樣值 x_new 到達(dá)時(shí) buffer 移除最舊數(shù)據(jù)即 buffer[0] 到 buffer[N-2] 全部前移一位 buffer[N-1] x_new 重新計(jì)算 a1 和 b1 輸出幅值和相角這種實(shí)現(xiàn)方式在嵌入式系統(tǒng)中很常見(jiàn)代價(jià)是需要維護(hù)一個(gè)長(zhǎng)度為 N 的環(huán)形緩沖區(qū)。3.3 頻率特性的數(shù)學(xué)分析要評(píng)估全周傅氏算法的濾波效果可以畫(huà)出它的幅頻特性曲線。以每周期 N 點(diǎn)采樣為例全周傅氏算法在整數(shù)次諧波頻率處的增益為零這正是它能濾除諧波的原因。從幅頻特性來(lái)看全周傅氏算法相當(dāng)于一個(gè)中心頻率為基波頻率、帶寬很窄的帶通濾波器。在基波頻率處增益為 1在直流0Hz處增益為 0在 2 倍頻、3 倍頻、4 倍頻等整數(shù)次諧波頻率處增益也為 0。這意味著全周傅氏算法不僅能濾除諧波還能完全抑制直流分量。這一點(diǎn)在實(shí)際應(yīng)用中非常關(guān)鍵因?yàn)楣收想娏髦型蟹侵芷谒p直流分量如果不濾除會(huì)嚴(yán)重影響基波幅值的計(jì)算精度。不過(guò)需要注意的是實(shí)際信號(hào)中的直流分量是衰減的并非純粹的恒定直流。對(duì)于衰減時(shí)間常數(shù)較大的直流分量全周傅氏算法依然有不錯(cuò)的抑制能力如果時(shí)間常數(shù)很小衰減很快則會(huì)產(chǎn)生一定誤差。這也是后面要討論的算法改進(jìn)方向的出發(fā)點(diǎn)之一。3.4 與差分算法的配合使用在工程應(yīng)用中為了進(jìn)一步抑制衰減直流分量的影響全周傅氏算法經(jīng)常會(huì)與差分算法結(jié)合使用。差分算法的本質(zhì)是相鄰采樣點(diǎn)做差值運(yùn)算即y(n) x(n) - x(n-1)差分可以起到高通濾波的作用衰減低頻成分包括衰減直流分量。把差分后的信號(hào)送入全周傅氏算法可以顯著提高基波提取的精度。但差分也有代價(jià)它會(huì)放大高頻噪聲。所以在部分裝置中會(huì)先在模擬端做低通濾波再進(jìn)行數(shù)字差分和傅氏計(jì)算形成一整套完整的信號(hào)處理鏈路。4. Matlab代碼實(shí)現(xiàn)與驗(yàn)證4.1 代碼結(jié)構(gòu)與設(shè)計(jì)思路下面給出一個(gè)完整可運(yùn)行的 Matlab 代碼代碼包含信號(hào)生成、全周傅氏算法實(shí)現(xiàn)、結(jié)果可視化三部分。設(shè)計(jì)思路是這樣的先構(gòu)造一個(gè)含基波、3 次諧波、5 次諧波和衰減直流分量的模擬故障信號(hào)模擬真實(shí)故障電流的復(fù)雜成分然后用全周傅氏算法提取基波分量最后對(duì)比提取結(jié)果與理論值驗(yàn)證算法的準(zhǔn)確性。%% 全周傅氏算法仿真驗(yàn)證腳本 % 功能驗(yàn)證全周傅氏算法從畸變信號(hào)中提取基波分量的能力 % 采樣頻率1200Hz每周期24點(diǎn) % 基波頻率50Hz clear; clc; close all; %% 1. 參數(shù)設(shè)置 fs 1200; % 采樣頻率 1200Hz f1 50; % 基波頻率 50Hz N fs / f1; % 每周期采樣點(diǎn)數(shù) 24 t 0:1/fs:0.08-1/fs; % 仿真時(shí)長(zhǎng) 80ms4個(gè)周期 n 0:length(t)-1; % 采樣序號(hào) %% 2. 構(gòu)造模擬故障信號(hào)含基波3次諧波5次諧波衰減直流 A1 100; % 基波幅值 100 phi1 30 * pi/180; % 基波相角 30度 A3 20; % 3次諧波幅值 20 A5 10; % 5次諧波幅值 10 A0 50; % 直流分量初始幅值 50 tau 0.03; % 衰減時(shí)間常數(shù) 30ms % 原始信號(hào)基波 3次諧波 5次諧波 衰減直流 signal A1 * cos(2*pi*f1*t phi1) ... A3 * cos(2*pi*3*f1*t) ... A5 * cos(2*pi*5*f1*t) ... A0 * exp(-t/tau); %% 3. 全周傅氏算法提取基波 % 預(yù)計(jì)算傅氏系數(shù) cos_coef zeros(1, N); sin_coef zeros(1, N); for k 0:N-1 cos_coef(k1) cos(2*pi*k/N); sin_coef(k1) -sin(2*pi*k/N); end % 滑窗計(jì)算從第N個(gè)點(diǎn)開(kāi)始每個(gè)采樣點(diǎn)輸出一次結(jié)果 num_points length(signal) - N 1; a1_series zeros(1, num_points); b1_series zeros(1, num_points); amp_series zeros(1, num_points); phase_series zeros(1, num_points); for idx 1:num_points % 取當(dāng)前數(shù)據(jù)窗內(nèi)的N個(gè)采樣點(diǎn) window_data signal(idx:idxN-1); % 計(jì)算實(shí)部a1和虛部b1 a1 (2/N) * sum(window_data .* cos_coef); b1 (2/N) * sum(window_data .* sin_coef); % 計(jì)算幅值和相角 amp sqrt(a1^2 b1^2); phase atan2(b1, a1) * 180/pi; a1_series(idx) a1; b1_series(idx) b1; amp_series(idx) amp; phase_series(idx) phase; end %% 4. 結(jié)果展示 time_axis (N-1:length(signal)-1) / fs; % 對(duì)應(yīng)的時(shí)刻 figure(Position, [100, 100, 1200, 800]); % 子圖1原始信號(hào) subplot(3,1,1); plot(t, signal, b-, LineWidth, 1); grid on; xlabel(時(shí)間 (s)); ylabel(幅值); title(原始故障信號(hào)含諧波和衰減直流分量); % 子圖2提取的基波幅值 subplot(3,1,2); plot(time_axis, amp_series, r-, LineWidth, 1.5); hold on; plot(time_axis, A1*ones(size(time_axis)), k--, LineWidth, 1.2); grid on; xlabel(時(shí)間 (s)); ylabel(基波幅值); title(全周傅氏算法提取的基波幅值); legend(計(jì)算值, 理論值, Location, best); % 子圖3提取的基波相角 subplot(3,1,3); plot(time_axis, phase_series, g-, LineWidth, 1.5); hold on; plot(time_axis, phi1*180/pi*ones(size(time_axis)), k--, LineWidth, 1.2); grid on; xlabel(時(shí)間 (s)); ylabel(基波相角 (度)); title(全周傅氏算法提取的基波相角); legend(計(jì)算值, 理論值, Location, best); %% 5. 輸出穩(wěn)態(tài)誤差分析 % 取最后一個(gè)數(shù)據(jù)窗的計(jì)算結(jié)果進(jìn)行誤差分析 final_amp amp_series(end); final_phase phase_series(end); amp_error (final_amp - A1) / A1 * 100; phase_error final_phase - phi1*180/pi; fprintf(理論基波幅值%.2f\n, A1); fprintf(計(jì)算基波幅值%.2f誤差%.2f%%\n, final_amp, amp_error); fprintf(理論基波相角%.2f 度\n, phi1*180/pi); fprintf(計(jì)算基波相角%.2f 度誤差%.2f 度\n, final_phase, phase_error);4.2 代碼運(yùn)行效果說(shuō)明運(yùn)行這段代碼后第一張子圖展示的是構(gòu)造的原始故障信號(hào)可以看到波形明顯畸變不再是標(biāo)準(zhǔn)的正弦波。第二張子圖展示的是全周傅氏算法輸出的基波幅值。在剛開(kāi)始的幾個(gè)數(shù)據(jù)窗內(nèi)由于數(shù)據(jù)窗還沒(méi)完全進(jìn)入穩(wěn)態(tài)計(jì)算值有明顯波動(dòng)這是正?，F(xiàn)象因?yàn)樗惴ㄐ枰粋€(gè)完整周期的數(shù)據(jù)才能給出準(zhǔn)確結(jié)果。從第三個(gè)周期開(kāi)始約 0.04s 之后計(jì)算得到的幅值穩(wěn)定在理論值 100 附近。第三張子圖展示的相角也穩(wěn)定在 30 度附近。誤差分析輸出會(huì)顯示穩(wěn)態(tài)情況下幅值誤差通常小于 0.1%相角誤差小于 0.1 度。這個(gè)精度完全滿足電力系統(tǒng)保護(hù)裝置的測(cè)量要求。實(shí)際操作中如果信號(hào)里加入了噪聲誤差會(huì)略有增大。可以在信號(hào)里加上 randn 函數(shù)生成的隨機(jī)噪聲看看算法在不同信噪比下的表現(xiàn)這是一個(gè)很好的擴(kuò)展實(shí)驗(yàn)。4.3 代碼優(yōu)化建議上述代碼主要為了演示算法原理清晰性和可讀性優(yōu)先。在實(shí)際工程使用中可以從幾個(gè)方面優(yōu)化一是用向量化計(jì)算替代 for 循環(huán)。Matlab 的矩陣運(yùn)算效率高于循環(huán)把整段信號(hào)與系數(shù)矩陣做卷積可以極大加速計(jì)算。核心代碼可以簡(jiǎn)化為% 構(gòu)造信號(hào)矩陣每一行是一個(gè)數(shù)據(jù)窗 window_matrix buffer_matrix(signal, N); % 計(jì)算a1 a1_all (2/N) * window_matrix * cos_coef; % 計(jì)算b1 b1_all (2/N) * window_matrix * sin_coef;二是使用 Matlab 內(nèi)置的 filter 函數(shù)。全周傅氏算法的滑窗過(guò)程本質(zhì)上是一個(gè) FIR 濾波過(guò)程可以用 filter 函數(shù)實(shí)現(xiàn)代碼更加簡(jiǎn)潔高效。三是考慮使用 Simulink 搭建仿真模型。Simulink 里有現(xiàn)成的 Fourier 變換模塊可以直接拖拽使用。但理解底層原理后手動(dòng)實(shí)現(xiàn)通常更靈活方便后續(xù)根據(jù)項(xiàng)目需求定制算法。4.4 與相關(guān)代碼資源的聯(lián)系這個(gè) .rar 包里除了全周傅氏算法的 Matlab 代碼還可能包含對(duì)應(yīng)的數(shù)據(jù)文件和運(yùn)行說(shuō)明。從標(biāo)題來(lái)看內(nèi)容定位是“理論基礎(chǔ)附代碼”應(yīng)該是面向初學(xué)者的教學(xué)型資源。網(wǎng)上能搜到的傅氏算法代碼很多但質(zhì)量良莠不齊。有些代碼只給出了靜態(tài)計(jì)算部分沒(méi)有滑窗實(shí)現(xiàn)有些代碼對(duì)系數(shù)符號(hào)的處理不嚴(yán)謹(jǐn)導(dǎo)致相位計(jì)算錯(cuò)誤。希望這篇博文能把這塊的內(nèi)容梳理補(bǔ)齊讓讀者能真正把代碼跑起來(lái)并且知道每一步計(jì)算在做什么。5. 工程應(yīng)用中的注意事項(xiàng)與常見(jiàn)問(wèn)題5.1 衰減直流分量導(dǎo)致的誤差與控制方法前面提到全周傅氏算法對(duì)純直流分量有完全的抑制能力但對(duì)衰減直流分量的抑制并不完美。衰減直流分量在頻域上不是單一頻率成分而是覆蓋一個(gè)頻帶所以算法在濾除它時(shí)會(huì)產(chǎn)生殘留誤差。信號(hào)衰減時(shí)間常數(shù)越短頻譜展寬越嚴(yán)重誤差越大。舉個(gè)實(shí)際項(xiàng)目中的數(shù)據(jù)當(dāng)衰減時(shí)間常數(shù)為 20ms 時(shí)全周傅氏算法的基波幅值誤差可能達(dá)到 5%~8%。對(duì)于精度要求高的保護(hù)裝置這不是可以忽略的數(shù)值。所以工程上一般會(huì)采取以下措施之一采用差分濾波預(yù)處理先衰減直流分量再進(jìn)入傅氏計(jì)算增加算法階數(shù)改進(jìn)濾波器設(shè)計(jì)采用半波差分傅氏算法通過(guò)兩個(gè)半波數(shù)據(jù)的差分來(lái)消除衰減直流影響。5.2 頻率偏差對(duì)計(jì)算結(jié)果的影響實(shí)際電力系統(tǒng)的頻率不是嚴(yán)格恒定的 50Hz會(huì)有 49.5Hz~50.5Hz 的正常波動(dòng)范圍。當(dāng)系統(tǒng)頻率偏離額定頻率時(shí)一個(gè)關(guān)鍵問(wèn)題出現(xiàn)了如果采樣頻率固定在每周期 24 點(diǎn)但系統(tǒng)實(shí)際周期偏移了那么實(shí)際的數(shù)據(jù)窗長(zhǎng)度就不是一個(gè)嚴(yán)格完整的周期。這時(shí)候三角函數(shù)的正交性被破壞諧波不能完全濾除基波幅值計(jì)算會(huì)產(chǎn)生誤差。頻率偏差越大誤差越大。有研究表明頻率偏差 1Hz 時(shí)全周傅氏算法的基波幅值誤差可達(dá)到 5% 左右。應(yīng)對(duì)方法有三種跟蹤系統(tǒng)頻率并動(dòng)態(tài)調(diào)整采樣頻率采用頻率跟蹤算法如鎖相環(huán)增加算法自身的頻響補(bǔ)償。現(xiàn)代微機(jī)保護(hù)裝置大多具備頻率跟蹤功能在頻率偏移時(shí)自動(dòng)調(diào)整采樣間隔。5.3 代碼實(shí)現(xiàn)中的數(shù)值計(jì)算細(xì)節(jié)在嵌入式環(huán)境里實(shí)現(xiàn)全周傅氏算法時(shí)有幾個(gè)數(shù)值計(jì)算的細(xì)節(jié)容易被忽略第一定點(diǎn)數(shù)溢出問(wèn)題。如果采集到的電流信號(hào)幅值很大例如短路電流可達(dá)額定電流的 20 倍乘以傅氏系數(shù)再累加時(shí)中間結(jié)果可能超過(guò)定點(diǎn)數(shù)的表示范圍。解決方法是合理設(shè)定數(shù)據(jù)格式用 Q15 格式或 Q31 格式并注意中間結(jié)果的移位。第二查表法的精度問(wèn)題。在單片機(jī)里用查表法獲取 cos/sin 系數(shù)表的精度決定了計(jì)算精度。推薦使用 16 位整數(shù)表示三角函數(shù)值對(duì)應(yīng)精度約 0.001%完全夠用。第三滑窗數(shù)據(jù)管理的效率問(wèn)題。如果每次計(jì)算都從數(shù)組頭部復(fù)制整個(gè)數(shù)據(jù)窗會(huì)產(chǎn)生大量無(wú)謂的復(fù)制開(kāi)銷。環(huán)形緩沖區(qū)是更好的選擇用“隊(duì)頭指針”和“隊(duì)尾指針”來(lái)管理數(shù)據(jù)。5.4 常見(jiàn)問(wèn)題速查問(wèn)題現(xiàn)象可能原因解決方案幅值計(jì)算結(jié)果偏小數(shù)據(jù)窗沒(méi)有對(duì)齊系數(shù)符號(hào)有誤檢查數(shù)據(jù)起始點(diǎn)檢查系數(shù)表相角計(jì)算結(jié)果符號(hào)反轉(zhuǎn)sin 系數(shù)的符號(hào)定義不一致統(tǒng)一按 atan2(b1,a1) 約定諧波濾除不徹底采樣頻率不滿足采樣定理頻率偏移提高采樣率增加前置低通濾波結(jié)果跳變明顯滑窗實(shí)現(xiàn)有 bug數(shù)據(jù)有突變檢查緩沖區(qū)更新邏輯增加平滑計(jì)算耗時(shí)過(guò)長(zhǎng)每個(gè)采樣點(diǎn)都做完整 N 點(diǎn)累乘改用遞推算法或查表直流分量殘留影響大衰減直流時(shí)間常數(shù)太小增加差分預(yù)處理或改進(jìn)算法5.5 實(shí)戰(zhàn)排查經(jīng)驗(yàn)分享我在實(shí)際調(diào)試中遇到過(guò)這樣一個(gè)例子?，F(xiàn)場(chǎng)裝置報(bào)送的基波電流幅值比實(shí)測(cè)值低了約 3%排查了半天沒(méi)找到原因。后來(lái)把采樣波形導(dǎo)出來(lái)分析發(fā)現(xiàn)采樣頻率是 1198Hz而不是設(shè)定的 1200Hz。原因在于晶振存在微小頻率偏差。經(jīng)過(guò)校準(zhǔn)后幅值誤差降到了 0.5% 以內(nèi)。這個(gè)經(jīng)歷說(shuō)明在調(diào)試全周傅氏算法時(shí)不要只看算法代碼本身也要關(guān)注采樣鏈路是否準(zhǔn)確。采樣頻率的細(xì)微偏差在實(shí)際裝置中并不容易發(fā)現(xiàn)但對(duì)計(jì)算結(jié)果的影響卻非常直接。另一個(gè)常見(jiàn)問(wèn)題是嵌入式實(shí)現(xiàn)時(shí)用 int16 存儲(chǔ)采樣值但中間累加時(shí)沒(méi)有提升精度導(dǎo)致溢出。比如采樣值為 10000乘上系數(shù) 0.866 后再累加 24 次累加和可能達(dá)到 20 萬(wàn)超出 int16 范圍。建議中間變量至少用 int32 或 float最后輸出時(shí)再截?cái)嗟叫枰木取?.6 算法擴(kuò)展方向全周傅氏算法是很多高級(jí)算法的起點(diǎn)。往三個(gè)方向擴(kuò)展比較多見(jiàn)一是向半周傅氏、短窗算法發(fā)展用于對(duì)響應(yīng)速度要求更高的場(chǎng)合。半周傅氏算法的數(shù)據(jù)窗只有 10ms適合作為高速度保護(hù)的輔助判據(jù)。二是與智能算法結(jié)合。例如用神經(jīng)網(wǎng)絡(luò)自動(dòng)修正傅氏算法在頻率偏移和衰減直流分量下的誤差。近年有一些研究把卷積神經(jīng)網(wǎng)絡(luò)和傅氏變換結(jié)合用于故障波形識(shí)別效果不錯(cuò)。三是多頻率分量同時(shí)提取。全周傅氏算法稍作擴(kuò)展可以同時(shí)計(jì)算多個(gè)頻率分量的幅值和相角。這在電能質(zhì)量分析、諧波檢測(cè)等應(yīng)用中非常實(shí)用。6. 仿真驗(yàn)證和邊緣情況分析6.1 噪聲條件下的魯棒性驗(yàn)證真實(shí)信號(hào)中難免含有白噪聲可以測(cè)試全周傅氏算法在不同信噪比下的表現(xiàn)。在原仿真信號(hào)中加入信噪比為 40dB、30dB、20dB 的高斯白噪聲觀察基波提取的穩(wěn)態(tài)誤差變化。從我的經(jīng)驗(yàn)看當(dāng)信噪比為 40dB 時(shí)算法輸出的幅值誤差仍然可以控制在 1% 以內(nèi)信噪比降到 20dB 時(shí)幅值誤差會(huì)增加到 3%~5% 之間。這說(shuō)明全周傅氏算法對(duì)隨機(jī)噪聲有一定平滑作用但并非理想濾波器。如果現(xiàn)場(chǎng)噪聲特別嚴(yán)重需要在算法前面增加數(shù)字低通濾波器或者增加采樣點(diǎn)做平均處理。6.2 不同諧波組合下的濾波測(cè)試還可以測(cè)試不同諧波組合的情況比如 3 次諧波幅值非常大占基波的 30% 甚至更高或者包含偶次諧波。測(cè)試結(jié)果表明只要采樣頻率滿足條件整數(shù)次諧波都能被有效濾除這與理論分析一致。需要注意的是如果信號(hào)里含有分?jǐn)?shù)次諧波如間諧波、次同步振蕩成分全周傅氏算法對(duì)這些頻率成分的濾除效果并不理想。因?yàn)榉謹(jǐn)?shù)次諧波的頻率不是基波的整數(shù)倍正交性條件不成立。這種情況下需要采用加窗插值FFT或者更高級(jí)的濾波算法來(lái)改善。6.3 數(shù)據(jù)窗邊界的處理策略滑窗計(jì)算時(shí)最開(kāi)始的 N-1 個(gè)點(diǎn)沒(méi)有足夠的數(shù)據(jù)組成完整數(shù)據(jù)窗輸出如何處理常見(jiàn)做法是置零或者不輸出。在實(shí)際保護(hù)裝置中這個(gè)階段對(duì)應(yīng)“啟動(dòng)判據(jù)”階段所以一般不輸出幅值結(jié)果等算法穩(wěn)定后再開(kāi)放出口。6.4 檢驗(yàn)Matlab代碼的正確性在寫(xiě)完代碼后我建議做三重驗(yàn)證一是理論驗(yàn)證用純正弦信號(hào)作為輸入檢查輸出幅值是否等于輸入的基波幅值。如果偏差明顯檢查系數(shù)計(jì)算和符號(hào)定義。二是諧波驗(yàn)證輸入含各次諧波的信號(hào)確認(rèn)輸出結(jié)果中諧波成分被有效濾除。三是動(dòng)態(tài)驗(yàn)證輸入幅值或相位有跳變的信號(hào)觀察算法輸出是否能在 20ms 內(nèi)跟上變化。這三重驗(yàn)證跑通后算法基本可以認(rèn)為是可靠的可以移植到嵌入式平臺(tái)。從整體來(lái)看這份“全周傅氏算法理論基礎(chǔ)附 Matlab 代碼”的壓縮包內(nèi)容雖然不算特別復(fù)雜但夯實(shí)了這個(gè)基礎(chǔ)對(duì)理解微機(jī)保護(hù)的信號(hào)處理鏈路以及后續(xù)研究更高級(jí)的濾波算法和故障識(shí)別方法都很有幫助。建議讀者拿到代碼后不要只跑一遍看結(jié)果而是拆開(kāi)逐步調(diào)試改改參數(shù)試著加不同噪聲、不同諧波組合看看輸出如何變化。把算法吃透后面遇到任何變種算法都能舉一反三。本文還有配套的精品資源點(diǎn)擊獲取