用python计算thd

Question

我试着计算交流电压的总谐波失真值。我正在使用Arduino以超过8千赫的速率采样电压数据，并将这些数据存储到一个文本文件中。然后，我尝试使用以下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%）。我做错什么了？有没有其他方法可以找出这个问题？在

Answer 1

尽管这是一个长时间的沉默，但对于任何像我这样遇到这篇文章的人来说：OP方法有几个问题。在

1）FFT返回的量值包括0频点的幅值，因此，如果信号中存在任何直流偏压，那么假设max（abs_data）是与基频相对应的幅值是不正确的。这是个问题

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

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

2）下半部分的abs_数据应该抛出，它是第一部分的“镜像”反映。这是由于傅里叶变换的性质。在

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

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) ]))

Answer 2

很可能是测量过程本身增加了额外的失真。在

如果将Arduino ADC与高级测量设备进行比较，Arduino的值很可能会更差。至少你需要一个非常稳定和无抖动的时钟。在

此外，数据的输出（我猜是通过UART）可能会干扰ADC测量的定时。在