numerical differentiation of oscillating singal

Hi all,
  I am looking for a numerical method for differentiation of
oscillating signal. I have tried several method like forward/central
difference method, fourier differentiator but no help. Does anyone
have any idea?

"Rex_chaos" <rex_chaos@21cn.com> wrote in message
news:f7a7417.0310280538.3d3322dc@posting.google.com...
> Hi all,
>   I am looking for a numerical method for differentiation of
> oscillating signal. I have tried several method like forward/central
> difference method, fourier differentiator but no help. Does anyone
> have any idea?

If the signal is properly sampled, then it respresents a continuous
bandlimited signal.  If you want the derivative of this signal, there are
nifty ways to get it, depending on exactly how, for the purpose of
extracting its derivative, you want that signal to be reconstructed.

Do you want the derivative of the ideally interpolated signal, the signal
that shows up at the output of some system, or the signal as interpolated
with some other interpolation filter?

Use

  B=REMEZ(N,F,A,'differentiator') and B=REMEZ(N,F,A,W,'differentiator')

For bandlimited differentiation. The problem with differentiation is the
noise at high frequencies.

HTH
-- 
-----------------------------------------------------------
Dr. Izhak

"Rex_chaos" <rex_chaos@21cn.com> wrote in message
news:f7a7417.0310280538.3d3322dc@posting.google.com...
> Hi all,
>   I am looking for a numerical method for differentiation of
> oscillating signal. I have tried several method like forward/central
> difference method, fourier differentiator but no help. Does anyone
> have any idea?

"Matt Timmermans" <mt0000@sympatico.nospam-remove.ca> wrote in message news:<R8vnb.1971$Nz5.205470@news20.bellglobal.com>...
> "Rex_chaos" <rex_chaos@21cn.com> wrote in message
> news:f7a7417.0310280538.3d3322dc@posting.google.com...
> > Hi all,
> >   I am looking for a numerical method for differentiation of
> > oscillating signal. I have tried several method like forward/central
> > difference method, fourier differentiator but no help. Does anyone
> > have any idea?
> 
> If the signal is properly sampled, then it respresents a continuous
> bandlimited signal.  If you want the derivative of this signal, there are
> nifty ways to get it, depending on exactly how, for the purpose of
> extracting its derivative, you want that signal to be reconstructed.
> 
> Do you want the derivative of the ideally interpolated signal, the signal
> that shows up at the output of some system, or the signal as interpolated
> with some other interpolation filter?

> 
Actually, I am processing a signal f(t) and in each step of time, 

 f(t+dt) = G[D[f(t)]]

here D means the differentiation of f(t) and G is a function to form
D[f(t)] to a new signal. However, f(t) is a signal with a rapidly
oscillating tail. A common algorithm (such as forward difference) to
perform numerical differentiation will lead to wrong result. For this
purpose, I adopt a high-order difference formula. It's weird that the
high-order formula is even worse than the central difference formula.
I also take the fourier method but no help

D[f(t)] = iff(i*w*F(w));

w denotes for the angular frequency and ifft denotes the inverse
fourier transformation. F(w) is the fourier spectrum of f(t).

>   B=REMEZ(N,F,A,'differentiator') and B=REMEZ(N,F,A,W,'differentiator')
> 
> For bandlimited differentiation. The problem with differentiation is the
> noise at high frequencies.
What is remez? Is it a function provided by matlab or other dsp package?

Thanks.

Rex_chaos wrote:
>>  B=REMEZ(N,F,A,'differentiator') and B=REMEZ(N,F,A,W,'differentiator')
>>
>>For bandlimited differentiation. The problem with differentiation is the
>>noise at high frequencies.
> 
> What is remez? Is it a function provided by matlab or other dsp package?
> 

from matlab:

  REMEZ Parks-McClellan optimal equiripple FIR filter design.
     B=REMEZ(N,F,A) returns a length N+1 linear phase (real, symmetric
     coefficients) FIR filter which has the best approximation to the
     desired frequency response described by F and A in the minimax
     sense. F is a vector of frequency band edges in pairs, in ascending
     order between 0 and 1. 1 corresponds to the Nyquist frequency or half
     the sampling frequency. A is a real vector the same size as F
     which specifies the desired amplitude of the frequency response of the
     resultant filter B. The desired response is the line connecting the
     points (F(k),A(k)) and (F(k+1),A(k+1)) for odd k; REMEZ treats the
     bands between F(k+1) and F(k+2) for odd k as "transition bands" or
     "don't care" regions. Thus the desired amplitude is piecewise linear
     with transition bands. The maximum error is minimized.

     For filters with a gain other than zero at Fs/2, e.g., highpass
     and bandstop filters, N must be even.  Otherwise, N will be
     incremented by one.

     B=REMEZ(N,F,A,W) uses the weights in W to weight the error. W has one
     entry per band (so it is half the length of F and A) which tells
     REMEZ how much emphasis to put on minimizing the error in each band
     relative to the other bands.

     B=REMEZ(N,F,A,'Hilbert') and B=REMEZ(N,F,A,W,'Hilbert') design filters
     that have odd symmetry, that is, B(k) = -B(N+2-k) for k = 1, ..., N+1.
     A special case is a Hilbert transformer which has an approx. amplitude
     of 1 across the entire band, e.g. B=REMEZ(30,[.1 .9],[1 1],'Hilbert').

     B=REMEZ(N,F,A,'differentiator') and B=REMEZ(N,F,A,W,'differentiator')
     also design filters with odd symmetry, but with a special weighting
     scheme for non-zero amplitude bands. The weight is assumed to be equal
     to the inverse of frequency times the weight W. Thus the filter has a
     much better fit at low frequency than at high frequency. This designs
     FIR differentiators.

     B=REMEZ(...,{LGRID}), where {LGRID} is a one-by-one cell array 
containing
     an integer, controls the density of the frequency grid. The 
frequency grid
     size is roughly LGRID*N/2*BW, where BW is the fraction of the total 
band
     interval [0,1] covered by F. LGRID should be no less than its 
default of 16.
     Increasing LGRID often results in filters which are more exactly equi-
     ripple, at the expense of taking longer to compute.

     [B,ERR]=REMEZ(...) returns the maximum ripple height ERR.

     [B,ERR,RES]=REMEZ(...) returns a structure RES of optional results
     computed by REMEZ, and contains the following fields:

        RES.fgrid: vector containing the frequency grid used in
                   the filter design optimization
          RES.des: desired response on fgrid
           RES.wt: weights on fgrid
            RES.H: actual frequency response on the grid
        RES.error: error at each point on the frequency grid (desired - 
actual)
        RES.iextr: vector of indices into fgrid of extremal frequencies
        RES.fextr: vector of extremal frequencies

     See also CREMEZ, FIRLS, FIR1, FIR2, BUTTER, CHEBY1, CHEBY2, ELLIP,
     FREQZ, FILTER and GREMEZ in the Filter Design toolbox.

Thanks. I try that function. Maybe I misuse it, the effect of the
function it not like a differentiator. Would you please give me an
example to illustrate how to differentiate a signal numerical with a
filter?

Thanks a lot

Rex_chaos wrote:

> Thanks. I try that function. Maybe I misuse it, the effect of the
> function it not like a differentiator. Would you please give me an
> example to illustrate how to differentiate a signal numerical with a
> filter?
> 
> Thanks a lot

Here's a silly-simple example:

Suppose the differentiator function is the simple difference between the
current sample and the previous one. As a difference function,
x[i] - x[i-1]. This is also a transversal FIR whose tap values are
1, -1. (Latest value first.) Any function whose value depends only on
past inputs can be represented as a filter.

Jerry

P.S. That structure is also a silly-simple Hilbert transformer, giving a
90-degree phase shift as a differentiator must. It's main trouble is
that its output is half an interval delayed with respect to the
undifferentiated signal. We can align signal and derivative by delaying
the signal. Spanning three samples in the FIFO is simpler than creating
a half-period delay. The differentiator (and Hilbert) output is computed
with the taps 1, 0, -1, while the merely delayed output uses 0, 1, 0.

Math fills in details that intuition doesn't reveal, but using intuition
at first helps us to see where the math is leading.
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

"Rex_chaos" <rex_chaos@21cn.com> wrote in message
news:f7a7417.0310282119.35acf7bd@posting.google.com...
> Actually, I am processing a signal f(t) and in each step of time,
>
>  f(t+dt) = G[D[f(t)]]
>
> here D means the differentiation of f(t) and G is a function to form
> D[f(t)] to a new signal. [...]

Oh, you're doing a simulation?  In that case, you can't really get an
accurate derivative from amplitude samples, because the derivative depends
on past *and future* samples.  If you don't have the future samples, you can
only guess at the derivative, and and make your time step small enough to
provide sufficient accuracy.  If your function is oscillating quickly with
respect to the time step, then it won't be very accurate at all.

Now, your particular case looks very broken.  The process you wrote doesn't
seem to converge to a smooth function as you decrease the time step, which
would explain why you're getting such wierd results.  Do you really mean for
those f's to be the same signal?

> Oh, you're doing a simulation?  
yes.

>In that case, you can't really get an
> accurate derivative from amplitude samples, because the derivative depends
> on past *and future* samples.  If you don't have the future samples, you can
> only guess at the derivative, and and make your time step small enough to
> provide sufficient accuracy.  If your function is oscillating quickly with
> respect to the time step, then it won't be very accurate at all.
Of course, I did all you mentioned above.
 
> Now, your particular case looks very broken.  The process you wrote doesn't
> seem to converge to a smooth function as you decrease the time step, which
> would explain why you're getting such wierd results.  Do you really mean for
> those f's to be the same signal?
As long as the time step is small enough, the signal can be considered
as a smooth-vaired one, isn't it?