Спектральная плотность мощности от jTransforms DoubleFFT_1D - PullRequest
12 голосов
/ 16 февраля 2011

Я использую 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 прямое ДПФ реальных данных, оставляя результат в.Физическая структура выходных данных выглядит следующим образом:

если 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 на выходе этого метода.

Параметры: a - данные для преобразования

Теперь используем методы, описанные выше: (поскольку длина моего массива данных равна 10, используются методы "n is even")

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 - реальные и мнимые компоненты в выводе FFT.для бункера 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 Гц

.

Ответы [ 3 ]

9 голосов
/ 16 февраля 2011

Я думаю, что вам нужно интерпретировать выходные данные следующим образом:

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

[Извините, что теперь у меня есть два отдельных ответа - я думаю, что два связанных вопроса были объединены, пока я работал над новым ответом]

1 голос
/ 14 февраля 2011

Каждая запись в преобразовании представляет (комплексную) величину частоты в выборке.

Плотность мощности в данной частоте является просто величиной комплексной амплитуды преобразования в этой частоте.величина комплексного числа вычисляется из его компонентов, и у вас не должно возникнуть проблем с получением этого

. Каждый столбец представляет амплитуды для возрастающих частот, начиная с 0 (первая запись), затем 2 Pi / T (гдеT - длина вашей выборки), до последней выборки 2 * Pi * N / T (где N - количество выборок)

. Существуют другие соглашения, в которых преобразование возвращается для -Pi * N./ T частота до Pi * N / T, а 0 частотная составляющая находится в середине массива

надеюсь, это поможет

0 голосов
/ 14 февраля 2011

Чтобы получить величину bin k, вам нужно вычислить sqrt (re * re + im * im), где, im, im - реальные и мнимые компоненты в выходе FFT для bin 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]))
...