4
$\begingroup$

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):   
    '''
    x is an 1-D array, 
    sig is the input signal and a function of x.
    '''
    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 the link below.
Boundary effects

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:
jumps in the smooth func

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

Also I tried removing the zero padding in the fft source code and replaced it with a constant value, but it gave undesired results.
Can you suggest a way of removing this zero padding in the scipy fft source code?

$\endgroup$
  • $\begingroup$ So what would it do at the edges? To me it seems like you simply want to keep the original signal at the edges. In image processing you would mitigate the boundary issues by simply preserving the boundary at a half the width of the kernel. Is this not what you want to do with your 1D signal? Do your FFT based convolution and replace the edges with the original signal at half the width of the gaussian kernel. $\endgroup$ – Matt Esch Apr 4 '12 at 9:40
  • $\begingroup$ Dupe of math.stackexchange.com/q/127875/2206 $\endgroup$ – endolith Aug 13 '13 at 18:23
3
$\begingroup$

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

Yes there can be, but it all depends what you want to occur at the boundary. If you don't want it to be smoothed then simply replace the edge with the original signal for a sensible number of samples, depending on the width of your gaussian. You gaussian might be > 0 for say 10 samples, so replace the 5 samples at the start and end of the signal with the original samples.

Can you suggest a way of removing this zero padding in the scipy fft source code?

Yes. Pad it yourself. If you want to consider the wrap around, as if your signal is truly periodic, then fabricate the wraparound effect for a power of 2 in size by carefully padding your own signal. Just as a simple idea, multiply your window by 2 and find the nearest power of 2, then fill the middle of the copy with zeros to pad it out

[1234567890] [12345] - zero pad here - [67890]

Now when you apply the convolution, you should get the same effect for the first window's worth of the result as if the convolution kernel had wrapped around the original signal. To illustrate, take your gaussian and zero pad it at the end, and slide through the convolution to see for yourself that it yields the same thing

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[xgaussianx 00000 000000000000 00000]

so the value for the right most sample comes from this:

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[sianx00000 00000 000000000000 xgaus]

Which is equivalent to what we actually want

[1234567890]
[xsianxgaus] (gaussian shifted by 5 samples and wrapped around)

and the value for the left most sample comes from

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[00000xgaus sianx 000000000000 00000]

Which is also equivalent to what we actually want

[1234567890]
[sianxxgaus] (gaussian shifted by -5 samples and wrapped around)

All this relies on deciding what it is you actually want to happen. If you understand the mathematical process of convolution, and you understand exactly what it is you want to calculate, you can fabricate the conditions under power of 2 scenarios to ensure you get fast computation. Generally speaking most people just opt for the replace because their signal windows are not periodic and have sharp discontinuities at the edges.

$\endgroup$
3
$\begingroup$

I know this is old, but I thought I would add something here. You can smoothly blend the part the filtered and unfiltered segments with a gradual function such as a sinusoidal window function. That will make this a little more agnostic to the type of input function since it will retain the original, unprocessed function at the boundaries yet leave less visually unsavory artifacts where the original and processed segments join. In this particular case, it will result in a perfect match assuming you keep the processed signal at 0 during the boundary transient. You would need to start blending in the processed signal

You window two segments: Y) A segment of data you have processed X) The segment you wish to leave unprocessed with the inverse of the window.

Output = window*X + (1.0-window)*Y One possible window function could be:

window =
0.0 for 0 < n < M

0.5 + 0.5*cos(pi*(n-M)/N) for M <= n < (M+N)

1.0 for (M+N) <= n < S

S = Size of entire Data Set

M = Length of transient (size if filter kernel)

N = Number of samples over which you want to blend X and Y (raise from 0 to 1.0)

n = 0 ... S

That way it will smoothly transition the two functions. This is not necessarily better than just rectangular windowing and appending the two together as has been suggested. This could be considered equivalent of the above considering N-->0. Thinking of it that way, "N" is the blending factor. The larger N, the more gradual the transition form X to Y.

This principle is used by Paul Nasca in the ZynAddSubFX synthesizer filter routine to reduce aliasing due to rapidly changing filter coefficients. Two chunks of audio are processed through filters with two different cutoff frequencies where the second filter must allow the boundary transient to decay before it can be made audible. By gradually blending the audio chunks from first filter (steady state) with the second filter, the transient due to the discontinuous filter coefficients change will be masked by the initial steady state filter until the transient has settled. This is performed over a long enough audio block to minimize the "zipper" noise that you would hear if it were to be stair-stepped (short blending period).

I don't consider this application to be very much different from the audio processing example I used.

$\endgroup$
  • $\begingroup$ Wow! An amazingly detailed answer. $\endgroup$ – Naresh Apr 18 '13 at 4:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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