Спектральная плотность мощности от jTransforms DoubleFFT_1D
-
14-11-2019 - |
Вопрос
Я использую Java-библиотеку Jtransforms для анализа заданного набора данных.
Пример данных следующий:
980,988,1160,1080,928,1068,1156,1152,1176,1264
Я использую функцию DoubleFFT_1D в jTransforms.Вывод данных следующий:
10952, -152, 80.052, 379.936, -307.691, 12.734, -224.052, 427.607, -48.308, 81.472
У меня возникли проблемы с интерпретацией вывода.Я понимаю, что первый элемент выходного массива — это сумма 10 входов (10952).Его
другие элементы выходного массива, которые я не понимаю.В конечном итоге я хочу отобразить спектральную плотность мощности входных данных на графике и найти значения в диапазоне от 0 до 0,5 Гц.
В документации по функциям jTransform указано (где a — набор данных):
public void realForward(double[] a)
Вычисляет 1D вперед DFT реальных данных, оставляя результат в A.Физическая планировка выходных данных заключается в следующем:если n четно, то
a[2*k] = Re[k], 0 <= k < n / 2 a[2*k+1] = Im[k], 0 < k < n / 2 a[1] = Re[n/2]
если n нечетно, то
a[2*k] = Re[k], 0 <= k < (n+1)/2 a[2*k+1] = Im[k], 0 < k< (n-1)/2 a[1] = Im[(n-1)/2]
Этот метод вычисляет только половину элементов реального преобразования.Другая половина удовлетворяет условию симметрии.Если вы хотите полного реального преобразования вперед, используйте Realforwardfull.Чтобы вернуть исходные данные, используйте RealInverse на выводе этого метода.
Параметры:а - данные для преобразования
Теперь, используя методы выше:(поскольку длина моего массива данных равна 10, используются методы «n четно»)
Re[0] = 10952
Re[1] = 80.052
Re[2] = -307.691
Re[3] = -224.052
Re[4] = -48.308
Re[5] = 12.734
Im[0] = -152
Im[1] = 379.936
Im[2] = 12.734
Im[3] = 427.607
Im[4] = 81.472
Итак, несколько вопросов:Этот вывод выглядит корректно?Мне кажется, что Re[0] не должно быть 10952, что является суммой всех элементов исходного массива.
Похоже, вывод следует немного подкорректировать:(Я ошибаюсь?)
Re[0] = 80.052
Re[1] = -307.691
Re[2] = -224.052
Re[3] = -48.308
Re[4] = -152
Im[0] = 379.936
Im[1] = 12.734
Im[2] = 427.607
Im[3] = 81.472
Теперь воспользуюсь следующим методом, опубликованным на форуме:
Чтобы получить величину bin k, вам нужно вычислить sqrt(re * re + im * im)
, где re, im — действительные и мнимые компоненты выходного сигнала БПФ для интервала k.
Для вашего конкретного БПФ re[k] = a[2*k] and im[k] = a[2*k+1]
.Поэтому для расчета спектра мощности:
for k in 0 to N/2 - 1
{
spectrum[k] = sqrt(sqr(a[2*k]) + sqr(a[2*k+1]))
}
Таким образом:
spectrum[0] = 388.278
spectrum[1] = 307.955
spectrum[2] = 482.75
spectrum[3] = 94.717
Некоторые вопросы.Эти данные выглядят корректно?Я на правильном пути?Будут ли эти данные спектра отображать что-то вроде этого:
388.278 at .125 Hz
307.955 at .25 Hz
482.75 at .375 Hz
94.717 at .5 Hz
Я далеко?Моя цель — создать гистограмму спектральной плотности мощности от 0 до 0,5 Гц.
Решение
Я думаю, вам нужно интерпретировать выходные данные следующим образом:
10952 Re[0] = sum of all inputs = DC component
-152 Re[5] - see note about a[1] being special - there is no Im[0]
80.052 Re[1]
379.936 Im[1]
-307.691 Re[2]
12.734 Im[2]
-224.052 Re[3]
427.607 Im[3]
-48.308 Re[4]
81.472 Im[4]
Таким образом, величины таковы:
spectrum[0] = 10952
spectrum[1] = sqrt(80.052^2 + 379.936^2) = 388.278
spectrum[2] = sqrt(-307.691^2 + 12.734^2) = 307.427
spectrum[3] = sqrt(-224.052^2 + 427.607^2) = 482.749
spectrum[4] = sqrt(-48.308^2 + 81.472^2) = 94.717
[Извините, что сейчас у меня есть два отдельных ответа - я думаю, что два связанных вопроса объединились, пока я работал над новым ответом]
Другие советы
Каждая запись в преобразовании представляет (комплексную) величину частоты в выборке.
плотность мощности на данной частоте - это просто величина комплексной амплитуды преобразования на этой частоте.величина комплексного числа вычисляется из его компонентов, и у вас не должно возникнуть проблем с ее получением.
В каждом столбце представлены амплитуды возрастающих частот, начиная с 0 (первая запись), затем 2 Pi/T (где T — длина вашей выборки), до последней выборки 2*Pi*N /T (где N — число образцов)
существуют и другие соглашения, в которых преобразование возвращается для частоты -Pi * N /T до Pi * N / T, а частотный компонент 0 находится в середине массива
надеюсь это поможет
Чтобы получить величину интервала k, вам необходимо вычислить sqrt(re * re + im * im), где re, im — действительные и мнимые компоненты в выходных данных БПФ для интервала k.
Для вашего конкретного БПФ re[k] = a[2*k]
и im[k] = a[2*k+1]
.Поэтому для расчета спектра мощности:
for k in 0 to N/2 - 1
spectrum[k] = sqrt(sqr(a[2*k]) + sqr(a[2*k+1]))