Спектральная плотность мощности от jTransforms DoubleFFT_1D

StackOverflow https://stackoverflow.com/questions/5010261

  •  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]))
Лицензировано под: CC-BY-SA с атрибуция
Не связан с StackOverflow
scroll top