[SciPy-User] Boxcar smoothing of 1D data array...?

[tho...@nbi.ku.dk]

Hello list;

This seems like it should be a simple task, but I couldn't seem to find
anything in the docs about it - or rather, what I found seems to be from
the Numeric/Numarray days and not valid anymore.

As the subject line suggests, I have a 1D array that I want to
smooth/convolve with a Boxcar kernel of a certain width. In IDL there's
simply a function to do this, and there might be something people have
hacked together out there to do it too - but isn't there a simple way to
do it using built-in NumPy and SciPy tools?

Cheers;
Emil

_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[Benjamin Root]

Emil,

You can find various windowing functions like boxcar, hamming and hanning in scipy.signals module.

Ben Root

[David Baddeley]

Alternatively you could just use scipy.convolve with a tophat kernel ie (for a filter of length N & signal y):

scipy.convolve(y, ones(N)/N)

see the docs for scipy.convolve for more info (you might want to specify how it handles the ends, for example)

cheers,

David

---

From: Benjamin Root <ben....@ou.edu>
To: SciPy Users List <scipy...@scipy.org>
Sent: Wed, 16 June, 2010 3:46:37 AM
Subject: Re: [SciPy-User] Boxcar smoothing of 1D data array...?

[Davide Lasagna]

You could use the Savitzky-Golay smoothing filter function present in the cookbook. It is very well suited for such operation.

My two cents,

Davide

[Sturla Molden]

Den 16.06.2010 02:26, skrev David Baddeley:

Alternatively you could just use scipy.convolve with a tophat kernel ie (for a filter of length N & signal y):

scipy.convolve(y, ones(N)/N)

see the docs for scipy.convolve for more info (you might want to specify how it handles the ends, for example)

You should not use convolution for boxcar filtering. It can be solved using a recursive filter, basically

    y[n] = y[n-1] + x[n] - x[n-m]

then normalize y by 1/m.

Sturla

[Anne Archibald]

How does the numerical stability of this compare to a FIR
implementation (with or without a Fourier transform)?

Anne

> Sturla

[David Baddeley]

Out of curiosity, are there any reasons other than performance (which might be moot if you have to implement the recursive filter as a python loop) for not using a convolution?

cheers,

David

---

From: Sturla Molden <stu...@molden.no>
To: scipy...@scipy.org
Sent: Fri, 18 June, 2010 11:22:45 PM

Subject: Re: [SciPy-User] Boxcar smoothing of 1D data array...?

[Sturla Molden]

Den 18. juni 2010 kl. 16.51 skrev Anne Archibald <aarc...@physics.mcgill.ca
>:

>>
>>
>>
>> y[n] = y[n-1] + x[n] - x[n-m]
>>
>> then normalize y by 1/m.
>
> How does the numerical stability of this compare to a FIR
> implementation (with or without a Fourier transform)?
>
>>

For practical purposes, x will be a digital signal (from an ADC) or a
digital image. Thus the recursive boxcar can be implemented with
integer maths. Stability is excellent as numerical error is 0. :-)

You just have to make sure that y does not overflow (e.g. let y be 32
bit if x is 16 bit).

[Anne Archibald]

On 19 June 2010 06:00, Sturla Molden <stu...@molden.no> wrote:
>
>
> Den 18. juni 2010 kl. 16.51 skrev Anne Archibald <aarc...@physics.mcgill.ca
>  >:
>
>>>
>>>
>>>
>>>     y[n] = y[n-1] + x[n] - x[n-m]
>>>
>>> then normalize y by 1/m.
>>
>> How does the numerical stability of this compare to a FIR
>> implementation (with or without a Fourier transform)?
>>
>>>
>
> For practical purposes, x will be a digital signal (from an ADC) or a
> digital image. Thus the recursive boxcar can be implemented with
> integer maths. Stability is excellent as numerical error is 0. :-)
>
> You just have to make sure that y does not overflow (e.g. let y be 32
> bit if x is 16 bit).

Heh. You have a point there. But I should say that in the application
in which we use boxcar filtering (searching for single pulses in radio
pulsar search data), the data has already been processed sufficiently
that we can't use integers any more, and in fact we use 32-bit floats
rather than doubles. It's kind of moot for us in any case since we
plan to modify the code to do matched filtering with a different
filter, so convolution will be necessary.

It's also worth checking: while scipy.signal does implement IIR
filters, I don't think it takes advantage of zero coefficients to
avoid arithmetic, so using it to implement a boxcar is probably worse
than using even a non-FFT convolution. Is this right?

Anne

[Charles R Harris]

Note that the transfer functions are the same if the zero at 1 in the numerator exactly cancels the zero in the denominator. However, the output has to be correctly initialized, since the initial value is otherwise carried along, which you can see by setting the inputs to zero. So the algebraic cancellation doesn't remove all the side effects of using a recursive filter.

Chuck