FPGA Central

Hi all,

I'm trying to postprocess a signal that contains a low-frequency sinusoid
but my DSP skills are mediocre (I've had formal classes in continuous tim
systems). I have a very long data series sampled at 100 Hz, and I kno
that the signal of interest is roughly sinusoidal and in the range of 0.2
to 0.5 Hz (it is a physical oscillation of a plant). The frequency ma
shift slightly within that band over time, and the amplitude definitel
changes. I know there is a small random walk and some white noise, an
possibly some well-separated-in-frequency sinusoidal noise components.

MAIN POINT: I need a good (ideally, best) estimate of the actual value o
the sinuosoid at each sample point, along with an estimate of th
derivative.

I have coded a Savitzky-Golay filter, which essentially fits polynomials
but it is not clear that that is the best model (or maybe I'm choosing th
parameters too haphazardly). The other wrinkle is that perhaps 0.1-1% o
the data points are missing, on average. I believe those sample points ar
taken, but dropped in communications, so the data is, in fact, evenl
sampled at 10ms intervals.

Thanks in advance,
Michael

On 29 Jun, 07:37, "Michael Plante" <michael.pla...@gmail.com> wrote:
> Hi all,
>
> I'm trying to postprocess a signal that contains a low-frequency sinusoid,
> but my DSP skills are mediocre (I've had formal classes in continuous time
> systems). *I have a very long data series sampled at 100 Hz,

What does 'very long' mean?

> and I know
> that the signal of interest is roughly sinusoidal and in the range of 0.25
> to 0.5 Hz (it is a physical oscillation of a plant). *The frequency may
> shift slightly within that band over time, and the amplitude definitely
> changes.

What's the typical time scale of these fluctuations?

>*I know there is a small random walk and some white noise, and
> possibly some well-separated-in-frequency sinusoidal noise components.

Overharmonics of the fundamental or unrelated with
the fundamental? What about aliasing? Are there
anti-alias filters prior to the sampling step?

> MAIN POINT: *I need a good (ideally, best) estimate of the actual valueof
> the sinuosoid at each sample point, along with an estimate of the
> derivative.

What 'value' is this? Frequency? Amplitude? Both?

Depending on the details, there are several possible
approaches:

- Adaptive filters.
- Phase-locked loops.
- Kalman filters.
- Frequency estimators.

You might get more useful answers if you describe both
the process and the goal for the analysis in more detail.

Rune

On Mon, 29 Jun 2009 00:37:38 -0500, Michael Plante wrote:

> Hi all,
>
> I'm trying to postprocess a signal that contains a low-frequency
> sinusoid, but my DSP skills are mediocre (I've had formal classes in
> continuous time systems). I have a very long data series sampled at 100
> Hz, and I know that the signal of interest is roughly sinusoidal and in
> the range of 0.25 to 0.5 Hz (it is a physical oscillation of a plant).
> The frequency may shift slightly within that band over time, and the
> amplitude definitely changes. I know there is a small random walk and
> some white noise, and possibly some well-separated-in-frequency
> sinusoidal noise components.
>
> MAIN POINT: I need a good (ideally, best) estimate of the actual value
> of the sinuosoid at each sample point, along with an estimate of the
> derivative.
>
> I have coded a Savitzky-Golay filter, which essentially fits
> polynomials, but it is not clear that that is the best model (or maybe
> I'm choosing the parameters too haphazardly). The other wrinkle is that
> perhaps 0.1-1% of the data points are missing, on average. I believe
> those sample points are taken, but dropped in communications, so the
> data is, in fact, evenly sampled at 10ms intervals.
>
> Thanks in advance,
> Michael

I was thinking that a plain old bandpass filter from 0.25 to .5Hz might
well be good enough. If throughput is a problem you could prefilter and
decimate, but I couldn't see needing to do that unless you're performing
the computation on some really wimpy embedded processor.

_If_ you can accurately identify the dropped samples you probably want to
fill them in with interpolated data -- with that big of a ratio between
your actual and desired signals a simple average of the surrounding two
points should be plenty good for all but the lowest SNR cases.

Like Rune said, you need to cough up more information. Bandpass
filtering (possibly using an FFT), picking out a signal logically, then
using some math will be a lot easier than trying to implement a Kalman
filter -- so you really want to know if you _must_ do something fancy
_before_ you go investing loads of time into it.

--
www.wescottdesign.com

>On 29 Jun, 07:37, "Michael Plante" <michael.pla...@gmail.com> wrote:
I have a very long data series sampled at 100 Hz,
>
>What does 'very long' mean?

Probably as much as half an hour. The operators tend to get tired afte
some time, and batteries die in a couple hours. However, I have abou
11 hours of data without the sinusoid of interest that I used to find a
Allan variance.


>> and I know
>> that the signal of interest is roughly sinusoidal and in the range o
0.2=
>5
>> to 0.5 Hz (it is a physical oscillation of a plant). =A0The frequenc
may
>> shift slightly within that band over time, and the amplitud
definitely
>> changes.
>
>What's the typical time scale of these fluctuations?

It's hard to say. That's one thing I'm going to look at. In reality, I'
doing this on each of six channels of
manufacturer-temperature-compensated (so they claim) MEMS IMU. Th
readings are from an underwater vehicle, and it is vertically stable, so i
oscillates when disturbed. From just looking at the signal, though, it'
hard to figure out...even when I've done my crude smoothing. It may b
possible to get you the time scale, but, in light of my clarifications
would it be helpful?


>>=A0I know there is a small random walk and some white noise, and
>> possibly some well-separated-in-frequency sinusoidal noise components.
>
>Overharmonics of the fundamental or unrelated with
>the fundamental? What about aliasing? Are there
>anti-alias filters prior to the sampling step?

I suspect the other components are unrelated. Possibly noise from fan
and such. There is an antialias filter at 40 Hz on each rate gyro, and 5
Hz on each accelerometer, and both are, I believe, analog 5-th orde
elliptic filters, but it's a very tiny OTS board, so I only have vague doc
from the manufacturer, as well as documentation on the chips they're using


>> MAIN POINT: =A0I need a good (ideally, best) estimate of the actua
value=
> of
>> the sinuosoid at each sample point, along with an estimate of the
>> derivative.
>
>What 'value' is this? Frequency? Amplitude? Both?

Neither...a cleaned-up time domain signal. Instantaneous amplitude, i
you will. Plus the derivative...


>Depending on the details, there are several possible
>approaches:
>
>- Adaptive filters.
>- Phase-locked loops.
>- Kalman filters.
>- Frequency estimators.
>
>You might get more useful answers if you describe both
>the process and the goal for the analysis in more detail.

Sorry...I wanted to keep the first post short. Again, the actua
frequency doesn't interest me so much as the fact that I suspect I can us
it to do a better job of filtering undesired components. I don't want t
use a Kalman filter because I prefer to keep my smoothing unbiased b
potentially-faulty system models at this point (more likely, I just don'
have appropriate coefficients yet). A phase-locked loop sounds attractive
now that you mention it, though I'll have to go study them some more. M
last attempt at one (for a different app) converged very slowly.





****
Tim Wescott wrote:
>> plain old bandpass filter

My problem with that (and maybe it's just my lack of experience wit
digital filters) is that I thought that would introduce a phase delay.
Again, I am more interested in the value at a given time. To a certai
extent, if it's exactly the same (phase-?) delay over the whole passband
I'd be okay, but I was hoping for some sort of sinusoidal interpolatio
scheme. A relative phase between the channels has physical meaning (e.g.
pseudoforces would be an input to help me figure out where the center o
mass is, I think). I have written dynamic equations, but I wish to clea
up the channels individually first. I would not expect a phase change ove
the passband to matter if it were a purely sinusoidal signal, but it'
probably not a "sharp" frequency peak, regardless of how you process it.

On 29 Jun, 19:19, "Michael Plante" <michael.pla...@gmail.com> wrote:
> >On 29 Jun, 07:37, "Michael Plante" <michael.pla...@gmail.com> wrote:

> >What 'value' is this? Frequency? Amplitude? Both?
>
> Neither...a cleaned-up time domain signal. *Instantaneous amplitude, if
> you will. *Plus the derivative...

Well, if this is a random-walk type signal, heck out the
Durbin-Koopman book on Kalman filters:

http://www.amazon.com/Analysis-Metho...6297738&sr=8-1

Within a couple of pages they come up with a motivational
example of a simple Kalman filter that works on an RW proces.
The rest of the book is then about expanding ondetails and
generalizing that example. One benefit of the Kalman filter
is that it easily handles missing samples.

Apart from that, don't dismiss the bandpass filter quite yet.
True, there might be some delays, but if you are a bit
cautious when designing the filter, you would be able to
correct for them.

Rune

>On 29 Jun, 19:19, "Michael Plante" <michael.pla...@gmail.com> wrote:
>> >On 29 Jun, 07:37, "Michael Plante" <michael.pla...@gmail.com> wrote:
>
>> >What 'value' is this? Frequency? Amplitude? Both?
>>
>> Neither...a cleaned-up time domain signal. =A0Instantaneous amplitude
if
>> you will. =A0Plus the derivative...
>
>Well, if this is a random-walk type signal, heck out the
>Durbin-Koopman book on Kalman filters:
>
>http://www.amazon.com/Analysis-Metho...nce/dp/019852=
>3548/ref=3Dsr_1_1?ie=3DUTF8&s=3Dbooks&qid=3D1246297738& sr=3D8-1

I'm kinda short on money right now, but I have about a half dozen books o
Kalman Filters, most of which I've read (can't quite make it throug
Maybeck, but I've read the other ones, including Gelb, Grewal&Andrews, an
a couple others). I appreciate the suggestion, but, given th
circumstances, is it fair to say I could find it somewhere other tha
Durbin-Koopman? If I took the KF route, there are two major caveats:

1) It would have to be a fixed-interval smoother, not a filter, per se.
This is not an objection, but I do think that got lost in the origina
reaction. (Tim had cautioned about "some really wimpy embedded processor.
This is post-processing on a desktop computer.) I don't want to thro
away the advantage of "noncausality" (if you choose to look at it tha
way...really just a gigantic delay).

2) This is part of the system ID stage, so I don't yet trust my syste
model (circular reasoning, and all that). Consequently, would I model i
as, say, a pair of states giving a (damped for stability?) sinusoid, with
third state representing the instantaneous frequency? I.e., somewha
similar to page 62 of Gelb, but adaptive? Would that even work?


>Within a couple of pages they come up with a motivational
>example of a simple Kalman filter that works on an RW proces.

Just to be clear, the RW is just the drifting bias of the sensor. It i
long term and nearly DC for the time periods and amplitudes of interest.
would add a fourth state to the above for the bias, and choose its proces
noise on the basis of my Allan Variance plots.


>One benefit of the Kalman filter
>is that it easily handles missing samples.

I had previously been averaging adjacent ones to fill them in, as Tim als
suggested (I have timestamped everything. They stick out.), but the KF i
definitely less ad-hoc in that regard.


>Apart from that, don't dismiss the bandpass filter quite yet.
>True, there might be some delays, but if you are a bit
>cautious when designing the filter, you would be able to
>correct for them.

Hrm. Okay, would you say the BPF is the best option mentioned so far?
Can I provide some other sort of information (a particular plot?) tha
would make the decision more clear-cut? I don't know if it's relevant, bu
I neglected to mention that the accelerometer channels show ver
significant quantization (the IMU is marketed for model aircraft, whic
have greater dynamics than my app).

On 29 Jun, 20:57, "Michael Plante" <michael.pla...@gmail.com> wrote:

> >Apart from that, don't dismiss the bandpass filter quite yet.
> >True, there might be some delays, but if you are a bit
> >cautious when designing the filter, you would be able to
> >correct for them.
>
> Hrm. *Okay, would you say the BPF is the best option mentioned so far?

Let's say the BPF or LP is the 'least bad' option so far.
There are a few too many loose ends to suggest other options.
As always, start with the simple, easy-to-use stuff first.
Only when that is insufficient, switch to more heavy-duty
stuff.

PLLs and Kalmans and what have you, might be what is needed
in the end, but it would just be too stupid starting up
those sorts of machineries whithout having at least tried
the simple filters first.

Rune

On Mon, 29 Jun 2009 12:30:36 -0700, Rune Allnor wrote:

> On 29 Jun, 20:57, "Michael Plante" <michael.pla...@gmail.com> wrote:
>
>> >Apart from that, don't dismiss the bandpass filter quite yet. True,
>> >there might be some delays, but if you are a bit cautious when
>> >designing the filter, you would be able to correct for them.
>>
>> Hrm. Â*Okay, would you say the BPF is the best option mentioned so far?
>
> Let's say the BPF or LP is the 'least bad' option so far. There are a
> few too many loose ends to suggest other options. As always, start with
> the simple, easy-to-use stuff first. Only when that is insufficient,
> switch to more heavy-duty stuff.
>
> PLLs and Kalmans and what have you, might be what is needed in the end,
> but it would just be too stupid starting up those sorts of machineries
> whithout having at least tried the simple filters first.
>
> Rune

The only thing that I may add is that if you really do have hours of data
collected at 100Hz, and if your frequencies of interest really only
extends to 0.5Hz, you may wish to consider some easy decimation before
you really start wringing out your data. Even on a PC 720000 data points
can take a while to process...

Having said that, I still think I'd do my initial investigation with
"raw" data -- I'd probably just do an FFT and wait for it, then plot the
spectrum up to a couple of Hz and see what's there.

--
www.wescottdesign.com