python - 在 python 中计算 thd

Question

我正在尝试计算提供的交流电压的总谐波失真值。我使用 Arduino 以超过 8 KHz 的速率对电压数据进行采样，并将这些数据存储到文本文件中。然后我尝试使用以下用 python 编写的代码片段计算 thd:

    import numpy as np
    import scipy.fftpack
    from scipy.fftpack import fft
    from numpy import genfromtxt

    sampled_data = genfromtxt('/../file.txt',delimiter=',')
    abs_yf=np.abs(fft(sampled_data))

    #As far as I know, THD=sqrt(sum of square magnitude of
    #harmonics+noise)/Fundamental value (Is it correct?)So I'm
    #just summing up square of all frequency data obtained from FFT,
    #sqrt() them and dividing them with fundamental frequecy value.

    def thd(abs_data):
    sq_sum=0.0
    for r in range(len(abs_data)):
       sq_sum=sq_sum+(abs_data[r])**2
    sq_harmonics=sq_sum-(max(abs_data))**2.0
    thd=100*sq_harmonics**0.5/max(abs_data)
    return thd

    print "Total Harmonic Distortion(in percent):"
    print thd(abs_yf)

问题是，在我的例子中，获得的 Thd 值在 5% 到 25% 之间变化。 (实际上它实际上不超过 5%)。我究竟做错了什么？还有其他方法可以找出thd吗？

Answer 1

虽然这很安静，但对于像我这样遇到此帖子的任何人:OP 方法存在一些问题。

1) FFT 返回的幅度包括 0 频率仓的幅度，因此如果信号中存在任何直流偏置，则 max(abs_data) 是对应于基频的幅度的假设是不正确的。这是行中的问题

thd = 100*sq_harmonics**0.5 / max(abs_data)

作为快速解决方案，可以忽略与 0 频率相关的振幅。

2) abs_data 的后半部分应该被丢弃，它是前半部分的“镜像”反射。这是由于傅里叶变换的性质。

这两个问题都可以通过更改函数的输入来解决，即通过替换

print thd(abs_yf)

与

print( thd(abs_yf[1:int(len(abs_yf)/2) ]) )

我们已将输入更改为不包含第一个或最后 N/2 个元素。

结果仍然不理想，因为窗口需要正好是整数个周期，正如上面前面的答案所指出的那样。使用带偏移量的纯正弦测试并调整窗口表明该方法工作得很好，但如果出现明显的窗口错误，则会严重失败。

t0=0
tf = 0.02  # integer number of cycles
dt = 1e-4
offset = 0.5
N = int((tf-t0)/dt)
time = np.linspace(0.0,tf,N )    #;

commandSigFreq = 100
Amplitude = 2

waveOfSin = Amplitude*np.sin(2.0*pi*commandSigFreq*time) + offset

abs_yf = np.abs(fft(waveOfSin))
#print("freq is" + str(scipy.fftpack.fftfreq(sampled_data, dt )  ))
#As far as I know, THD=sqrt(sum of square magnitude of
#harmonics+noise)/Fundamental value (Is it correct?)So I'm
#just summing up square of all frequency data obtained from FFT,
#sqrt() them and dividing them with fundamental frequency value.

def thd(abs_data):
    sq_sum=0.0
    for r in range( len(abs_data)):
       sq_sum = sq_sum + (abs_data[r])**2

    sq_harmonics = sq_sum -(max(abs_data))**2.0
    thd = 100*sq_harmonics**0.5 / max(abs_data)

    return thd

print("Total Harmonic Distortion(in percent):")
print(thd(abs_yf[1:int(len(abs_yf)/2) ]))