我有 1 个时间单位的信号历史记录。我的主导频率是 1/100 时间单位。当我使用 numpy 的 fft 函数时,我的分辨率受到信号历史范围的限制。 如何在不破坏信号的情况下提高频率梳的分辨率?
import numpy as np
import matplotlib.pyplot as plt
'''
I need to caputre a low-frequency oscillation with only 1 time unit of data.
So far, I have not been able to find a way to make the fft resolution < 1.
'''
timeResolution = 10000
mytimes = np.linspace(0, 1, timeResolution)
mypressures = np.sin(2 * np.pi * mytimes / 100)
fft = np.fft.fft(mypressures[:])
T = mytimes[1] - mytimes[0]
N = mypressures.size
# fft of original signal is limitted by the maximum time
f = np.linspace(0, 1 / T, N)
filteredidx = f > 0.001
freq = f[filteredidx][np.argmax(np.abs(fft[filteredidx][:N//2]))]
print('freq bin is is ', f[1] - f[0]) # 1.0
print('frequency is ', freq) # 1.0
print('(real frequency is 0.01)')
我认为我可以通过端到端粘贴信号并执行 fft 来人为地增加时间历史长度(从而减少频率梳的宽度)。由于某些我不明白的原因,这对我不起作用:
import numpy as np
import matplotlib.pyplot as plt
timeResolution = 10000
mytimes = np.linspace(0, 1, timeResolution)
mypressures = np.sin(2 * np.pi * mytimes / 100)
# glue data to itself to make signal articicially longer
timesby = 1000
newtimes = np.concatenate([mytimes * ii for ii in range(1, timesby + 1)])
newpressures = np.concatenate([mypressures] * timesby)
fft = np.fft.fft(newpressures[:])
T = newtimes[1] - newtimes[0]
N = newpressures.size
# fft of original signal is limitted by the maximum time
f = np.linspace(0, 1 / T, N)
filteredidx = f > 0.001
freq = f[filteredidx][np.argmax(np.abs(fft[filteredidx][:N//2]))]
print('freq bin is is ', f[1] - f[0]) # 0.001
print('frequency is ', freq) # 1.0
print('(real frequency is 0.01)')
最佳答案
您的目标是从“太短”(即 << sample_rate/frequency_of_interest, window)中恢复光谱信息,这似乎雄心勃勃。
即使在最简单的情况下(干净的正弦波,您的示例),数据看起来也很像一条直线(下方左面板)。只有在去除趋势之后,我们才能看到一点点曲率(下面的右面板,注意非常小的 y 值),这就是所有假设算法都可以通过的。特别是,据我所知,FT 将不起作用。
如果我们非常幸运,有一个出路:比较衍生品。
如果你有一个带有偏移量的正弦信号——比如 f = c + sin(om * t
´--- 那么一阶和三阶导数将是 om * cos(om * t)
和 -om^3 * cos(om * t)
´´。
如果信号足够简单和干净,则可以将其与稳健的数值微分一起用于恢复频率 omega。
在下面的演示代码中,我使用 SavGol 滤波器来获取导数,同时去除一些已添加到信号(橙色曲线)中的高频噪声(下面的蓝色曲线)。可能存在其他(更好的)数值微分方法。
样本运行:
Estimated freq clean signal: 0.009998
Estimated freq noisy signal: 0.009871
我们可以看到,在这个非常简单的情况下,频率恢复正常。
使用更多的导数和一些线性分解巫术可以恢复多个频率,但我不打算在这里探讨这个问题。
代码:
import numpy as np
import matplotlib.pyplot as plt
'''
I need to caputre a low-frequency oscillation with only 1 time unit of data.
So far, I have not been able to find a way to make the fft resolution < 1.
'''
timeResolution = 10000
mytimes = np.linspace(0, 1, timeResolution)
mypressures = np.sin(2 * np.pi * mytimes / 100)
fft = np.fft.fft(mypressures[:])
T = mytimes[1] - mytimes[0]
N = mypressures.size
# fft of original signal is limitted by the maximum time
f = np.linspace(0, 1 / T, N)
filteredidx = f > 0.001
freq = f[filteredidx][np.argmax(np.abs(fft[filteredidx][:N//2]))]
print('freq bin is is ', f[1] - f[0]) # 1.0
print('frequency is ', freq) # 1.0
print('(real frequency is 0.01)')
import scipy.signal as ss
plt.figure(1)
plt.subplot(121)
plt.plot(mytimes, mypressures)
plt.subplot(122)
plt.plot(mytimes, ss.detrend(mypressures))
plt.figure(2)
mycorrupted = mypressures + 0.00001 * np.random.normal(size=mypressures.shape)
plt.plot(mytimes, ss.detrend(mycorrupted))
plt.plot(mytimes, ss.detrend(mypressures))
width, order = 8999, 3
hw = (width+3) // 2
dsdt = ss.savgol_filter(mypressures, width, order, 1, 1/timeResolution)[hw:-hw]
d3sdt3 = ss.savgol_filter(mypressures, width, order, 3, 1/timeResolution)[hw:-hw]
est_freq_clean = np.nanmean(np.sqrt(-d3sdt3/dsdt) / (2 * np.pi))
dsdt = ss.savgol_filter(mycorrupted, width, order, 1, 1/timeResolution)[hw:-hw]
d3sdt3 = ss.savgol_filter(mycorrupted, width, order, 3, 1/timeResolution)[hw:-hw]
est_freq_noisy = np.nanmean(np.sqrt(-d3sdt3/dsdt) / (2 * np.pi))
print(f"Estimated freq clean signal: {est_freq_clean:10.6f}")
print(f"Estimated freq noisy signal: {est_freq_noisy:10.6f}")
关于numpy - fft 在短时间历史中找到低频,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/52579415/