电力系统谐波检测加汉宁窗插值算法Matlab实现

快速傅里叶变换是应用最广泛的一种谐波检测方法，但直接利用快速傅里叶变换进行谐波检测存在较大的误差，影响谐波分析结果的准确性。通过加汉宁窗及插值修正算法可以改善计算谐波频率、相位和幅值的准确度。简述了电力系统谐波检测非同步采样加汉宁窗插值算法的原理，并采用巴特沃斯低通滤波器滤除高频噪声。MATLAB仿真结果表明，加汉宁窗插值算法具有检测精度好，实现简单的优点。

% Hanning加窗插值算法
% cgz  2008.11.19

clear;
% 原始数据：直流:0V; 基波：49.8Hz100V10deg; HR2:0.5V40deg;....
hr0=0;f1=49.8;
hr（1）=1;deg（1）=pi/12;
hr（2）=2;deg（2）=pi/3;
hr（3）=3;deg（3）=pi/7;
hr（4）=9;deg（4）=pi/6;
hr（5）=4;deg（5）=pi/4;
hr（6）=6;deg（6）=2*pi/3;
hr（7）=8;deg（7）=pi;
M=7;f=[1:M]*f1; %设定频率

% 采样：N=1024 fs=5000Hz.
fs=5000;
N=1024;    
T=1/fs;
n=[0:N-1];t=n*T;
x=zeros（size（t））;
for k=1:M
    x=x+hr（k）*cos（2*pi*f（k）*t+deg（k）*pi/180）;
end

%加入高频噪声
y1=x+5*cos（2*pi*1000*t）;
%滤波
Y1=fft（y1）;Y1=fftshift（Y1）/1024;
N1=length（y1）;
w1=2*[0:N1-1]./N1-1;f_1=w1/（2*T）;
figure（1）
plot（f_1abs（Y1））;grid;

%滤波器
alphamax=1;alphamin=10;
wp=0.2;ws=0.22;
[N2wh]=buttord（wpwsalphamaxalphamin）;
[ba]=butter（N2wh）;
y2=filter（bay1）;
Y2=fft（y2）;Y2=fftshift（Y2）/512;
figure（2）
plot（f_1abs（Y2））;grid;

% 分析：
x=y2;
w=0.5-0.5*cos（2*pi*n/N）;    % hanning window。
Xk=fft（x.*w）;
amp=abs（Xk（1:N/2））/N*2;         %幅频
pha=angle（Xk（1:N/2））/pi*180;    %相频
for k=1:N/2
    if（amp（k）<0.01） pha（k）=0;  %当谐波<10mV时，其相位=0
    end
    if（pha（k）<0） pha（k）=pha（k）+360;%调整到0-360度
    end
end
fmin=fs/N;
xaxis=fmin*n（1:N/2）;    %横坐标为Hz
kx=round（[1:M]*50/fmin）;%各次谐波对应的下标（从0开始）
for m=1:M
    km（m）=searchpeaks（ampkx（m）+1）;  %km为谱峰（从1开始）
    if（amp（km（m）+1）        km（m）=km（m）-1;
    end
    beta（m）=amp（km（m）+1）./amp（km（m））;
    delta（m）=（2*beta（m）-1）./（1+beta（m））;
end
fx=（km-1+delta）*fmin;     %估计频率
hrx=amp（km）*2.*pi.*delta.*（1-delta.*delta）./sin（pi*delta）;  %估计幅度
degx=pha（km）-delta.*180/N*（N-1）; %估计相位
degx=mod（degx360）; %调整到0-360度
efx=（fx-f）./f*100;          %频率误差
ehr=（hrx-hr）./hr*100;       %幅度误差
edeg=（degx-deg）;   %相位误差  

% 结果输出：
figure（3）
subplot（221）;%画出采样序列
plot（tx）;
hold on;
plot（tx.*w‘r‘）;   %加窗波形
hold off;
xlabel（‘x（k）‘）;
title（‘采样序列‘）;
subplot（222）;%画出FFT分析结果
stem（xaxisamp‘.r‘）;
xlabel（‘normlized magnitude（*N/2）‘）;
title（‘FFT分析结果（幅频）‘）;
subplot（224）;
stem（xaxispha‘.r‘）;
xlabel（‘phase（deg）‘）;
title（‘FFT分析结果（相频）‘）;
subplot（223）;
plot（ehr）;
title（‘幅度误差（%）‘）;

%文本输出
disp ‘    分析结果的数值见result.txt文件.‘;
fid=fopen（‘result.txt‘‘w‘）;
fprintf（fid‘原始数据：f1=%6.1fHz N=%.f  fs=%.f \r\n\r\n‘f1Nfs）;
fprintf（fid‘谐波次数      1      2      3      4      5      6     7\r\n‘）;
fprintf（fid‘设定频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n‘f）;
fprintf（fid‘估计频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n‘fx）;
fprintf（fid‘误差（%%）  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n‘efx）;

fprintf（fid‘设定幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n‘hr）;
fprintf（fid‘估计幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n‘hrx）;
fprintf（fid‘误差（%%）  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n‘ehr）;

fprintf（fid‘设定相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n‘deg）;
fprintf（fid‘估计相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n‘degx）;
fprintf（fid‘误差（度）  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n\r\n‘edeg）;

%其他数据
fprintf（fid‘谱峰位置理论值：\r\n %6.4f %6.4f %6.4f

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件         933  2017-12-20 19:45  result.txt
     文件         422  2008-11-19 19:35  searchpeaks.m
     文件        3585  2017-12-20 19:45  hanning_win.m