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

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

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

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

1986, 1412.500,  17.75, 10:11:39, 13/7/2023, 1517844955
1987, 1410.833,  17.75, 10:11:59, 13/7/2023, 1517844975
1988, 1410.833,  17.75, 10:12:19, 13/7/2023, 1517844995
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) и получается результат фильтрации.

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

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

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

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

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

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

--Сложно себе представить работу такого фильтра в реальном масштабе времени. Ведь надо много вычислений O (n^2). В идеале хочется вычислять преобразование Фурье для последних N тактов, отбрасывать высокочастотные гармоники и вычислять обратное преобразование Фурье. Это еще O (n) действий. И всё это должно происходить за один такт. Но это запредельно много вычислений по сравнению с тем же FIR фильтром.

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

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

Вывод

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

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

Словарь

Акроним	Расшифровка
FIR	finite impulse response
DFT	Discrete Fourier Transform

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-chttps://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/