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
Instantaneous Frequency with I-Q demodulation
Started by ●March 5, 2009
Reply by ●March 5, 20092009-03-05
On Mar 5, 9:06�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 �Did you try this? f = angle(x(2:N).*conj(x(1:N-1)))/2/pi; John
Reply by ●March 5, 20092009-03-05
>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 inmatl=>ab >> is this: first of all i found the phase with atan2, i unwrapped thephase>> 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 somenthingsin>> frequency 'cause otherwise seems very noisy the istantaneous freq ( alot=> a >> spurious peaks). Is there a method to optimize this average or toavoid?>> (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 atanas>> 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
Reply by ●March 5, 20092009-03-05
On Mar 5, 10:05�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
Reply by ●March 5, 20092009-03-05
It sounds like you want to do FM discrimination. Try:>> diff(unwrap(angle(x)));Boo on you if this is a homework problem! Darrell
Reply by ●March 5, 20092009-03-05
On Mar 5, 10:36�am, Darrell <darrel...@gmail.com> wrote:> It sounds like you want to do FM discrimination. �Try: > > >> diff(unwrap(angle(x))); > > Boo on you if this is a homework problem! > > DarrellBoo 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.
Reply by ●March 5, 20092009-03-05
On Mar 6, 3:06�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 �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
Reply by ●March 5, 20092009-03-05
On Mar 5, 9:06�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 advanceFWIW. r b-j
Reply by ●March 12, 20092009-03-12
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 inmatl=>ab >> is this: first of all i found the phase with atan2, i unwrapped thephase>> 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 somenthingsin>> frequency 'cause otherwise seems very noisy the istantaneous freq ( alot=> a >> spurious peaks). Is there a method to optimize this average or toavoid?>> (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 atanas>> 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 <=3D1.> > >> Thanks in advance > >FWIW. > >r b-j >
Reply by ●March 13, 20092009-03-13
On Mar 12, 3:47�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