Как сделать оконный полосовой фильтр?

Question

[Antongard]

Здравствуйте.

Есть какой-то сложный сигнал. В нем есть очень много частот (от 0 до 100, например), это можно получить с помощью преобразования Фурье. Но требуется оставить только определенные частоты (например, от 10 до 50) и вывести сигнал обратно только уже с этими определенными частотами.

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

Недавно узнал про оконный полосовой фильтр. Не знаете, как сделать его на Python?

Спасибо)

Comments

[longclaps]

Бред сивой кобылы. Что значит получилось что-то не то? Покажите код.

[Армянское Радио]

А что у вас не так с Фурье-то получилось?

[Antongard]

longclaps, обновил код. Но это преобразование Фурье применяется ко всему сигналу, не разделяя его на части. Меня больше всего интересует именно оконный полосовой фильтр.

[Antongard]

Армянское Радио, обновил код

[Армянское Радио]

FFT = np.absolute(rfft(DATA)) - неправильно. Вы таким образом удалили из сигнала информацию о фазах, вот у вас и получилась фигня.

[Antongard]

Армянское Радио, freqs = rfftfreq(len(data), 1.0/DISCRETE) - вот так?

Однако, мы все-равно делаем Фурье на весь файл. Есть варианты именно с оконным полосовым фильтром?

[Армянское Радио]

Antongard, ну, делайте не на весь файл, а на окно, кто вам мешает? Если не использовать оконную функцию, у вас внесется погрешность, но для начала можно и без нее.

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

Вы какую вообще задачу решаете?

[Antongard]

Армянское Радио, Я разделял данные на части, производил все эти операции, но после этого концы частей не совпадали. Вот я и подумал, что есть другие методы. Для меня погрешность очень важна, так как я работаю с данными электроэнцефалограммы.
Сможете подсказать, как именно сделать оконный полосовой фильтр, ведь, как вы сказали, в нем математика оптимизирована?

[longclaps]

Antongard, а, так ты у нас энцефалограмщик! Щас дам код )))

[res2001]

Antongard,

после этого концы частей не совпадали

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

[longclaps]

res2001, вы не по адресу: я не Antongard.

[res2001]

longclaps, исправился.

[Antongard]

res2001, так как это сделать?

[res2001]

Antongard, Как-то так:
имеем 100 отсчетов сигнала, делим их по 10 отсчетов с 50% перекрытием:
1.начальный индекс: 0, конечный индекс: 10
2.начальный индекс: 5, конечный индекс: 15
3.начальный индекс: 10, конечный индекс: 20
...
Каждую часть прогоняете через вашу обработку, при склеивании результата, не забываете про перекрытие.
В итоге размер входного сигнала будет совпадать с размером результата (если в ходе обработки у вас не увеличиться размерность).

[Antongard]

res2001, спасибо. А в местах перекрытия брать среднее?

[res2001]

Antongard,

брать среднее

Что вы имеете ввиду?
50% перекрытие означает, что в следующей порции сигнала вы берете 50% предыдущей порции сигнала с конца.
20% перекрытия - 20% предыдущего сигнала с конца.
Т.е. это можно использовать как коэффициент приращения начального индекса, схематично так:

siglen = 10;  // длина порции сигнала
overlap = 0.2; // 20% перекрытие
step = (int)ceil(siglen * (1 - overlap)); // шаг приращения индекса
next_begin_idx = prev_begin_idx + step; // индекс следующей порции сигнала (включая)
next_end_idx = next_begin_idx + siglen; // индекс конца порции сигнала (не включая, как принято в C/C++/python)

[Antongard]

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

[res2001]

Antongard, не думаю, что это нужно (брать среднее), но можете поробовать.

[Алексей Узенгер]

Процедура удаления из спектра сигнала ненужных частот с последующим обратным Фурье преобразованием с точки зрения цифровой обработки сигналов является некорректной. Это отмечено у Лайонса в книгах по цифровой обработке сигналов.
Правильнее будет рассчитать полосовой цифровой фильтр (КИХ, БИХ, либо CIC) в математическом пакете, например Matlab (fdatools), и его использовать для обработки входного сигнала.

Answer 1

[longclaps]

import numpy as np, matplotlib.pyplot as plt
from scipy.fftpack import rfft, irfft

n = 2 ** 8
x = np.linspace(0, 2 * np.pi, n)
y = np.sin(x) + np.cos(x * 30)
f0 = rfft(y)
f1 = f0.copy()
f1[:n // 8] = 0
z = irfft(f1)
f, (spectra, signal) = plt.subplots(2, 1, sharey=False)
spectra.plot(x, f0, label='f0')
spectra.plot(x, f1, label='f1')
spectra.legend()
spectra.title.set_text('spectra')
signal.plot(x, y, label='y')
signal.plot(x, z, label='z')
signal.legend()
signal.title.set_text('signal')
plt.legend()
plt.show()

Comments

[Antongard]

Простите, а чем ваше решение отличается от моего?

[longclaps]

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

Answer 2

[Papayaved]

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

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

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

Comments

[Antongard]

Спасибо. А у вас есть пример реализации на python?

Answer 3

[polybook]

Необязательно делать через Фурье-преобразование, можно убрать тренд (НЧ) + сделать скользящее усреднение (сглаживание ВЧ). Параметры НЧ- и ВЧ-фильтра можно контролировать, делая Фурье-преобразование до и после. Подробнее: comma.polybook.ru/1.7.3.pdf