FFT найти низкие частоты в короткой истории - PullRequest
Новая
       18

FFT найти низкие частоты в короткой истории

0 голосов
/ 30 сентября 2018

У меня есть 1 единица времени истории сигналов.Моя доминирующая частота составляет 1/100 единиц времени.Когда я использую функцию fft numpy, я ограничен в разрешении размером истории сигналов. Как я могу увеличить разрешение моего частотного гребня, не повреждая мой сигнал? 

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

1 Ответ

0 голосов
/ 30 сентября 2018

Ваша цель - восстановить спектральную информацию из «слишком короткого» окна, т. Е. << sample_rate / interval_of_interest, окно выглядит амбициозным. </p>

Даже в самом простом случае (чистый синусоидальный сигнал, ваш пример) просмотр данныхочень похоже на прямую линию (левая панель внизу).Только после снятия тренда мы можем увидеть небольшую кривизну (справа внизу, обратите внимание на очень маленькие значения y), и это все, что может пройти любой гипотетический алгоритм.В частности, FT --- насколько я вижу - не будет работать.

enter image description here

Если нам очень повезет, есть один путьвыход: сравнение производных.Если у вас синусоидальный сигнал со смещением - как f = c + sin(om * t ´ ---, то первая и третья производные будут om * cos(om * t) и -om^3 * cos(om * t) ´´.Если сигнал достаточно простой и чистый, его можно использовать вместе с надежным числовым дифференцированием для восстановления частоты омега.

В демонстрационном коде ниже я использую фильтр SavGol для получения производных при избавлении от некоторой высокой частотышум (синяя кривая ниже), который был добавлен к сигналу (оранжевая кривая).Могут существовать и другие (лучшие) методы численного дифференцирования.

enter image description here

Пример выполнения: 

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}")
...