在对信号进行频谱分析时,我遇到了一个奇怪的问题,即所绘制的信号频率发生了偏移(或加倍)。这是一个展示我的方法的简单示例:以100kHz采样1kHz正弦信号。最后,信号仓出现在2kHz而不是1kHz。
import numpy as np
import matplotlib.pyplot as plt
time_step = 1.0/100e3
t = np.arange(0, 2**14) * time_step
sig = np.sin(2*np.pi*1e3*t)
sig_fft = np.fft.rfft(sig)
#calculate the power spectral density
sig_psd = np.abs(sig_fft) ** 2 + 1
#create the frequencies
fftfreq = np.fft.fftfreq(len(sig_psd), d=time_step)
#filter out the positive freq
i = fftfreq > 0
plt.plot(fftfreq[i], 10*np.log10(sig_psd[i]))
plt.xscale("log")
您使用错误的函数来计算频率。实信号的FFT变换具有负频率分量,该负频率分量是正信号的复共轭,即频谱是厄米对称的。 rfft()
利用了这一事实,不输出负频率,仅输出直流分量和正频率。因此,sig_psd
的长度是使用fft()
而不是rfft()
并将其传递给fftfreq()
时得到的结果的两倍。
解决方案:改用rfftfreq()
。
import numpy as np
import matplotlib.pyplot as plt
time_step = 1.0/100e3
t = np.arange(0, 2**14) * time_step
sig = np.sin(2*np.pi*1e3*t)
sig_fft = np.fft.rfft(sig)
#calculate the power spectral density
sig_psd = np.abs(sig_fft) ** 2 + 1
#create the frequencies
fftfreq = np.fft.rfftfreq(len(sig), d=time_step) # <-- note: len(sig), not len(sig_fft)
#filter out the positive freq
i = fftfreq > 0 # <-- note: Not really needed, this just removes the DC component
plt.plot(fftfreq[i], 10*np.log10(sig_psd[i]))
plt.xscale("log")