Единицы частоты при использовании БПФ в NumPy - PullRequest
7 голосов
/ 20 декабря 2011

Я использую функцию FFT в NumPy для обработки сигнала. У меня есть массив с именем signal который имеет одну точку данных для каждого часа и имеет в общей сложности 576 точек данных. Я использую следующий код на signal, чтобы посмотреть на его преобразование Фурье.

t = len(signal)
ft = fft(signal,n=t)
mgft=abs(ft)
plot(mgft[0:t/2+1])

Я вижу два пика, но я не уверен относительно того, что такое единицы измерения оси x, то есть как они отображаются на часы? Любая помощь будет оценена.

Ответы [ 3 ]

11 голосов
/ 20 декабря 2011

Учитывая частоту дискретизации FSample и размер блока преобразования N, вы можете рассчитать разрешение по частоте deltaF, интервал выборки deltaT и общее время захвата capT, используя отношения:

deltaT = 1/FSample = capT/N
deltaF = 1/capT = FSample/N

Помните также, что FFT возвращает значение от 0 до FSample или эквивалентно -FSample/2 до FSample/2.На вашем графике вы уже отбрасываете часть -FSample/2 до 0.NumPy включает вспомогательную функцию для расчета всего этого за вас: fftfreq .

Для ваших значений deltaT = 1 hour и N = 576 вы получаете deltaF = 0.001736 cycles/hour = 0.04167 cycles/day, от -0.5 cycles/hour до 0.5 cycles/hour.Так что, если у вас есть пик величины, скажем, в ячейке 48 (и ячейке 528), который соответствует частотной составляющей в 48*deltaF = 0.0833 cycles/hour = 2 cycles/day.

В общем, вы должны применить оконную функцию к вашим данным во временной области перед вычислением БПФ, чтобы уменьшить спектральную утечку .Окно Ханна почти никогда не является плохим выбором.Вы также можете использовать функцию rfft, чтобы пропустить часть -FSample/2, 0 вывода.Итак, ваш код будет:

ft = np.fft.rfft(signal*np.hanning(len(signal)))
mgft = abs(ft)
xVals = np.fft.fftfreq(len(signal), d=1.0) # in hours, or d=1.0/24 in days
plot(xVals[:len(mgft)], mgft)
0 голосов
/ 21 декабря 2011

Как правило, размерные единицы частоты из БПФ совпадают с размерными единицами частоты дискретизации, относящимися к данным, подаваемым в БПФ, например: на метр, на радиан, в секунду или в вашем случае в час.

Масштабируемые единицы измерения частоты для каждого индекса результирующего бункера FFT равны N / theSampleRate, с такими же размерными единицами, что и выше, где N - длина полного БПФ (вы можете составлять только половину этой длины в случае строго реальных данных).

Обратите внимание, что каждая пиковая ячейка результата БПФ представляет фильтр с ненулевой полосой пропускания, поэтому вы можете добавить некоторые неопределенности или границы ошибок к точкам результата, которые вы отображаете на значения частоты. Или даже используйте метод интерполяционной оценки, если это необходимо и подходит для исходных данных.

0 голосов
/ 20 декабря 2011

Результат преобразования fft отображается не на ЧАСЫ, а на частоты, содержащиеся в вашем наборе данных. Было бы полезно иметь ваш преобразованный график, чтобы мы могли видеть, где находятся пики.

Возможно, у вас есть всплеск в начале преобразованного буфера, так как вы не выполняли никаких операций с окнами.

...