Take the 2-minute tour ×
Programmers Stack Exchange is a question and answer site for professional programmers interested in conceptual questions about software development. 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 in this image

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.

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

share|improve this question
    
A couple of questions: What is the x argument in your smooth_func() function? And, what exactly does your direct convolution code do about the boundary conditions? –  comingstorm Apr 3 '12 at 19:08
    
x is a 1D array and signal is f(x). I don't know how to implement the boundary conditions properly that's why this problem is arising. –  Omkar Apr 4 '12 at 5:17

2 Answers 2

There are a variety of ways to try to resolve boundary conditions; it's a problem you need to deal with whether you are doing direct convolution or FFT-based convolution. Unfortunately, it's not clear exactly what your difficulty is (your graphs don't have enough context for me to tell where they come from or what they mean); in the following, I will have to guess.

The thing you need to need to remember about FFT-based convolution is that it treats its input as a periodic loop. Your question mentions a problem that "does not arise if we directly integrate to find convolution"; that difference is probably due to this property.

You claim that "the zero padding is responsible for the undesired boundary effects". However, in order for FFT convolution to match the results of direct convolution, you must ensure that there is sufficient zero padding added to the original data to keep the periodic nature of the FFT from interfering with the convolution. For this purpose, I believe that "sufficient zero padding" should be length of the filter you are convolving with. Try splitting the zero padding for your fftconvolve() evenly between the beginning and end of your data run.

If you are unsatisfied with the boundary effects of your direct convolution, I'm not sure what to tell you, since I don't know what your application is. I will advise that the result of your convolution will be as long as the original data plus the length of the filter, so you will need to provide for that.

share|improve this answer
    
Scipy's fftconvolve does the zero-padding for you; I think the problem isn't periodic boundary conditions, it's zero-padded boundary conditions. I bet the OP wants constant-padded boundaries. –  Andrew Jan 3 at 22:50

Use gaussian_filter instead of fftconvolve.

Compare the behavior of fftconvolve (with mode='same') to gaussian_filter (with mode='constant'):

import numpy as np
from scipy.signal import fftconvolve
from scipy.ndimage import gaussian_filter

x = np.linspace(-3, 3, 51)
y = np.sin(x)

blurring_kernel = np.zeros_like(x)
blurring_kernel[25] = 1
blurring_kernel = gaussian_filter(blurring_kernel, sigma=3)

a = fftconvolve(y, blurring_kernel, mode='same')
b = gaussian_filter(y, sigma=3, mode='constant')
print max(abs((a - b))) # a and b are identical
"""
You have to assume *something* outside the boundary.
This is probably what you want:
"""
c = gaussian_filter(y, sigma=3, mode='reflect')

You have to assume something outside the boundary of your signal. fftconvolve assumes zeros. gaussian_filter lets you choose from several different assumptions, and I find one of these is usually closer to my needs than assuming zeros.

For a quick fix, you could use gaussian_filter, or else pad your signal with something nonzero, to get the same effect at the boundary, perhaps using pad.

share|improve this answer

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.