在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频率仓的幅度，因此如果信号中存在任何DC偏置，则max(abs_data)是对应于基频的幅度的假设是不正确的。这是线路上的问题

thd = 100*sq_harmonics**0.5 / 最大值(abs_data)

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

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

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

打印 thd(abs_yf)

和

打印( 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) ]))