[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