[Math] How to remove the boundary effects arising due to zero padding in discrete fft

numerical methodssignal processing

I have made a python code to smoothen a given signal using the Weierstrass transform, which is basically the convolution of a normalised gaussian with a signal.

The code is as follows: 

#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')

If I run this code for say a step function, it smoothens the corner, but at the boundary it interprets another corner and smoothens that too, as a result giving unnecessary behaviour at the boundary. I explain this with a figure shown below.
enter image description here

This problem does not arise if we directly integrate to find convolution. Hence the problem is not in Weierstrass transform, and hence the problem is in the fftconvolve function of scipy.

To understand why this problem arises we first need to understand the working of fftconvolve in scipy.
The fftconvolve function basically uses the convolution theorem to speed up the computation.
In short it says:
convolution(int1,int2)=ifft(fft(int1)*fft(int2))
If we directly apply this theorem we dont get the desired result. To get the desired result we need to take the fft on a array double the size of max(int1,int2). But this leads to the undesired boundary effects. This is because in the fft code, if size(int) is greater than the size(over which to take fft) it zero pads the input and then takes the fft. This zero padding is exactly what is responsible for the undesired boundary effects.

Can you suggest a way to remove this boundary effects?

I have tried to remove it by a simple trick. After smoothening the function I am compairing the value of the smoothened signal with the original signal near the boundaries and if they dont match I replace the value of the smoothened func with the input signal at that point.
It is as follows: 

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

There is a problem with this method, because of using an epsilon there are small jumps in the smoothened function, as shown below:
enter image description here

Can there be any changes made in the above method to solve this boundary problem?

Best Answer

Raskolnikov's suggestion (to replace zero padding by boundary-value padding) is consistent with the mathematical practice, and will eliminate the undesired boundary effect. However, when the arrays have very different sizes (your Gaussian is much narrower than the signal you are smoothing), you may want to ditch fftconvolve in favor of simple convolve, which does not require padding the smaller array. Consider the following running times taken from StackOverflow:

2 arrays, size of $2^{20}$ and $2^4$:

  • numpy.convolve: $110$ ms
  • scipy.signal.convolve: $1.0$ s
  • scipy.signal.fftconvolve: $2.5$ s

2 longer arrays, size of $2^{22}$ and $2^{10}$.

  • numpy.convolve: $6.7$ s
  • scipy.signal.convolve: $221$ s
  • scipy.signal.fftconvolve: MemoryError
Related Question