Фильтрация пространства Фурье - PullRequest
Новая
       50

Фильтрация пространства Фурье

5 голосов
/ 23 сентября 2010

У меня есть действительный векторный временной ряд x длины T и фильтр h длины t << T. h - фильтр в пространстве Фурье, вещественный и симметричный.Это примерно 1 / f. </p>

Я хотел бы отфильтровать x с помощью h, чтобы получить y.

Предположим, что t == T и БПФ длины T могут поместиться в память (ни один из которых не являетсяправда).Чтобы получить мой отфильтрованный x в python, я бы сделал: 

import numpy as np
from scipy.signal import fft, ifft

y = np.real( np.ifft( np.fft(x) * h ) ) )

Поскольку условия не выполняются, я попробовал следующий метод взлома:

  1. Выберите размер отступа P
  2. Масштаб h с помощью сплайн-интерполяции был бы размером B + 2P> t (h_scaled)
  3. y =[];Цикл:
    • Возьмите блок длиной B + 2P из x (называемый x_b)
    • Выполните y_b = ifft (fft (x_b) * h_scaled)
    • Отбросьте отступ P из любогосторона y_b и объединить с y
    • Выбрать следующее перекрытие x_b с последним на P

Это хорошая стратегия?Как правильно выбрать отступ P?Как правильно это сделать?Я не знаю много обработки сигналов.Это хороший шанс для изучения.

Я использую cuFFT, чтобы ускорить процесс, и было бы здорово, если бы основная часть операций была FFT.Актуальная проблема - 3D.Кроме того, меня не беспокоят артефакты из акаузального фильтра.

Спасибо, Пол.

1 Ответ

6 голосов
/ 23 сентября 2010

Вы на правильном пути.Техника называется обработка с сохранением перекрытия .Является ли t достаточно коротким, чтобы БПФ такой длины помещались в памяти?Если это так, вы можете выбрать размер блока B, такой, что B > 2*min(length(x),length(h)), и обеспечивает быстрое преобразование.Затем при обработке вы отбрасываете первую половину y_b, а не с обоих концов.

Чтобы понять, почему вы отбрасываете первую половину, помните, что спектральное умножение такое же, как круговая сверткаобласть времени.Свертывание с добавлением нуля h создает странные сбои в первой половине результата, но ко второй половине все переходные процессы исчезли, потому что точка циклического переноса в x выровнена с нулевой частью h,Этому есть хорошее объяснение, с фотографиями, в «Теория и применение цифровой обработки сигналов» Лоуренса Рабинера и Бернарда Голда .

Важно, чтобы ваш фильтр во временной области сузился до 0 припо крайней мере, на одном конце, чтобы вы не получили звона артефактов.Вы упоминаете, что h является действительным в частотной области, что означает, что он имеет все 0 фаз.Обычно такой сигнал будет непрерывным только по кругу, и при использовании в качестве фильтра будет создавать искажения во всем диапазоне частот.Один простой способ создать приемлемый фильтр - это создать его в частотной области с 0 фазой, обратным преобразованием и поворотом.Например: 

def OneOverF(N):
    import numpy as np
    N2 = N/2; #N has to be even!
    x = np.hstack((1, np.arange(1, N2+1), np.arange(N2-1, 0, -1)))
    hf = 1/(2*np.pi*x/N2)
    ht = np.real(np.fft.ifft(hf)) # discard tiny imag part from numerical error
    htrot = np.roll(ht, N2)
    htwin = htrot * np.hamming(N)
    return ht, htrot, htwin

(я довольно плохо знаком с Python, пожалуйста, дайте мне знать, если есть лучший способ кодировать это).

Если вы сравните частотные характеристики ht, htrot и htwin, вы видите следующее (ось X - нормализованная частота до pi): 1/f filters with length 64

ht, вверху, имеет много пульсаций,Это связано с разрывом на краю.htrot, в середине, лучше, но все еще имеет пульсацию.htwin красиво и плавно, за счет выравнивания на чуть более высокой частоте.Обратите внимание, что вы можете увеличить длину отрезка прямой линии, используя большее значение для N.

Я писал о проблеме разрыва, а также написал пример Matlab / Octave в еще один вопрос SO , если вы хотите увидеть больше деталей.

...