Понимание вывода БПФ

87

Мне нужна помощь в понимании результатов вычисления ДПФ / БПФ.

Я опытный инженер-программист, и мне нужно интерпретировать некоторые показания акселерометра смартфона, например, найти основные частоты. К сожалению, пятнадцать лет назад я проспал большую часть уроков EE в колледже, но последние несколько дней я читал о ДПФ и БПФ (очевидно, без особой пользы).

Пожалуйста, никаких ответов на "иди и возьми урок EE". Я действительно планирую это сделать, если мой работодатель будет мне платить. :)

Итак, вот моя проблема:

Я записал сигнал с частотой 32 Гц. Вот 1-секундный образец из 32 точек, который я нанес в Excel.

введите описание изображения здесь

Затем я получил код БПФ, написанный на Java из Колумбийского университета (после того, как я выполнил предложения в сообщении « Надежное и быстрое БПФ на Java »).

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

Before: 

Re: [0.887  1.645  2.005  1.069  1.069  0.69  1.046  1.847  0.808  0.617  0.792  1.384  1.782  0.925  0.751  0.858  0.915  1.006  0.985  0.97  1.075  1.183  1.408  1.575  1.556  1.282  1.06  1.061  1.283  1.701  1.101  0.702  ]

Im: [0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ]

After: 

Re: [37.054  1.774  -1.075  1.451  -0.653  -0.253  -1.686  -3.602  0.226  0.374  -0.194  -0.312  -1.432  0.429  0.709  -0.085  0.0090  -0.085  0.709  0.429  -1.432  -0.312  -0.194  0.374  0.226  -3.602  -1.686  -0.253  -0.653  1.451  -1.075  1.774  ]

Im: [0.0  1.474  -0.238  -2.026  -0.22  -0.24  -5.009  -1.398  0.416  -1.251  -0.708  -0.713  0.851  1.882  0.379  0.021  0.0  -0.021  -0.379  -1.882  -0.851  0.713  0.708  1.251  -0.416  1.398  5.009  0.24  0.22  2.026  0.238  -1.474  ]

Итак, на данный момент я не могу разобраться в результатах. Я понимаю концепции DFT, такие как реальная часть - это амплитуды составляющих косинусоидальных волн, а мнимая часть - это амплитуды составляющих синусоидальных волн. Я также могу проследить эту диаграмму из замечательной книги « Руководство для ученых и инженеров по цифровой обработке сигналов »: введите описание изображения здесь

Итак, мои конкретные вопросы:

  1. Как я могу найти «наиболее часто встречающиеся частоты» на выходе БПФ? Это часть моего анализа данных акселерометра. Что мне следует читать: реальные (косинусные) или мнимые (синусоидальные) массивы?

  2. У меня есть 32-точечный ввод во временной области. Разве результат БПФ не должен быть массивом из 16 элементов для действительных чисел и массивом из 16 элементов для мнимых? Почему программа выдает как вещественные, так и воображаемые массивы размером 32?

  3. В связи с предыдущим вопросом, как мне анализировать индексы в выходных массивах? Учитывая мой ввод 32 отсчетов с частотой дискретизации 32 Гц, я понимаю, что выход массива из 16 элементов должен иметь свой индекс, равномерно распределенный до 1/2 частоты дискретизации (32 Гц), поэтому правильно ли я понимаю, что каждый элемент массива представляет (32 Гц * 1/2) / 16 = 1 Гц?

  4. Почему на выходе БПФ отрицательные значения? Я думал, что значения представляют собой амплитуды синусоиды. Например, выход Real [3] = -1,075 должен означать амплитуду -1,075 для косинусоидальной волны с частотой 3. Верно? Как амплитуда может быть отрицательной?

stackoverflowuser2010
источник
Что бы вы хотели вычислить по показаниям акселерометра: скорость, расстояние? Шум показаний акселерометра следует распределению Гаусса, и я не могу понять, как подгонка синусоиды может исправить это.
Али
2
тег java следует удалить, поскольку он более общий, чем для определенного языка
user3791372,
Если посмотреть на источник Колумбийского университета, это совсем не эффективно. Это простая неоптимизированная реализация Cooley-Tucky с поисковыми таблицами бабочек, а инверсия битов выполняется вручную вместо использования существующих библиотечных функций
Марк Джеронимус
@MarkJeronimus: Можете ли вы порекомендовать эффективную реализацию БПФ на Java? Если я правильно помню, причина, по которой я выбрал код Колумбийского университета, заключалась в том, что библиотека FFTW была слишком сложной для запуска на смартфоне Android.
stackoverflowuser2010
Я нашел несколько разрозненных «оптимизированных» реализаций, но в основном они представляют собой один алгоритм на размер N, поэтому, если вам нужен диапазон размеров, вам понадобятся все эти процедуры. На практике я в основном использовал Intel Integrated Performance Primitives (да, от Java до JNA), но это не бесплатно. Дома я использую в основном тот же алгоритм, который вы связали, но написанный с нуля в 2005 году с использованием учебника. Это просто БПФ (быстрое преобразование Фурье), ничего такого «быстрого» в нем, чтобы оправдать название «быстрое БПФ».
Марк Иероним

Ответы:

85
  1. Вам не следует искать действительную или воображаемую часть комплексного числа (это ваш реальный и воображаемый массив). Вместо этого вы хотите найти величину частоты, которая определяется как sqrt (real * real + imag * imag). Это число всегда будет положительным. Теперь все, что вам нужно найти, - это максимальное значение (игнорируйте первую запись в вашем массиве. Это ваше смещение постоянного тока и не несет никакой частотно-зависимой информации).

  2. Вы получаете 32 реальных и 32 мнимых выхода, потому что вы используете комплексное БПФ. Помните, что вы преобразовали свои 32 отсчета в 64 значения (или 32 комплексных значения), расширив его нулевыми мнимыми частями. В результате получается симметричный выходной сигнал БПФ, где результат частоты встречается дважды. После того, как вы будете готовы к использованию на выходах от 0 до N / 2, и после зеркального отображения на выходах от N / 2 до N. В вашем случае проще всего просто игнорировать выходы с N / 2 до N. Они вам не нужны, они просто артефакт того, как вы вычисляете свое БПФ.

  3. Частота для уравнения fft-bin равна (bin_id * freq / 2) / (N / 2), где freq - ваша частота дискретизации (также известная как 32 Гц, а N - размер вашего БПФ). В вашем случае это упрощается до 1 Гц на ячейку. Бины с N / 2 по N представляют отрицательные частоты (я знаю, странная концепция). Для вашего случая они не содержат какой-либо важной информации, потому что они просто зеркало первых N / 2 частот.

  4. Ваша действительная и мнимая части каждого бункера образуют комплексное число. Ничего страшного, если действительная и мнимая части отрицательны, а величина самой частоты положительна (см. Мой ответ на вопрос 1). Я предлагаю вам прочитать комплексные числа. Объяснение того, как они работают (и почему они полезны), выходит за рамки того, что можно объяснить в одном вопросе stackoverflow.

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

Нильс Пипенбринк
источник
1
Спасибо. По поводу 1: я видел эту страницу Matlab, которая показывает частотный спектр ( mathworks.com/help/techdoc/ref/fft.html ). На этой странице находится график с заголовком «Односторонний амплитудный спектр y (t)». Это график величины частоты, как вы предложили, sqrt (real ^ 2 + img ^ 2)? Относительно 3: я все еще не получаю результат 2 Гц на ячейку. В моем случае N = 32 и freq = 32, верно? Итак, имеется N / 2 = 32/2 = 16 бинов, а самая высокая частота (Найквиста) равна freq / 2 = 32/2 = 16 Гц, что дает 16 Гц на 16 бинов, что дает 1 Гц на интервал?
stackoverflowuser2010
1
Да, график показывает величину спектра - | Y (f) |. Полоски абсолютного значения означают величину. Ширина бина = частота дискретизации / размер БПФ. Ваша частота дискретизации - 32 Гц, размер БПФ - 32. Да, вы правы насчет ширины бина!
Мэтт Монтаг
Исправлена ​​частота бункера.
Андре Чалелла
1
Хороший ответ, спасибо! Извините, что я немного опаздываю на вечеринку, но, возможно, вы могли бы мне ответить, в какой единице обычно используется величина частоты (как упомянуто вами в пункте 1)? В моем случае по сигналу значений от акселерометра (м / с ^ 2). Я не могу понять этого.
Маркус Вюстенберг
Очаровательно! Полосы частот моей визуализации музыки выходили зеркально слева направо; ответ №2 объясняет это !! Псих!!
Ryan S
11

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

Пол Р
источник
Я ценю, что после этого ответа прошло довольно много времени ... Однако не могли бы вы уточнить, какие артефакты вы имеете в виду?
MattHusz
1
@MattHusz: общий термин, обозначающий происхождение этих артефактов - «спектральная утечка» - я добавил ссылку на ответ, которая объясняет это. Однако лучший способ описать эффект - это то, что ваш спектр будет «размазан» из-за неявного прямоугольного окна.
Paul R
6

1) Ищите в реальном массиве индексы с самыми высокими значениями, помимо первого (это компонент DC). Вероятно, вам понадобится частота дискретизации, значительно превышающая 32 Гц, и больший размер окна, чтобы получить значимые результаты.

2) Вторая половина обоих массивов является зеркалом первой половины. Например, обратите внимание, что последний элемент реального массива (1.774) совпадает со вторым элементом (1.774), а последний элемент мнимого массива (1.474) является отрицательным для второго элемента.

3) Максимальная частота, которую вы можете выбрать при частоте дискретизации 32 Гц, составляет 16 Гц ( предел Найквиста ), поэтому каждый шаг составляет 2 Гц. Как отмечалось ранее, помните, что первый элемент равен 0 Гц (т. Е. Смещение постоянного тока).

4) Конечно, отрицательная амплитуда имеет смысл. Это просто означает, что сигнал "перевернут" - стандартное БПФ основано на косинусе, который обычно имеет значение = 1 при t = 0, поэтому сигнал, который имел значение = -1 в момент времени = 0, будет иметь отрицательную амплитуду. .

сумрак-неактивный-
источник
Спасибо за ответ. (1) Вы имеете в виду, что я могу игнорировать воображаемый (синусоидальный) массив, и если да, то почему? Несомненно, важна синусоида? (2) Почему происходит это зеркальное отображение? Это просто результат алгоритма БПФ? Большинство людей игнорируют зеркальную половину? (3) Как вы рассчитали шаг 2 Гц? Я понимаю предел Найквиста в 16 Гц, поэтому, если есть 16 (не зеркальных) элементов массива, каждый элемент должен быть 16 Гц / 16 = 1 Гц каждый? (4) Чтобы найти основные частоты, нужно ли мне просто взять абсолютное значение значений амплитуды в выходных массивах?
stackoverflowuser2010
Вы не должны искать в реальном массиве максимальное значение, и вы не можете игнорировать массив синус / I. Вместо этого вам нужна величина составного комплексного вектора. Зеркальное отображение происходит потому, что половина входных данных (массив I) - это все нули, поэтому результат имеет половину степеней свободы. Вы можете игнорировать это, если ваши данные строго реальны.
hotpaw2
@duskwuff Большое спасибо: вы ответили на вопрос, который я собирался опубликовать, если бы я не нашел ваш ответ: как интерпретировать ВТОРОЙ части БПФ. Я хочу изменить данные и выполнить обратное, и я продолжал получать только половину результатов, потому что я изменил неправильные данные в этой части. Еще раз спасибо.
Martin
(3), значение step = 2Hz остается для меня неявным. У нас есть 16 ячеек, представленных массивом длины = 16. Нам нужно описать все частоты от 0 Гц до 16 Гц. Я предполагаю, что каждая корзина описывает часть этого диапазона, не так ли?
Krafter
@krafter Я думаю, что это уменьшилось вдвое, потому что вы не можете вывести частоту из одного значения (поскольку нет повторений).
JVE999,
5

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

hotpaw2
источник