基于MATLAB的电力谐波分析，在有频谱泄露、栅栏效应的前提下计算基波和谐波的频率、幅度、相位

基于MATLAB FFT的电力谐波分析，输入时域信号采样值，在有频谱泄露、栅栏效应的前提下计算基波和谐波的频率、幅度、相位。根据论文《基于加汉宁窗的FFT高精度谐波检测改进算法》实现。函数[f_HA,A_HA,phi_HA]=harmonic_analysis(x,fs,N_HA)实现频谱分析的功能，主函数test_FFT产生仿真采样信号并调用频谱分析函数，计算电力质量指标。代码中有详细的注释。

function test_FFT

%信号采样的参数
N_HA=9;%待分析的谐波个数
fs=1024;%采样频率，大于谐波最高频率乘以2，至少N_HA*50*2或者N_HA*60*2
Nsample=2048;%102420484096
Nbit=12;%ADC采样的位数取值81012
kADC=280*sqrt（2）/2^（Nbit-1）;%ADC的比例系数，允许最大电压为正负280*sqrt（2）

%电力质量数据缓存
Sizebuffer=10000;
U1buffer=zeros（1Sizebuffer）;%存储U1的缓存
f1buffer=zeros（1Sizebuffer）;%存储基波频率的缓存
harmonic_all_buffer=zeros（1Sizebuffer）;%总谐波比例
harmonic_even_buffer=zeros（1Sizebuffer）;%偶次谐波比例
harmonic_odd_buffer=zeros（1Sizebuffer）;%奇次谐波比例

for i=1:Sizebuffer
    
    %用采集卡采集电压数据
    %采样结果放到数组x中，x的数据点数为Nsample
    %TODO: 添加从采集卡中读取数据程序
    
    
    
    %软件测试时用matlab产生测试数据
    %电压信号的参数
    f1=50+randn（）*0.3;%基波的频率，5049.650.1
    %谐波（包含基波和高次谐波）的频率、幅度有效值（V）、初始相位（角度）
    f_HA_groundtruth=f1*（1:9）;
    U1=220+randn（）*5;
    Uharmonic=[5.33.2291.1810.5];%谐波的幅度
    Uharmonic（1:2:end）=Uharmonic（1:2:end）*rand（）*1;%偶次谐波乘以系数，*0.5 *1 *3
    Uharmonic（2:2:end）=Uharmonic（2:2:end）*rand（）*1;%奇次谐波乘以系数，*0.5 *1 *3
    A_HA_groundtruth=[U1Uharmonic];
    phi_HA_groundtruth=[010983086088];%有谐波信号
    %A_HA_groundtruth=[22000000000];phi_HA_groundtruth=[000000000];%标准信号
    %生成采样数据
    t=（0:Nsample-1）/fs;%采样时刻数组
    x_contious=zeros（1Nsample）;
    for i_harmonic=1:numel（f_HA_groundtruth）
        x_contious=x_contious+A_HA_groundtruth（i_harmonic）*sqrt（2）*cos（f_HA_groundtruth（i_harmonic）*2*pi*t+phi_HA_groundtruth（i_harmonic）/180*pi）;
    end
    %建模ADC采样的幅度离散化
    ADCvalue=round（x_contious/kADC）;
    x=ADCvalue*kADC;
    
    
    

    %谐波分析
    [f_HAA_HAphi_HA]=harmonic_analysis（xfsN_HA）;
    %比较计算结果
    %[f_HA;A_HA;phi_HA]
    %[f_HA_groundtruth;A_HA_groundtruth;phi_HA_groundtruth]
    U1buffer（i）=A_HA（1）;
    f1buffer（i）=f_HA（1）;
    harmonic_all_buffer（i）=norm（A_HA（2:end））/A_HA（1）;
    harmonic_even_buffer（i）=norm（A_HA（2:2:end））/A_HA（1）;
    harmonic_odd_buffer（i）=norm（A_HA（3:2:end））/A_HA（1）;
    
    %结果分析
    figure（1）
    subplot（231）
    plot（U1buffer（1:i））;
    ylabel（‘基波电压有效值（V）‘）;
    subplot（232）
    plot（f1buffer（1:i））;
    ylabel（‘基波频率（Hz）‘）;
    subplot（234）
    plot（harmonic_all_buffer（1:i）*100）;
    ylabel（‘总谐波比例（%）‘）;
    subplot（235）
    plot（harmonic_even_buffer（1:i）*100）;
    ylabel（‘偶次谐波比例（%）‘）;
    subplot（236）
    plot（harmonic_odd_buffer（1:i）*100）;
    ylabel（‘奇次谐波比例（%）‘）;
    
    pause（0.5）;%延时，使用数据采集卡时延时可以长些，模拟产生测试数据时减小延时
end


end


%电力谐波分析函数[f_HAA_HAphi_HA]=harmonic_analysis（xfsN_HA）
%根据论文基于加汉宁窗的FFT高精度谐波检测改进算法实现。
%在中心谱线估计