<p>尽管这是一个长时间的沉默,但对于任何像我这样遇到这篇文章的人来说:OP方法有几个问题。在</p>
<p>1)FFT返回的量值包括0频点的幅值,因此,如果信号中存在任何直流偏压,那么假设max(abs_data)是与基频相对应的幅值是不正确的。这是个问题</p>
<blockquote>
<p>thd = 100*sq_harmonics**0.5 / max(abs_data)</p>
</blockquote>
<p>作为一个快速的解决方案,与0频率相关的振幅可以忽略不计。在</p>
<p>2)下半部分的abs_数据应该抛出,它是第一部分的“镜像”反映。这是由于傅里叶变换的性质。在</p>
<p>这两个问题都可以通过更改函数的输入来解决,即通过替换</p>
<blockquote>
<p>print thd(abs_yf)</p>
</blockquote>
<p>与</p>
<blockquote>
<p>print( thd(abs_yf[1:int(len(abs_yf)/2) ]) )</p>
</blockquote>
<p>其中我们更改了输入,不包括第一个或最后一个N/2个元素。在</p>
<p>结果仍然不理想,因为窗口需要正好是前面提到的答案的整数个周期。使用带偏移量的纯正弦曲线进行测试并调整窗口,结果表明该方法工作良好,但如果窗口出现严重错误,则会严重失败。在</p>
<pre><code>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) ]))
</code></pre>