aabzel 23 июл 2023 в 04:04

Ряд Фурье как Фильтр Нижних Частот

Простой

4 мин

7.8K

Математика*Визуализация данных*Промышленное программирование*Производство и разработка электроники*Физика

Туториал

Спустя 10 лет опыта работы по специальности после завершения института настал тот самый первый случай, когда понадобилось вспомнить преобразование Фурье из курса спец разделов математического анализа.

Допустим у нас есть сигнал с очень зашумленной высокочастотной составляющей. Для конкретики пусть это будет график естественной освещенности от времени. Вот например такой:

При этом интерес представляет самая низкочастотная составляющая этого сигнала. Та, которая с периодом около 24х часов. Надо определить амплитуду, фазу и желательно точную частоту самой первой гармоники.

Вот этот сигнал измерен реальным физическим датчиком освещенности BH1750 и записан в текстовый *.csv файл.

...
1989, 1410.000,  17.75, 10:12:39, 13/7/2023, 1517845015
1990, 1413.333,  17.75, 10:12:59, 13/7/2023, 1517845035
1991, 1415.833,  17.75, 10:13:19, 13/7/2023, 1517845055
1992, 1420.833,  17.75, 10:13:39, 13/7/2023, 1517845075
1993, 1424.167,  17.75, 10:13:59, 13/7/2023, 1517845095
...

Ось X-это 6й столбец, Ось Y-2ой столбец. При этом стоит задача выделить фазу этого сигнала. Проще говоря, надо найти момент времени, когда происходит максимум усреднённого сигнала.

Как же выделить фазу зашумленного сигнала?

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

Если мы пропустим сигнал через цифровой фильтр скользящего среднего (FIR фильтр), то выходной сигнал сместится по фазе вправо и реальная фаза окажется искажена. Нужен другой способ.

Из курса высшей математики известно, что любую периодическую и кусочно-непрерывную функцию можно представить в виде суммы синусов разной амплитуды и разных фаз согласно формуле (1)

$f(t)= A_{0}+A_{1}sin(\omega t-\phi_1)+A_{2}sin(2\omega t-\phi_2)+A_{3}sin(3\omega t-\phi_3)+... \;\;\;\;\;\;\;(1)$

Это ряд Фурье. Если повышать количество слагаемых (гармоник), то можно уменьшить среднеквадратическую ошибку между исходным сигналом и функцией составленной при помощи этих синусов.

Если оставить только первые низкочастотные гармоники, то получится, что как-будто сигнал пропустили через абстрактный фильтр нижних частот.

Как же вычислить коэффициенты ряде Фурье?

Надо воспользоваться дискретным преобразованием Фурье (DFT). Это чисто численная процедура. Даешь массив исходного сигнала произвольного размера N, указываешь частоту дискретизации (для данного эксперимента это 20 секунд), делаешь вычисления и получаешь массив комплексных чисел той же размерности N.

bool dft_calc(const double* const signal, uint32_t n_big, double complex* const dft_out, uint32_t out_len,
              double t_big) {
    LOG_INFO(DFT, "Calc %u Samples", n_big);
    bool res = false;
    double measure_interval_s = ((double)n_big) * t_big;
    if(signal) {
        res = true;
        uint32_t k = 0;
        for(k = 0; k < n_big; k++) {
            dft_out[k] = 0.0 + 0.0 * I;
            uint32_t n = 0;
            for(n = 0; n < n_big; n++) {
                dft_out[k] += ((double)signal[n]) * (cos(TWO_PI_VAL * ((double)k * n) / ((double)n_big)) -
                                                     sin(TWO_PI_VAL * ((double)k * n) / ((double)n_big)) * I);
            }
            dft_out[k] = 2.0 * dft_out[k] / ((double)n_big);
        }
    }
    return res;
}

В данной реализации временная сложность DFT составляет O(N^2) где N количество элементов в выборке сигнала.

Вот спектр данного сигнала. Тут период в часах.

Модуль комплексного числа будет амплитудой синуса, аргумент комплексного числа будет фазой синуса. Дале просто подставляем эти числа в формулу (1) и получается результат фильтрации.

Вот несколько разложений в ряд Фурье до первой гармоники. Как видно на глаз фаза определилась достаточно достоверно.

К слову преобразование Фурье это настоящий трудяга в IT. Он отлично подходит для сжатия периодических данных. Вместо огромного файла лога периодического сигнала Вы просто сохраняете массив из нескольких пар вещественных коэффициентов для амплитуды и фазы данной конкретной гармоники ряда Фурье. Далее по этим коэффициентам функция примерно восстанавливается до той степени точности которая вам нужна. Однако происходит потеря качества исходного сигнала. То есть ряд Фурье можно использовать для компрессии экспериментальных данных.

Также преобразование Фурье используется для демодуляции радио сигналов в которых происходит фазовая, частотная и амплитудная модуляции.

Достоинства фильтра низких частот на преобразовании Фурье:

1. Если вычислить DFT, обнулить высокие частоты и вычислить обратное DFT, то не меняется фаза тех низкочастотных гармоник, которые мы оставили.

Фильтр на основе преобразования Фурье отлично подходит для пост обработки на PC.

Недостатки фильтра низких частот на преобразовании Фурье:

Сложно себе представить работу такого фильтра в реальном масштабе времени. Ведь надо много вычислений O(n^2). В идеале хочется вычислять преобразование Фурье для последних N тактов, отбрасывать высокочастотные гармоники и вычислять обратное преобразование Фурье. Это еще O(n) действий. И всё это должно происходить за один такт. Но это запредельно много вычислений по сравнению с тем же FIR фильтром.
Надо заранее знать самую низкую частоту сигнала. Получившаяся формула ряда Фурье будет отрисовывать график сигнала с периодом равным размеру выборки умноженную на период дискретизации. Таким образом преобразование Фурье заводится только "с толкача". Преобразованию Фурье надо на самом первом шаге давать подсказку (период первой гармоники).
Преобразование Фурье не очень подходит для непериодических сигналов. Для непериодических сигналов надо смотреть в сторону Вейвлет преобразований.
Преобразование Фурье показывает вклад частот f=1/T, 2f, 3f. Но преобразование Фурье не покажет вклад частоты 0.5f, 1.5f, 2.5f. Вернее покажет, если вы в два раза увеличите выборку 2T и выполните новое отдельное преобразование обработав при этом в 2 раза больше данных.

Вывод

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

Теория рядов Фурье отлично подходит для пост обработки периодического сигнала, если нужно выделить какие-то отдельные частоты в этом сигнале.

Вообще фильтровать сигналы преобразованием Фурье это нетипичный use-сase. Обычно настраивают цифровые фильтры для фильтрации. Но это только подтверждает цитату Барнса Уоллеса (Barnes Wallis)

Чем меньше Вы знаете о предмете, тем выше Ваша способность для представления оригинальных идей

В общем, вступайте в ряды Фурье! Сходимость! Равенство! Гильбертово пространство!

Словарь

Акроним	Расшифровка
FIR	finite impulse response
DFT	Discrete Fourier Transform
PC	Personal computer
FFT	Fast Fourier transform
FIR	Finite Impulse Response

Links

http://www.mathprofi.ru/ryady_furie_primery_reshenij.html
https://www.youtube.com/watch?v=i7dPkwKHFr0
http://latex.codecogs.com/eqneditor/editor.php
https://code911.top/howto/how-to-run-python-script-from-https://pyneng.readthedocs.io/ru/latest/book/05_basic_scripts/args.html
https://ezgif.com/maker
https://habr.com/ru/articles/269991/
https://habr.com/ru/articles/196374/
https://habr.com/ru/articles/324152/

https://www.youtube.com/watch?v=ds0cmAV-Yek

Только зарегистрированные пользователи могут участвовать в опросе. Войдите, пожалуйста.

Вы использовали в своих программах преобразование Фурье?

75.34% да55

24.66% нет18

Проголосовали 73 пользователя. Воздержались 11 пользователей.

Только зарегистрированные пользователи могут участвовать в опросе. Войдите, пожалуйста.

Вы использовали преобразование Фурье для фильтрации сигнала?

45.95% да34

54.05% нет40

Проголосовали 74 пользователя. Воздержались 10 пользователей.

Теги:

Хабы:

Ряд Фурье как Фильтр Нижних Частот

Вывод

Словарь

Links

Публикации

Истории

Ближайшие события