[SciPy-user] Re: data smoothing: interpolate.splrep ignores s parameter

george young gry at ll.mit.edu
Thu Feb 26 13:20:36 CST 2004

On Thu, 26 Feb 2004 09:15:58 -0500
Paul Barrett <barrett at stsci.edu> threw this fish to the penguins:

> george young wrote:
> > [[reposting with "rough" data file appended, sorry]]
> > [SciPy-0.2.0_alpha_200.4161, Numeric-23.1, Python 2.3.3, x86 Linux]
> > My goal really is to smooth some noisy measurement data without messing
> > up it's *shape*.  My first attempt was 1d splines.  I did:
> > 
> > Splines may not be the right method anyway, since they tend to warp the
> > shape of the curve, and I need to get the data's derivatives.  Is there
> > some way to fashion a low pass filter?  It seems like fft should be useful
> > here, but I have very little experience with fft's.
> > 
> [snip, snip]
> Have you considered wavelet smoothing or denoising?
> The 'a-trous' wavelet method is a simple, understandable technique that can be 
> done by convolving a function with the data.  The difference between your raw 
> data and your smoothed (convolved) data is the wavelet at that particular 
> spatial resolution.  You can filter the wavelet values using thresholding and 
> then add them back into the smoothed data to get your smoothed or denoised 
> result. This can be done at different spatial resolutions depending on the what 
> scale you wish to smooth.
> If you set the filtering threshold correctly, you'll be able to keep the 
> significant features of the data without loosing resolution as is common with 
> averaging techniques, or adding spurious features as is common with FFT 
> techniques.  A B3-spline is good kernel function and usually does a good job in 
> this situation.  Its 1-dimensional representation is [1, 4, 6, 4, 1]/16.

That works very nicely, thanks!  It decreases the noise locally without
warping the overall shape.  I just did: 

    ytmp = zeros(len(data[:, 1]), Float)
    xmin=0; xmax=len(data)-1
    fn= [[-2,1], [-1,4], [0,6], [1,4], [2,1]]
    for i in range(len(data)):
        y = w = 0.0
        for f in fn:
            if xmin <= (i + f[0]) <= xmax:
                y += f[1] * data[i + f[0], 1]
                w += f[1]
        ytmp[i] = y / w

    data[:, 1] = ytmp

I know nothing of wavelets or dsp.  Is there something I could read in a few hours
to understand a little of "a-trous" and "kernels"?  

Is it safe to play with the numbers, e.g. [1,5,8,5,1]/20 or is there 
subtlety that will bite me?

I assume there's functionality in scipy for doing this instead of hard
coding it?  I noticed a convolve function, but not much documentation.

Thanks very much again,
	George Young
"Are the gods not just?"  "Oh no, child.
What would become of us if they were?" (CSL)

More information about the SciPy-user mailing list