Базовое использование FFTW 3.3.3 с реальными данными

Я новичок в БПФ, и меня попросили найти способ анализа/обработки определенного набора данных, собранных буровыми установками. В собранных данных много шума из-за движений буровой установки (например, вверх и вниз из-за приливов и волн). Меня попросили очистить собранные данные с помощью FFT=>filtering=>IFFT.

Я использую C++ и библиотеку FFTW 3.3.3.

Пример лучше, чем что-либо еще, так что:

У меня есть БД с, например, селевым потоком (литры в минуту). Селевой поток собирается каждые 5 секунд, для каждого измерения в БД есть временная метка (например, 1387411235).

Таким образом, мои IN_data для моего БПФ представляют собой пару меток времени/селей (например, 1387456630/3955,94, 1387456635/3954,92 и т. д.).

Отображение этих данных действительно выглядит как зашумленный звуковой сигнал, и соответствующие события могут быть замаскированы шумом.

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

Я действительно не знаю, как приступить к этой работе, какую версию процедуры FFTW использовать (c2c, r2c и т. д.), есть ли какая-либо предварительная обработка данных и / или постобработка. . Есть много примеров и руководств, которые я читал в Интернете, но я француз (извините за мои ошибки здесь), и это не всегда имеет для меня смысл, особенно с единицами OUT_data, типом OUT_data, входными и выходными данными. размер массива, работа с окнами (кстати, что это такое), короче говоря, я потерялся...

Я полагаю, что моя проблема была бы довольно простой для тех, кто привык к FFTW, но для меня сейчас это очень сложно.

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

  • Какую процедуру FFTW использовать в обоих направлениях (FFT и IFFT) (какой вид, тип и размер массива для IN_data и OUT_data).
  • Как интерпретировать полученный массив (какие единицы вернет FFTW).

На данный момент краткий пример того, что я сделал:

fftw_plan p;
p  = (fftw_plan)fftw_plan_dft_1d(size,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);

с «in» и «out» как fftw_complex (сложный элемент моего массива In_data установлен в 1 для каждых данных, действительно не знаю, почему, но в учебнике сказано, чтобы сделать это).

Этот код основан на примере, найденном в Интернете, но мое отсутствие знаний/понимания сильно тормозит, и мне было интересно, есть ли здесь кто-нибудь, кто мог бы дать мне объяснения/рабочий процесс/понимание/ссылки о том, как это сделать.

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


person Arnaud    schedule 24.12.2013    source источник
comment
Не извиняйтесь за то, что вы француз ;)   -  person Nobody moving away from SE    schedule 24.12.2013
comment
Я извиняюсь больше за отсутствие навыков/знаний, мой друг... но все равно спасибо...   -  person Arnaud    schedule 24.12.2013
comment
FFTW — довольно сложный пакет для полного новичка DSP — я предлагаю начать с чего-то гораздо более простого, например. KissFFT. Кроме того, вы должны быть уверены, что ваши данные имеют однородную выборку, иначе подход на основе БПФ не будет работать.   -  person Paul R    schedule 24.12.2013
comment
Спасибо Полу за ваш ответ, измерение выполняется каждые 5 секунд, поэтому я предполагаю, что выборка является единообразной (в редких случаях данные не собираются, но в БД всегда есть соответствующая метка времени).   -  person Arnaud    schedule 24.12.2013
comment
ОК, я имел в виду, что для данных с равномерной выборкой между выборками должно быть ровно 5 секунд без пропущенных точек данных, иначе вам понадобится более сложный метод, чем БПФ.   -  person Paul R    schedule 24.12.2013
comment
Я всегда могу продублировать данные из предыдущей записи во входном массиве, если найду отсутствующие точки данных, чтобы унифицировать свой набор данных. Мы не отправляем людей в космос здесь, у меня есть небольшая погрешность, я думаю...   -  person Arnaud    schedule 24.12.2013


Ответы (1)


Это довольно амбициозный проект для тех, кто совсем не знаком с DSP, но вы можете начать с прочтения overlap-add, который, по сути, является методом, который вам нужен для вашего подхода FFT-filter-IFFT к очистке этих данных. Вам также следует посетить сайт DSP StackExchange dsp.stackexchange.com, где теоретические основы и применение фильтрации в частотной области рассматриваются в нескольких похожих статьях. Вопросы и Ответы.

person Paul R    schedule 24.12.2013
comment
Еще раз спасибо, Пол, я знаю, что это амбициозно, но, как однажды сказал один известный человек, мы выбираем идти... не потому, что это легко, а потому, что это трудно. Я рассмотрю KissFFT, но FFTW выглядит очень эффективным, и обработка этого сигнала должна быть такой же, как обработка аудиосигнала (хотя и с гораздо большей частотой дискретизации). Здесь есть много вопросов, касающихся обработки звука, и часто ваши ответы, кажется, помогают потерянным... Могу ли я использовать звуковую нить, чтобы погрузиться в эту проблему с FFTW? - person Arnaud; 24.12.2013
comment
KissFFT также довольно эффективен - для сейсмических данных с его низкой частотой дискретизации вам, вероятно, не нужно слишком беспокоиться об эффективности, поэтому просто заставить что-то работать, вероятно, более приоритетно. Вы даже можете подумать о прототипировании этого сначала в MATLAB (или его бесплатном эквиваленте Octave), так как это будет очень быстро реализовать, и вы многому научитесь в процессе и обретете уверенность, прежде чем переходить к реализации на C/C++ со своим Библиотека БПФ по выбору. - person Paul R; 24.12.2013
comment
Да, я искал kissFFT, и это кажется более простым способом. Я сделаю это, ты опытный, и я должен тебя слушать. - person Arnaud; 24.12.2013
comment
Хороший выбор — у KissFFT также нет потенциально сложных проблем с лицензированием, в отличие от FFTW. И вы всегда можете изменить библиотеку БПФ позже, если решите, что сделали неправильный выбор - это действительно черный ящик, если речь идет о вашей реализации. Также обратите внимание, что автор KissFFT участвует в StackOverflow. - person Paul R; 24.12.2013
comment
И последний вопрос (на сегодня ;-)), должен ли я изучать концепции окон, чтобы реализовать то, что меня просят? - person Arnaud; 24.12.2013
comment
В этом контексте у окон есть два аспекта: во-первых, это выбор размера окна, который представляет собой количество точек, которые вы будете обрабатывать для каждого последующего блока, например. 1024 или 2048 — выбор будет зависеть от необходимого вам разрешения по частоте. Во-вторых, это применение оконной функции (например, Хэмминга), которая обычно требуется при выполнении спектрального анализа, но не требуется для операции частотной фильтрации. - person Paul R; 24.12.2013
comment
Спасибо, Пол, счастливого Рождества тебе и твоим близким. Не волнуйся, я вернусь! ;-) - person Arnaud; 24.12.2013
comment
Спасибо и удачи с вашим проектом - я собираюсь путешествовать до начала января, поэтому, возможно, меня не будет рядом, чтобы ответить на вопросы в течение недели или около того, но я постараюсь наверстать упущенное, когда вернусь. А пока: Joyeux Noël et Bonne Année! ???? - person Paul R; 24.12.2013
comment
Вам необходимо найти баланс между разрешением по времени и разрешением по частоте. Это отчасти связано с неопределенностью Гейзенберга. Вам нужно длинное окно, чтобы быть уверенным в частоте события, но короткое окно, чтобы сохранить уверенность в местонахождении временных событий. -- По этой причине были изобретены преобразование Габора и более поздние вейвлет-преобразования. Использование скользящего окна с функцией непрерывного окна похоже на идею преобразования Габора. -- Вам нужно представить, что в БПФ конец окна обтекает начало, если получится скачок, БПФ будет содержать искусственные дефекты. - person Lutz Lehmann; 28.12.2013
comment
Спасибо LutzL за ваш ответ. Функция, которую я кодирую, на данный момент не работает в реальном времени. Конечный пользователь установит диапазон данных, к которым он хочет применить фильтрацию БПФ. На данный момент мы не хотим, чтобы на наших серверах выполнялась обработка FFT в реальном времени (возможно, позже, когда я буду более знаком со всем этим DSP). Итак, прочитав много о методе добавления с перекрытием (как предложил Пол Р.), я все еще задаюсь вопросом, правильный ли это путь для меня, но я должен признать, что это очень интересно и дает мне более глобальное понимание все поле, с которым я имею дело. - person Arnaud; 31.12.2013
comment
Среднее количество точек, которые я буду обрабатывать, составляет 1024, что соответствует 1,42 часам данных. Я думал начать с окна размером 1024. Я также слышал о преобразовании вейвлетов и нашел много библиотек, доступных в Интернете, но я не уверен, действительно ли они здесь нужны. KissFFT запущен и работает, но, как сказал PAUL R, это просто черный ящик, и настоящая работа до И после FFT/IFFT. Конечно, я все еще борюсь со всеми этими концепциями, но день за днем, пост за постом, я становлюсь все более уверенным и знакомым. - person Arnaud; 31.12.2013
comment
Мне просто нужна небольшая помощь для начала, и я надеюсь, что смогу справиться с этим самостоятельно для остальных. Всех с новым годом, да прибудет с нами сила в 2014 году. - person Arnaud; 31.12.2013