Как убрать граничные эффекты, возникающие из-за нулевого заполнения в scipy/numpy fft?
Я сделал код Python, чтобы сгладить данный сигнал, используя преобразование Вейерштрасса, которое в основном является сверткой нормализованного гауссиана с сигналом.
Код выглядит следующим образом:
#Importing relevant libraries
from __future__ import division
from scipy.signal import fftconvolve
import numpy as np
def smooth_func(sig, x, t= 0.002):
N = len(x)
x1 = x[-1]
x0 = x[0]
# defining a new array y which is symmetric around zero, to make the gaussian symmetric.
y = np.linspace(-(x1-x0)/2, (x1-x0)/2, N)
#gaussian centered around zero.
gaus = np.exp(-y**(2)/t)
#using fftconvolve to speed up the convolution; gaus.sum() is the normalization constant.
return fftconvolve(sig, gaus/gaus.sum(), mode='same')
Если я запускаю этот код, скажем, для функции шага, он сглаживает угол, но на границе он интерпретирует другой угол и сглаживает его, что приводит к ненужному поведению на границе. Я объясняю это цифрой, показанной по ссылке ниже.
Граничные эффекты
Эта проблема не возникает, если мы напрямую интегрируем, чтобы найти свертку. Следовательно, проблема не в преобразовании Вейерштрасса, и, следовательно, проблема в функции fftconvolve scipy.
Чтобы понять, почему возникает эта проблема, сначала нам нужно понять, как работает fftconvolve в scipy.
Функция fftconvolve в основном использует теорему свертки для ускорения вычислений.
Короче говоря это говорит:
Свертка (int1,int2)= ОБПФ (FFT (int1)* FFT (int2))
Если мы непосредственно применим эту теорему, мы не получим желаемого результата. Чтобы получить желаемый результат, нам нужно взять значение fft в массиве, в два раза больше max(int1,int2). Но это приводит к нежелательным граничным эффектам. Это потому, что в коде FFT, если размер (int) больше, чем размер (над которым взять FFT), он обнуляет ввод, а затем принимает FFT. Это заполнение нулями именно то, что отвечает за нежелательные граничные эффекты.
Можете ли вы предложить способ устранения этого граничного эффекта?
Я попытался удалить его с помощью простого трюка. После сглаживания функции я сравниваю значение сглаженного сигнала с исходным сигналом около границ, и если они не совпадают, я заменяю значение сглаженного функционала входным сигналом в этой точке.
Это выглядит следующим образом:
i = 0
eps=1e-3
while abs(smooth[i]-sig[i])> eps: #compairing the signals on the left boundary
smooth[i] = sig[i]
i = i + 1
j = -1
while abs(smooth[j]-sig[j])> eps: # compairing on the right boundary.
smooth[j] = sig[j]
j = j - 1
Есть проблема с этим методом, потому что при использовании эпсилона в сглаженной функции есть небольшие скачки, как показано ниже:
прыгает в гладкой функции
Могут ли быть внесены какие-либо изменения в вышеуказанный метод для решения этой краевой задачи?
2 ответа
То, что ядро симметричного фильтра создает на концах, зависит от того, что вы предполагаете, что данные находятся за пределами.
Если вам не нравится внешний вид текущего результата, который принимает нули за пределами обоих концов, попробуйте расширить данные с другим предположением, скажем, отражением данных или продолжением полиномиальной регрессии. Расширьте данные на обоих концах как минимум на половину длины ядра фильтра (за исключением случаев, когда ваше расширение содержит нули, которые предоставляются бесплатно с существующим заполнением нулями, необходимым для некруговой свертки). Затем удалите добавленные конечные расширения после фильтрации и посмотрите, нравится ли вам внешний вид вашего предположения. Если нет, попробуйте другое предположение. Или, еще лучше, используйте фактические данные за гранью, если у вас есть такие.
Лучший подход, вероятно, использовать mode = 'valid'
:
The output consists only of those elements that do not
rely on the zero-padding.
Если вы не можете обернуть свой сигнал или обработанный сигнал является выдержкой из более крупного сигнала (в этом случае: обработайте полный сигнал, а затем обрежьте интересующую область), вы всегда будете иметь краевые эффекты при выполнении свертки. Вы должны выбрать, как вы хотите иметь с ними дело. С помощью mode = valid
просто обрезать их, что является довольно хорошим решением. Если вы знаете, что сигнал всегда "ступенчатый", вы можете расширить фронт и конец обработанного сигнала соответствующим образом.