Python改进半波傅氏算法求基波幅值（适用无偶次谐波）

        改进半波傅氏算法能有效滤除故障电流中的直流分量和奇数次谐波，正确提取基波分量，若有偶次谐波误差很大。

"""
改进半波傅里叶算法：提取基波幅值，该例中29
"""
import numpy as np

# 定义常数
w1 = 100 * np.pi #基频f=50
T=2*np.pi/w1 #基波周期
N=64 #一个周期采样点数
Ts = T/N #采样周期
r=30 #直流衰减
# 定义函数故障电流函数i(t)：有直流衰减分量，有基波3次5次谐波，基频f=50
def i(t):
    return 70 * np.exp(-r * t) + 29 * np.sin(w1 * t + np.pi / 6) + 67 * np.sin(3 * w1 * t) + 46 * np.sin(5 * w1 * t)
def cosk(k):
    return np.cos(k*2*np.pi/N)
def sink(k):
    return np.sin(k*2*np.pi/N)

def improved_half_fourier():
    # 数据窗:t1不需从0开始，只要时间间隔为Ts即可
    t1= np.arange(0,T/2,Ts)
    t2=t1+Ts
    t3=t2+Ts
    k=np.arange(0,N/2,1)+1
    ka=np.cos(w1*Ts)
    kb=np.sin(w1*Ts)

    ##实际数据由设备采样数据
    an1=4/N*np.sum(i(t1)*cosk(k))
    bn1=4/N*np.sum(i(t1)*sink(k))

    an2=4/N*np.sum(i(t2)*cosk(k))
    bn2=4/N*np.sum(i(t2)*sink(k))

    an3=4/N*np.sum(i(t3)*cosk(k))
    bn3=4/N*np.sum(i(t3)*sink(k))

    Q=an2-ka*an1+kb*bn1
    R=bn2-ka*bn1-kb*an1
    X=an3-2*ka*an2+an1
    Y=bn3-2*ka*bn2+bn1

    wa=X*(Q*Y-X*R)/(kb*(X**2+Y**2))
    wb=Y*(Q*Y-X*R)/(kb*(X**2+Y**2))

    anc=an1-wa
    bnc=bn1-wb
    I=np.sqrt(anc**2+bnc**2)

    print("估计值：",I)


improved_half_fourier()
"""估计值： 28.99999999999958"""