Numerical Differentiation in C++

First some info about what I'm trying to do:

Basically I am working with EEG brain signals which are being digitized
through a DAQ device (at around 6Khz) and output in a real-time audio
program (Max/MSP at 48 or 44.1Khz).  

I am trying to code an external object in C++ that will compute time
domain slope descriptor parameters as defined by Horth [1970]

The equations for the parameters are given on page 3 here
http://www.eurasip.org/Proceedings/Eusipco/Eusipco2007/Papers/C1L-E03.pdf

Basically:
1) Activity = variance of input signal x
2) Mobility = sqrt[var(x')/var(x)]
3) Complexity = sqrt[Mobility(x')/Mobility(x)]

I have calculated the variance easily enough but the derivatives are
giving me problems.

At the moment I have implemented the central difference method using the
following approach in the signal loop:

		X = InputSample;
		1stDeriv = (X - PrevPrevX)/(2*(1/SampleRate));
			
		PrevPrevX = PrevX;
		PrevX = X;

		OutputSample = 1stDeriv;

This works well for the 1st derivative (a 10Hz Sin passed through gives a
phase shifted version with amplitude of roughly 10*2*Pi).

However when I pass the function through a second time to get the 2nd
derivative the amplitude ends up being huge! (around 18000!)

Also I've read that high frequency noise can result from this type of
numerical differentiation but his shouldn't be too much of an issue as the
only frequencies of interest to me are < 100 Hz.

Is there a more accurate way to calculate these derivatives on a sample by
sample basis in the way I've already used?  
(I need to use the results in real-time as the signals are coming in)

Thanks.
J

On Dec 10, 6:58 am, "Jam" <JamKrom...@gmail.com> wrote:
> First some info about what I'm trying to do:
>
> Basically I am working with EEG brain signals which are being digitized
> through a DAQ device (at around 6Khz) and output in a real-time audio
> program (Max/MSP at 48 or 44.1Khz).
>
> I am trying to code an external object in C++ that will compute time
> domain slope descriptor parameters as defined by Horth [1970]
>
> The equations for the parameters are given on page 3 herehttp://www.eurasip.org/Proceedings/Eusipco/Eusipco2007/Papers/C1L-E03...
>
> Basically:
> 1) Activity = variance of input signal x
> 2) Mobility = sqrt[var(x')/var(x)]
> 3) Complexity = sqrt[Mobility(x')/Mobility(x)]
>
> I have calculated the variance easily enough but the derivatives are
> giving me problems.
>
> At the moment I have implemented the central difference method using the
> following approach in the signal loop:
>
>                 X = InputSample;
>                 1stDeriv = (X - PrevPrevX)/(2*(1/SampleRate));
>
>                 PrevPrevX = PrevX;
>                 PrevX = X;
>
>                 OutputSample = 1stDeriv;
>
> This works well for the 1st derivative (a 10Hz Sin passed through gives a
> phase shifted version with amplitude of roughly 10*2*Pi).
>
> However when I pass the function through a second time to get the 2nd
> derivative the amplitude ends up being huge! (around 18000!)
>
> Also I've read that high frequency noise can result from this type of
> numerical differentiation but his shouldn't be too much of an issue as the
> only frequencies of interest to me are < 100 Hz.
>
> Is there a more accurate way to calculate these derivatives on a sample by
> sample basis in the way I've already used?
> (I need to use the results in real-time as the signals are coming in)
>
> Thanks.
> J

These are sometimes usefull.

http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_smoothing_filter

On Dec 10, 7:11 am, stanp <stan.p...@gmail.com> wrote:
> On Dec 10, 6:58 am, "Jam" <JamKrom...@gmail.com> wrote:
>
>
>
> > First some info about what I'm trying to do:
>
> > Basically I am working with EEG brain signals which are being digitized
> > through a DAQ device (at around 6Khz) and output in a real-time audio
> > program (Max/MSP at 48 or 44.1Khz).
>
> > I am trying to code an external object in C++ that will compute time
> > domain slope descriptor parameters as defined by Horth [1970]
>
> > The equations for the parameters are given on page 3 herehttp://www.eurasip.org/Proceedings/Eusipco/Eusipco2007/Papers/C1L-E03...
>
> > Basically:
> > 1) Activity = variance of input signal x
> > 2) Mobility = sqrt[var(x')/var(x)]
> > 3) Complexity = sqrt[Mobility(x')/Mobility(x)]
>
> > I have calculated the variance easily enough but the derivatives are
> > giving me problems.
>
> > At the moment I have implemented the central difference method using the
> > following approach in the signal loop:
>
> >                 X = InputSample;
> >                 1stDeriv = (X - PrevPrevX)/(2*(1/SampleRate));
>
> >                 PrevPrevX = PrevX;
> >                 PrevX = X;
>
> >                 OutputSample = 1stDeriv;
>
> > This works well for the 1st derivative (a 10Hz Sin passed through gives a
> > phase shifted version with amplitude of roughly 10*2*Pi).
>
> > However when I pass the function through a second time to get the 2nd
> > derivative the amplitude ends up being huge! (around 18000!)
>
> > Also I've read that high frequency noise can result from this type of
> > numerical differentiation but his shouldn't be too much of an issue as the
> > only frequencies of interest to me are < 100 Hz.
>
> > Is there a more accurate way to calculate these derivatives on a sample by
> > sample basis in the way I've already used?
> > (I need to use the results in real-time as the signals are coming in)
>
> > Thanks.
> > J
>
> These are sometimes usefull.
>
> http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_smoothing_filter
We use the Savitsky Golay algorithm to smooth graphical data to show
velocities.  The Wikipedia article is great for showing what the
filter does but in Numerical Recipes in C the the code is provided
that allows one to choose the order of the filter, the order of the
derivative and where the "center point" is.   We have found that using
the Savitsky Golay algorithm as a filter doesn't work near as well as
it does as an after the fact smoother.  If one tries to use the SG
algorithm as a filter with 7 or 8 proceeding data points it tends to
over exaggerate the velocity or acceleration at the latest point.
Try it.  It doesn't take long to implement but I don't think you will
like the results.  If you make the filter order too low the SG filter
will round the pulses.  If you make the filter order too high the
filter will exaggerate the current point.

It is very hard to calculate accelerations in signals that have pulses
and no models.  There needs to be some way of predicting and expected
value and updating that with the real value.  When we implemented our
graphs we did not have all the modeling techniques we do now.  When we
get a chance I would like to replace the SG algorithm with an observer
based filter.  I know the observer based filter would be far superior
but then I have a model.

Peter Nachtwey

On Dec 10, 6:58 am, "Jam" <JamKrom...@gmail.com> wrote:
> First some info about what I'm trying to do:
>
> Basically I am working with EEG brain signals which are being digitized
> through a DAQ device (at around 6Khz) and output in a real-time audio
> program (Max/MSP at 48 or 44.1Khz).
>
> I am trying to code an external object in C++ that will compute time
> domain slope descriptor parameters as defined by Horth [1970]
>
> The equations for the parameters are given on page 3 herehttp://www.eurasip.org/Proceedings/Eusipco/Eusipco2007/Papers/C1L-E03...
>
> Basically:
> 1) Activity = variance of input signal x
> 2) Mobility = sqrt[var(x')/var(x)]
> 3) Complexity = sqrt[Mobility(x')/Mobility(x)]
>
> I have calculated the variance easily enough but the derivatives are
> giving me problems.
>
> At the moment I have implemented the central difference method using the
> following approach in the signal loop:
>
>                 X = InputSample;
>                 1stDeriv = (X - PrevPrevX)/(2*(1/SampleRate));
>
>                 PrevPrevX = PrevX;
>                 PrevX = X;
>
>                 OutputSample = 1stDeriv;
>
> This works well for the 1st derivative (a 10Hz Sin passed through gives a
> phase shifted version with amplitude of roughly 10*2*Pi).
>
> However when I pass the function through a second time to get the 2nd
> derivative the amplitude ends up being huge! (around 18000!)
>
> Also I've read that high frequency noise can result from this type of
> numerical differentiation but his shouldn't be too much of an issue as the
> only frequencies of interest to me are < 100 Hz.
>
> Is there a more accurate way to calculate these derivatives on a sample by
> sample basis in the way I've already used?
> (I need to use the results in real-time as the signals are coming in)
>
> Thanks.
> J

You could use a more complicated FIR differentiator. Since your
frequencies of interest are so much smaller than your sample rate, you
probably wouldn't need a very long one. As always, you have a tradeoff
between filter length and how accurate the differentiation operation
is over the band of interest. If you do this, make sure you take the
delay of the filter into account, though, if you're going to make
comparisons between the signal and its derivatives at what should be
the same time instant.

Jason

On Dec 11, 5:40 pm, cincy...@gmail.com wrote:
> On Dec 10, 6:58 am, "Jam" <JamKrom...@gmail.com> wrote:
>
>
>
> > First some info about what I'm trying to do:
>
> > Basically I am working with EEG brain signals which are being digitized
> > through a DAQ device (at around 6Khz) and output in a real-time audio
> > program (Max/MSP at 48 or 44.1Khz).
>
> > I am trying to code an external object in C++ that will compute time
> > domain slope descriptor parameters as defined by Horth [1970]
>
> > The equations for the parameters are given on page 3 herehttp://www.eurasip.org/Proceedings/Eusipco/Eusipco2007/Papers/C1L-E03...
>
> > Basically:
> > 1) Activity = variance of input signal x
> > 2) Mobility = sqrt[var(x')/var(x)]
> > 3) Complexity = sqrt[Mobility(x')/Mobility(x)]
>
> > I have calculated the variance easily enough but the derivatives are
> > giving me problems.
>
> > At the moment I have implemented the central difference method using the
> > following approach in the signal loop:
>
> >                 X = InputSample;
> >                 1stDeriv = (X - PrevPrevX)/(2*(1/SampleRate));
>
> >                 PrevPrevX = PrevX;
> >                 PrevX = X;
>
> >                 OutputSample = 1stDeriv;
>
> > This works well for the 1st derivative (a 10Hz Sin passed through gives a
> > phase shifted version with amplitude of roughly 10*2*Pi).
>
> > However when I pass the function through a second time to get the 2nd
> > derivative the amplitude ends up being huge! (around 18000!)
>
> > Also I've read that high frequency noise can result from this type of
> > numerical differentiation but his shouldn't be too much of an issue as the
> > only frequencies of interest to me are < 100 Hz.
>
> > Is there a more accurate way to calculate these derivatives on a sample by
> > sample basis in the way I've already used?
> > (I need to use the results in real-time as the signals are coming in)
>
> > Thanks.
> > J
>
> You could use a more complicated FIR differentiator. Since your
> frequencies of interest are so much smaller than your sample rate, you
> probably wouldn't need a very long one. As always, you have a tradeoff
> between filter length and how accurate the differentiation operation
> is over the band of interest. If you do this, make sure you take the
> delay of the filter into account, though, if you're going to make
> comparisons between the signal and its derivatives at what should be
> the same time instant.
>
> Jason
The SG algorithm basically provides the coefficients for a FIR filter
that provides a least squares fit for the given order and number of
points you choose.

Peter Nachtwey

On Mon, 10 Dec 2007 05:58:12 -0600, "Jam" <JamKromium@gmail.com>
wrote:

>First some info about what I'm trying to do:
>
   (snipped)
>
>I have calculated the variance easily enough but the derivatives are
>giving me problems.
>
>At the moment I have implemented the central difference method using the
>following approach in the signal loop:
>
>		X = InputSample;
>		1stDeriv = (X - PrevPrevX)/(2*(1/SampleRate));
>			
>		PrevPrevX = PrevX;
>		PrevX = X;
>
>		OutputSample = 1stDeriv;
>
>This works well for the 1st derivative (a 10Hz Sin passed through gives a
>phase shifted version with amplitude of roughly 10*2*Pi).
>
>However when I pass the function through a second time to get the 2nd
>derivative the amplitude ends up being huge! (around 18000!)
>
>Also I've read that high frequency noise can result from this type of
>numerical differentiation but his shouldn't be too much of an issue as the
>only frequencies of interest to me are < 100 Hz.
>
>Is there a more accurate way to calculate these derivatives on a sample by
>sample basis in the way I've already used?  
>(I need to use the results in real-time as the signals are coming in)
>
>Thanks.
>J

Hello J,
  Jason's post seems very sensible to me.
However, a "central difference" differentiator 
does NOT amplify high frequency noise nearly 
as much as a simple "first difference" differentiator.
In fact, the "central difference" differentiator 
has infinite attenuation as Fs/2 (half the sample 
rate).

I don't know why your 2nd derivative is giving such 
large amplitudes, but in theory the 2nd derivative 
amplitudes should merely be the square of your 
1st derivative amplitudes. 

Good Luck,
[-Rick-]

On Wed, 12 Dec 2007 14:03:24 GMT, R.Lyons@_BOGUS_ieee.org (Rick Lyons)
wrote:


J,
 I forgot to mention, you might find
the following article to be interesting.

   http://electronicdesign.com/Articles/ArticleID/13358/13358.html

[-Rick-]
 

Thanks for the replies, I will start looking into the information
provided.

I have since modeled the three equations in my original post in Simulink
(using central difference) and they work as desired.  

I am limited to float type (no double unfortunately) in the C++ system and
using the derivatives as initially used, the 3 equations give constant
NaNs.    

I don't have much experience with FIR filters but will look into using
them here.

Rick, I read your article "A Differentiator With a Difference" and was
wondering why you haven't divided the difference by (1/sample rate) as in
most numerical methods?  
I tried implementing your equations in my code but they also give very
large values for the 2nd derivative.

J

On Wed, 12 Dec 2007 10:46:06 -0600, "Jam" <JamKromium@gmail.com>
wrote:

>Thanks for the replies, I will start looking into the information
>provided.
>
>I have since modeled the three equations in my original post in Simulink
>(using central difference) and they work as desired.  
>
>I am limited to float type (no double unfortunately) in the C++ system and
>using the derivatives as initially used, the 3 equations give constant
>NaNs.    
>
>I don't have much experience with FIR filters but will look into using
>them here.
>
>Rick, I read your article "A Differentiator With a Difference" and was
>wondering why you haven't divided the difference by (1/sample rate) as in
>most numerical methods?  
>I tried implementing your equations in my code but they also give very
>large values for the 2nd derivative.
>
>J

Hi Jam,
  well, ... I didn't worry about any "scaling" of the 
differentiator's output by a constant.  I only worried 
about (1) minimizing the arithmetic computations, 
and (2) over what frequency range was the differentiator's 
frequency response was reasonably linear.

See Ya',
[-Rick-]

On Wed, 12 Dec 2007 10:46:06 -0600, "Jam" <JamKromium@gmail.com>
wrote:

>Thanks for the replies, I will start looking into the information
>provided.
>
>I have since modeled the three equations in my original post in Simulink
>(using central difference) and they work as desired.  
>
>I am limited to float type (no double unfortunately) in the C++ system and
>using the derivatives as initially used, the 3 equations give constant
>NaNs.    
>
>I don't have much experience with FIR filters but will look into using
>them here.
>
>Rick, I read your article "A Differentiator With a Difference" and was
>wondering why you haven't divided the difference by (1/sample rate) as in
>most numerical methods?  
>I tried implementing your equations in my code but they also give very
>large values for the 2nd derivative.
>
>J

Whoops.  I forgot to say, in my previous post, 
that the freq magnitude plots in that "A Differentiator 
With a Difference" blog were "normalized" so that 
we could compare the freq ranges of linear operation 
of the various differentiators.  In reality, that 
"A Differentiator With a Difference" has a magnitude 
gain of roughly 1.05 at a freq of Fs/10.  At that same 
freq two of those differentiators cascaded should 
have a magnitude gain of roughly 1.12.

(At a freq of Fs/4, two of those differentiators cascaded 
should have a magnitude gain of roughly 4.5.)

So I'm not sure why you are obtaining very
large values for the 2nd derivative.

See Ya,
[-Rick-]