Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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?

share|improve this question
    
Don't "zero pad". It's that simple. If you have to pad the boundaries by some value, then use the boundary value. –  Raskolnikov Apr 4 '12 at 7:09
    
Well you see that method dosent work, because if you pad it with some value other than zero, then at the boundaries you will be taking a convolution of a constant and a gaussian, and you get even horrible behaviour at the boundary. –  Omkar Apr 4 '12 at 7:16
    
Why? The stepfunction you are using is essentially constant and convolution with a gaussian just gives a constant when you are far from jumps. Is it because of the discreteness of the fft? –  Raskolnikov Apr 4 '12 at 7:19
    
Ya what you are saying is right, I had tried doing the same thing in the scipy fft source code, but it gave undesired boundaries. But I will try again maybe something else in that code is creating that behaviour. Thanks for the reply –  Omkar Apr 4 '12 at 7:27
2  
You could try posting here: dsp.stackexchange.com –  Vitaliy Kaurov Apr 4 '12 at 8:05
show 1 more comment

1 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
share|improve this answer
    
fftconvolve should produce identical output to convolve, so while using convolve for smaller sizes may speed up the computation, it doesn't solve the OP's problem with boundaries –  endolith Aug 13 '13 at 16:21
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.