Instantaneous Frequency with I-Q demodulation

Hello everybody,

I have a few questions, maybe you can help me. I have 2 signal I,Q @ 40khz
and the task is to obtain the instantaneous frequency. What i did in matlab
is this: first of all i found the phase with atan2, i unwrapped the phase
and I derivated in order to have the frequency (with command diff).

A first issue is here: I had to average the phase to see somenthings in
frequency 'cause otherwise seems very noisy the istantaneous freq ( a lot a
spurious peaks). Is there a method to optimize this average or to avoid?
(now i use a very long moving average filter)

Second issue: i would like to implement it on a 16 bit microcontroller. in
order to avoid phase unwrapping i read is possibile to use only atan as
indicated here:

http://www.dsprelated.com/showmessage/19121/2.php

but I don't understand why. If anyone could explain me... would be
awesome.

Thanks in advance 

Marco  

On Mar 5, 9:06&#4294967295;am, "freedreamer" <freedrea...@email.it> wrote:
> Hello everybody,
>
> I have a few questions, maybe you can help me. I have 2 signal I,Q @ 40khz
> and the task is to obtain the instantaneous frequency. What i did in matlab
> is this: first of all i found the phase with atan2, i unwrapped the phase
> and I derivated in order to have the frequency (with command diff).
>
> A first issue is here: I had to average the phase to see somenthings in
> frequency 'cause otherwise seems very noisy the istantaneous freq ( a lot a
> spurious peaks). Is there a method to optimize this average or to avoid?
> (now i use a very long moving average filter)
>
> Second issue: i would like to implement it on a 16 bit microcontroller. in
> order to avoid phase unwrapping i read is possibile to use only atan as
> indicated here:
>
> http://www.dsprelated.com/showmessage/19121/2.php
>
> but I don't understand why. If anyone could explain me... would be
> awesome.
>
> Thanks in advance
>
> Marco &#4294967295;

Did you try this?

f = angle(x(2:N).*conj(x(1:N-1)))/2/pi;

John

>On Mar 5, 9:06=A0am, "freedreamer" <freedrea...@email.it> wrote:
>> Hello everybody,
>>
>> I have a few questions, maybe you can help me. I have 2 signal I,Q @
40kh=
>z
>> and the task is to obtain the instantaneous frequency. What i did in
matl=
>ab
>> is this: first of all i found the phase with atan2, i unwrapped the
phase
>> and I derivated in order to have the frequency (with command diff).
>>
>> A first issue is here: I had to average the phase to see somenthings
in
>> frequency 'cause otherwise seems very noisy the istantaneous freq ( a
lot=
> a
>> spurious peaks). Is there a method to optimize this average or to
avoid?
>> (now i use a very long moving average filter)
>>
>> Second issue: i would like to implement it on a 16 bit microcontroller.
i=
>n
>> order to avoid phase unwrapping i read is possibile to use only atan
as
>> indicated here:
>>
>> http://www.dsprelated.com/showmessage/19121/2.php
>>
>> but I don't understand why. If anyone could explain me... would be
>> awesome.
>>
>> Thanks in advance
>>
>> Marco =A0
>
>Did you try this?
>
>f =3D angle(x(2:N).*conj(x(1:N-1)))/2/pi;
>
>John
>

mmhh what should be x? in that way you obtain same result as atan i guess
, doesn't it?
maybe i didn't understand, Could you explain your mental construction?

thanks 

Marco

On Mar 5, 10:05&#4294967295;am, "freedreamer" <freedrea...@email.it> wrote:
> >On Mar 5, 9:06=A0am, "freedreamer" <freedrea...@email.it> wrote:
> >> Hello everybody,
>
> >> I have a few questions, maybe you can help me. I have 2 signal I,Q @
> 40kh=
> >z
> >> and the task is to obtain the instantaneous frequency. What i did in
> matl=
> >ab
> >> is this: first of all i found the phase with atan2, i unwrapped the
> phase
> >> and I derivated in order to have the frequency (with command diff).
>
> >> A first issue is here: I had to average the phase to see somenthings
> in
> >> frequency 'cause otherwise seems very noisy the istantaneous freq ( a
> lot=
> > a
> >> spurious peaks). Is there a method to optimize this average or to
> avoid?
> >> (now i use a very long moving average filter)
>
> >> Second issue: i would like to implement it on a 16 bit microcontroller.
> i=
> >n
> >> order to avoid phase unwrapping i read is possibile to use only atan
> as
> >> indicated here:
>
> >>http://www.dsprelated.com/showmessage/19121/2.php
>
> >> but I don't understand why. If anyone could explain me... would be
> >> awesome.
>
> >> Thanks in advance
>
> >> Marco =A0
>
> >Did you try this?
>
> >f =3D angle(x(2:N).*conj(x(1:N-1)))/2/pi;
>
> >John
>
> mmhh what should be x? in that way you obtain same result as atan i guess
> , doesn't it?
> maybe i didn't understand, Could you explain your mental construction?
>
> thanks
>
> Marco- Hide quoted text -
>
> - Show quoted text -

x is the complex-valued input signal that you want the IF of.
N is length(x).
angle is atan2.
pi is the ratio of a circle's circumference to its diameter.

John

It sounds like you want to do FM discrimination.  Try:

>> diff(unwrap(angle(x)));

Boo on you if this is a homework problem!

Darrell

On Mar 5, 10:36&#4294967295;am, Darrell <darrel...@gmail.com> wrote:
> It sounds like you want to do FM discrimination. &#4294967295;Try:
>
> >> diff(unwrap(angle(x)));
>
> Boo on you if this is a homework problem!
>
> Darrell

Boo on anyone who uses diff() to take a derivative without considering
how good an approximation it is to an actual derivative function.  If
you are not highly oversampled, and your input signal is not
bandlimited you got problems.

On Mar 6, 3:06&#4294967295;am, "freedreamer" <freedrea...@email.it> wrote:
> Hello everybody,
>
> I have a few questions, maybe you can help me. I have 2 signal I,Q @ 40khz
> and the task is to obtain the instantaneous frequency. What i did in matlab
> is this: first of all i found the phase with atan2, i unwrapped the phase
> and I derivated in order to have the frequency (with command diff).
>
> A first issue is here: I had to average the phase to see somenthings in
> frequency 'cause otherwise seems very noisy the istantaneous freq ( a lot a
> spurious peaks). Is there a method to optimize this average or to avoid?
> (now i use a very long moving average filter)
>
> Second issue: i would like to implement it on a 16 bit microcontroller. in
> order to avoid phase unwrapping i read is possibile to use only atan as
> indicated here:
>
> http://www.dsprelated.com/showmessage/19121/2.php
>
> but I don't understand why. If anyone could explain me... would be
> awesome.
>
> Thanks in advance
>
> Marco &#4294967295;

You don't need to use atan2 numerically to find this rate of change of
phase. Instead differntiate analytically atan Q/I. Then use this
result instead.


Hardy

On Mar 5, 9:06&#4294967295;am, "freedreamer" <freedrea...@email.it> wrote:
> Hello everybody,
>
> I have a few questions, maybe you can help me. I have 2 signal I,Q @ 40khz
> and the task is to obtain the instantaneous frequency. What i did in matlab
> is this: first of all i found the phase with atan2, i unwrapped the phase
> and I derivated in order to have the frequency (with command diff).
>
> A first issue is here: I had to average the phase to see somenthings in
> frequency 'cause otherwise seems very noisy the istantaneous freq ( a lot a
> spurious peaks). Is there a method to optimize this average or to avoid?
> (now i use a very long moving average filter)
>
> Second issue: i would like to implement it on a 16 bit microcontroller. in
> order to avoid phase unwrapping i read is possibile to use only atan as
> indicated here:
>
> http://www.dsprelated.com/showmessage/19121/2.php
>
> but I don't understand why. If anyone could explain me... would be
> awesome.

rather than compute the angle for every sample (which would be between
-pi and +pi), unwrapping it (which is an icky "if-then-else" operation
which some DSP chips are not so good with), and then diff-ing it, a
better method is to compute from the two complex numbers what the
angle difference is, which should be small.

let
     v[n] = x[n] + j*y[n] = |v[n]|*exp(j*arg{v[n]})

     arg{v[n]} - arg{v[n-1]} = arg{ v[n]/v[n-1] }

          = arg{ (x[n]+j*y[n])/(x[n-1]+j*y[n-1]) }

          = arg{ ((x[n]+j*y[n])*(x[n-1]-j*y[n-1]))/((x[n-1]+j*y[n-1])*
(x[n-1]-j*y[n-1])) }

          = arg{ (x[n]*x[n-1]+y[n]*y[n-1]) + j*(y[n]*x[n-1]-x[n]*y
[n-1]) }

now, if we can assume that the angular difference is small between
adjacent samples, we can assume we are always in the right half-plane
and apply directly the arctan() function.

     arg{v[n]} - arg{v[n-1]} = arctan( (y[n]*x[n-1]-x[n]*y[n-1])/(x[n]
*x[n-1]+y[n]*y[n-1]) )

note the division in the above.

here is a pretty good approximation to arctan that needs no additional
division:

    arctan(u) ~= u * p(u^2)       for -1 <= u <= 1

where

    p(x) =        1.0
             +    -3.31327349395871e-01  x
             +     1.82590340927122e-01  x^2
             +    -8.73985945234805e-02  x^3
             +     2.15337663896783e-02  x^4

i can't remember what the error is, but it's small.

if the magnitude of v[n] is roughly constant, and you absolutely hate
doing division, this might be useful:

     1/(1+x) ~=   1.0
             +    -9.99351503499687e-01    *  x
             +     9.86903179099794e-01    *  x^2
             +    -9.07341721411229e-01    *  x^3
             +     6.73561921455864e-01    *  x^4
             +    -3.26278125828937e-01    *  x^5
             +     7.25062501841954e-02    *  x^6

relative error no greater than  1.05725067401401e-05  for 0 <= x <= 1.


> Thanks in advance

FWIW.

r b-j

thanks, now is much clear. :D 

>On Mar 5, 9:06=A0am, "freedreamer" <freedrea...@email.it> wrote:
>> Hello everybody,
>>
>> I have a few questions, maybe you can help me. I have 2 signal I,Q @
40kh=
>z
>> and the task is to obtain the instantaneous frequency. What i did in
matl=
>ab
>> is this: first of all i found the phase with atan2, i unwrapped the
phase
>> and I derivated in order to have the frequency (with command diff).
>>
>> A first issue is here: I had to average the phase to see somenthings
in
>> frequency 'cause otherwise seems very noisy the istantaneous freq ( a
lot=
> a
>> spurious peaks). Is there a method to optimize this average or to
avoid?
>> (now i use a very long moving average filter)
>>
>> Second issue: i would like to implement it on a 16 bit microcontroller.
i=
>n
>> order to avoid phase unwrapping i read is possibile to use only atan
as
>> indicated here:
>>
>> http://www.dsprelated.com/showmessage/19121/2.php
>>
>> but I don't understand why. If anyone could explain me... would be
>> awesome.
>
>rather than compute the angle for every sample (which would be between
>-pi and +pi), unwrapping it (which is an icky "if-then-else" operation
>which some DSP chips are not so good with), and then diff-ing it, a
>better method is to compute from the two complex numbers what the
>angle difference is, which should be small.
>
>let
>     v[n] =3D x[n] + j*y[n] =3D |v[n]|*exp(j*arg{v[n]})
>
>     arg{v[n]} - arg{v[n-1]} =3D arg{ v[n]/v[n-1] }
>
>          =3D arg{ (x[n]+j*y[n])/(x[n-1]+j*y[n-1]) }
>
>          =3D arg{ ((x[n]+j*y[n])*(x[n-1]-j*y[n-1]))/((x[n-1]+j*y[n-1])*
>(x[n-1]-j*y[n-1])) }
>
>          =3D arg{ (x[n]*x[n-1]+y[n]*y[n-1]) + j*(y[n]*x[n-1]-x[n]*y
>[n-1]) }
>
>now, if we can assume that the angular difference is small between
>adjacent samples, we can assume we are always in the right half-plane
>and apply directly the arctan() function.
>
>     arg{v[n]} - arg{v[n-1]} =3D arctan( (y[n]*x[n-1]-x[n]*y[n-1])/(x[n]
>*x[n-1]+y[n]*y[n-1]) )
>
>note the division in the above.
>
>here is a pretty good approximation to arctan that needs no additional
>division:
>
>    arctan(u) ~=3D u * p(u^2)       for -1 <=3D u <=3D 1
>
>where
>
>    p(x) =3D        1.0
>             +    -3.31327349395871e-01  x
>             +     1.82590340927122e-01  x^2
>             +    -8.73985945234805e-02  x^3
>             +     2.15337663896783e-02  x^4
>
>i can't remember what the error is, but it's small.
>
>if the magnitude of v[n] is roughly constant, and you absolutely hate
>doing division, this might be useful:
>
>     1/(1+x) ~=3D   1.0
>             +    -9.99351503499687e-01    *  x
>             +     9.86903179099794e-01    *  x^2
>             +    -9.07341721411229e-01    *  x^3
>             +     6.73561921455864e-01    *  x^4
>             +    -3.26278125828937e-01    *  x^5
>             +     7.25062501841954e-02    *  x^6
>
>relative error no greater than  1.05725067401401e-05  for 0 <=3D x <=3D
1.
>
>
>> Thanks in advance
>
>FWIW.
>
>r b-j
>

On Mar 12, 3:47&#4294967295;pm, "freedreamer" <freedrea...@email.it> wrote:
> thanks, now is much clear. :D
>
>

yer welcome.

in case you haven't, you need to unwrap some lines of ASCII math
because Google Groups wrapped it at unexpected places.

the power series for arctan(x) or 1/(1+x) might be more precision than
you need.  may or may not be useful.

r b-j