Resolving frequencies in an oversampled spectrum

Again, thank you all for the input that you have provided.  I've been
working with the ESPRIT algorithm, and it seems to be doing an
admirable job of separating the frequency components.  I just have a
few last second thoughts that I would appreciate a gut check on (if
you're not familiar with the term, it's a "yeah, that's about right" or
a "You're out of your mind!").
1.  I'm assuming in the esprit method one obtains the highest accuracy
if you choose the number of shifts so that the matrix which you first
perform SVD on is as nearly square as possible.
2.  I'm also assuming that the penalty associated is slower processing
of the SVD (I work in matlab, btw).
3.  I have not found any reference to a version of esprit which
generates amplitudes so I am doing that as a separate step after esprit
gives me the frequencies.
4.  I am also assuming that esprit cannot provide phase information,
likewise I'm extracting that from amplitude info (atan of overlap of
raw data with sin and cos of matched frequencies respectively.)  Does
anybody know of a better way to get the phase?  (You can even choose
whatever arbitrary reference you want.  See I can be generous!)

Again, thank you all so much for the assistance.

Dan

> I haven't been following, but you might try Peter K's frequency
> estimation page for resources:
>
>   http://home.comcast.net/~kootsoop/freqalgs.htm
>
> Best of luck!
>
> --Randy
>
> >
> >
> > Thanks,
> > Dan
> >
> >
> > cf wrote:
> >> Vladimir Vassilevsky wrote:
> >> >
> >> >
> >> > Dan K wrote:
> >> >
> >> >> I have a time series which is highly oversampled.  After low pass
> >> >> filtering, I have about 17,000 data points.  In that span I have a very
> >> >> low frequency signal (about 1/50th of a period), two more moderate
> >> >> frequency signals (can be between 1 and 8 periods) and the sum of those
> >> >> two signals.  From sweep to sweep the two signals will both move
> >> >> together in frequency space (by very small amounts) and their relative
> >> >> intensities can change dramatically.  What I need is a high accuracy
> >> >> method of determining the frequency of these two signals.
> >> >
> >> > If you can describe your signals as the models with the unknown
> >> > parameters, you can apply the Proni algorithm or some other method of
> >> > the parametric estimation.
> >> >
> >> > For example, if your signal is described by a model Y = A sin(Wt + Fi)
> >> > with A, W, Fi as the unknown parameters, you can solve for the values
> >> > which provide the best fit with your data.
> >> >
> >> > This is a very accurate method although it takes a lot of computation.
> >> >
> >> >
> >> > Vladimir Vassilevsky
> >> >
> >> > DSP and Mixed Signal Design Consultant
> >> >
> >> > http://www.abvolt.com
> >>
> >> Hi Dan,
> >>
> >> I agree with Vladimir.  If your signal model is sinusoidal, you can
> >> probably try using a high resolution parameteric estimation method like
> >> ESPRIT (Estimation of Signal Parameters via Rotational Invariance
> >> Technique) where the parameters you are trying to estimate are the
> >> different frequencies.  It basically involves partitioning your data
> >> into two overlap sets, e.g. x0 - x10000 and x9000 - x17000.  Then you
> >> will see that the second set of data is just a phase shifted version of
> >> the first set.  From this relationship, you can obtain a general eigen
> >> equation such that the frequency estimates are contained inside the
> >> eigenvalue matrix.  I found that this technique works reasonably well
> >> (measure by the variance of the estimates) starting at about SNR = 15 dB
> >> and improve as the SNR increases.  For detail, I suggest you look at
> >> some advanced DSP book like:
> >>
> >> Manolakis, Ingle, Kogon, "Statistical and adaptive signal processing :
> >> spectral estimation, signal modeling, adaptive filtering, and array
> >> processing", 2005.
> >>
> >> This is an e-book so you should be able to view it from an university's
> >> library website.
> >>
> >> Hope this helps,
> >>
> >> cf
> >
>
> --
> %  Randy Yates                  % "...the answer lies within your soul
> %% Fuquay-Varina, NC            %       'cause no one knows which side
> %%% 919-577-9882                %                   the coin will fall."
> %%%% <yates@ieee.org>           %  'Big Wheels', *Out of the Blue*, ELO
> http://home.earthlink.net/~yatescr

"Dan" <dkominsky@primephotonics.com> writes:

> Again, thank you all for the input that you have provided.  I've been
> working with the ESPRIT algorithm, and it seems to be doing an
> admirable job of separating the frequency components.  

Hi Dan,

I don't know what the ESPRIT algorithm does. Perhaps you could educate
me and possibly others in this group by giving us a quick summary 
and/or a few pointers/references. 
-- 
%  Randy Yates                  % "Ticket to the moon, flight leaves here today 
%% Fuquay-Varina, NC            %  from Satellite 2"
%%% 919-577-9882                % 'Ticket To The Moon' 
%%%% <yates@ieee.org>           % *Time*, Electric Light Orchestra
http://home.earthlink.net/~yatescr

Sorry Randy,
The best explanation I've found of ESPRIT so far can be found on the
following website:
http://cnx.org/content/m10588/latest/

If you use MatLab you can find many of the esprit modules associated
with what seems to be the standard book at:
http://www.ee.uwa.edu.au/~roberto/teach/as413/course/sasptoolbox/

Hope this helps you as much as it has me.

Dan

(BTW, much of this thread is also taking place in comp.dsp, which is
how I goofed on refering to something that was said in only one group.
oops.)


Randy Yates wrote:
> "Dan" <dkominsky@primephotonics.com> writes:
>
> > Again, thank you all for the input that you have provided.  I've been
> > working with the ESPRIT algorithm, and it seems to be doing an
> > admirable job of separating the frequency components.
>
> Hi Dan,
>
> I don't know what the ESPRIT algorithm does. Perhaps you could educate
> me and possibly others in this group by giving us a quick summary
> and/or a few pointers/references.
> --
> %  Randy Yates                  % "Ticket to the moon, flight leaves here today
> %% Fuquay-Varina, NC            %  from Satellite 2"
> %%% 919-577-9882                % 'Ticket To The Moon'
> %%%% <yates@ieee.org>           % *Time*, Electric Light Orchestra
> http://home.earthlink.net/~yatescr

"Dan K" <dkominsky@primephotonics.com> writes:

> Sorry Randy,
> The best explanation I've found of ESPRIT so far can be found on the
> following website:
> http://cnx.org/content/m10588/latest/

Written by our very own Julius Kusuma!

> If you use MatLab you can find many of the esprit modules associated
> with what seems to be the standard book at:
> http://www.ee.uwa.edu.au/~roberto/teach/as413/course/sasptoolbox/
>
> Hope this helps you as much as it has me.

Excellent! Thanks Dan!

--Randy

>
> Dan
>
> (BTW, much of this thread is also taking place in comp.dsp, which is
> how I goofed on refering to something that was said in only one group.
> oops.)
>
>
> Randy Yates wrote:
>> "Dan" <dkominsky@primephotonics.com> writes:
>>
>> > Again, thank you all for the input that you have provided.  I've been
>> > working with the ESPRIT algorithm, and it seems to be doing an
>> > admirable job of separating the frequency components.
>>
>> Hi Dan,
>>
>> I don't know what the ESPRIT algorithm does. Perhaps you could educate
>> me and possibly others in this group by giving us a quick summary
>> and/or a few pointers/references.
>> --
>> %  Randy Yates                  % "Ticket to the moon, flight leaves here today
>> %% Fuquay-Varina, NC            %  from Satellite 2"
>> %%% 919-577-9882                % 'Ticket To The Moon'
>> %%%% <yates@ieee.org>           % *Time*, Electric Light Orchestra
>> http://home.earthlink.net/~yatescr
>

-- 
%  Randy Yates                  % "Bird, on the wing,
%% Fuquay-Varina, NC            %   goes floating by
%%% 919-577-9882                %   but there's a teardrop in his eye..."
%%%% <yates@ieee.org>           % 'One Summer Dream', *Face The Music*, ELO
http://home.earthlink.net/~yatescr

Hi Dan,

1. I don't think making the augmented data matrix (the matrix you 
perform SVD on) as square as possible will make your estimates more 
accurate.  To make your estimates more accurate, 1) you should collect 
as much data as you can to construct augmented data matrix.  However, as 
I recall, you only have a fixed number of samples to work with, so this 
might not help you.  2) You can try "stacking" (stacking horizontally) 
more than two overlapping data sets (instead of just one (vertically) as 
described by many textbooks).  This should make your noise subspace 
bigger, thereby "absorbing" more of the noise from your signal (I don't 
have a good technically explanation for this, but I think this will 
work).  I have done this before for a related problem so you can try it 
here as well.  Pls. let me know how it goes.

2. Yeah, doing SVD incurs O(n^3) computational complexity (for a square 
n X n matrix), so it's pretty slow, esp. for big matrices.

3 and 4. I haven't seen any esprit based algorithm for amplitude and 
phase estimation either, but I have a feeling this can be done.  Perhaps 
I will look into it :) .  For amplitude estimation, I have seen a least 
squares based algorithm.  If you have access to IEEE papers, you can 
look at:

P. Stoica, H. Li and J. Li, "Amplitude Estimation of Sinusoidal Signals: 
Survey, New Results, and an Application", vol. 48(2), pp. 338-352, Feb. 
2000.

Hope this helps,

cf


Dan wrote:
> Again, thank you all for the input that you have provided.  I've been
> working with the ESPRIT algorithm, and it seems to be doing an
> admirable job of separating the frequency components.  I just have a
> few last second thoughts that I would appreciate a gut check on (if
> you're not familiar with the term, it's a "yeah, that's about right" or
> a "You're out of your mind!").
> 1.  I'm assuming in the esprit method one obtains the highest accuracy
> if you choose the number of shifts so that the matrix which you first
> perform SVD on is as nearly square as possible.
> 2.  I'm also assuming that the penalty associated is slower processing
> of the SVD (I work in matlab, btw).
> 3.  I have not found any reference to a version of esprit which
> generates amplitudes so I am doing that as a separate step after esprit
> gives me the frequencies.
> 4.  I am also assuming that esprit cannot provide phase information,
> likewise I'm extracting that from amplitude info (atan of overlap of
> raw data with sin and cos of matched frequencies respectively.)  Does
> anybody know of a better way to get the phase?  (You can even choose
> whatever arbitrary reference you want.  See I can be generous!)
> 
> Again, thank you all so much for the assistance.
> 
> Dan
> 

cf,

Thanks for the info.  I'm not entirely clear as to what you mean by
stacking overlapping data sets horizontally.  I believed that a
requirement of this method was that the data had to be contiguous so
that any adjacent columns of the augmented data matrix would be exactly
one step out of phase with each other.  Are you suggesting taking two
full data sets, creating the augmented matrix for each of them, and
then concatenating them horizontally to make a double size augmented
data matrix?  I would say that I clearly do not fully understand how
the math is working here yet.

I am taking a look at the paper that you referenced.  It may be a while
before I have the time to work all the way through it, but it does look
promising.  Thanks,

Dan

cf wrote:
> Hi Dan,
>
> 1. I don't think making the augmented data matrix (the matrix you
> perform SVD on) as square as possible will make your estimates more
> accurate.  To make your estimates more accurate, 1) you should collect
> as much data as you can to construct augmented data matrix.  However, as
> I recall, you only have a fixed number of samples to work with, so this
> might not help you.  2) You can try "stacking" (stacking horizontally)
> more than two overlapping data sets (instead of just one (vertically) as
> described by many textbooks).  This should make your noise subspace
> bigger, thereby "absorbing" more of the noise from your signal (I don't
> have a good technically explanation for this, but I think this will
> work).  I have done this before for a related problem so you can try it
> here as well.  Pls. let me know how it goes.
>
> 2. Yeah, doing SVD incurs O(n^3) computational complexity (for a square
> n X n matrix), so it's pretty slow, esp. for big matrices.
>
> 3 and 4. I haven't seen any esprit based algorithm for amplitude and
> phase estimation either, but I have a feeling this can be done.  Perhaps
> I will look into it :) .  For amplitude estimation, I have seen a least
> squares based algorithm.  If you have access to IEEE papers, you can
> look at:
>
> P. Stoica, H. Li and J. Li, "Amplitude Estimation of Sinusoidal Signals:
> Survey, New Results, and an Application", vol. 48(2), pp. 338-352, Feb.
> 2000.
>
> Hope this helps,
>
> cf
>
>
> Dan wrote:
> > Again, thank you all for the input that you have provided.  I've been
> > working with the ESPRIT algorithm, and it seems to be doing an
> > admirable job of separating the frequency components.  I just have a
> > few last second thoughts that I would appreciate a gut check on (if
> > you're not familiar with the term, it's a "yeah, that's about right" or
> > a "You're out of your mind!").
> > 1.  I'm assuming in the esprit method one obtains the highest accuracy
> > if you choose the number of shifts so that the matrix which you first
> > perform SVD on is as nearly square as possible.
> > 2.  I'm also assuming that the penalty associated is slower processing
> > of the SVD (I work in matlab, btw).
> > 3.  I have not found any reference to a version of esprit which
> > generates amplitudes so I am doing that as a separate step after esprit
> > gives me the frequencies.
> > 4.  I am also assuming that esprit cannot provide phase information,
> > likewise I'm extracting that from amplitude info (atan of overlap of
> > raw data with sin and cos of matched frequencies respectively.)  Does
> > anybody know of a better way to get the phase?  (You can even choose
> > whatever arbitrary reference you want.  See I can be generous!)
> >
> > Again, thank you all so much for the assistance.
> > 
> > Dan
> >

Hi Dan,

What I mean by stacking the columns horizontally is that giving the data 
matrix X, you may want to try something like:
     [   |         |        |        ]
     [   |         |        |        ]
X = [  x(1:5)   x(2:6)   x(3:7) ....]
     [   |         |        |        ]
     [   |         |        |        ]

where x(i:j) is a column vector (containing your data) with elements 
from index 1 to 5 (using matlab notation here).  Then the next column is 
a shift version of the previous one.  Such a matrix is called a Hankel 
matrix, where the elements on the antidiagonal is the same.  If you have 
something like this, this will increase the noise subspace.  After doing 
this, you might try to factorize this matrix, perhaps something like
X=ABC, such that the signal subspace vectors of X spans the same space 
as the columns in matrix A (that space being the signal subspace or the 
column space of X).  Then you might be able to obtain the matrix pencil 
relation you got before in the original ESPRIT algorithm.  It sounds 
very complicated, and I don't have all the math worked out yet, but you 
can take a stab at it if you want.

About the paper, I am glad you actually found it because I just realize 
I forgot to tell you which journal it is published in (IEEE Trans on 
Signal Processing).  But I guess you probably googled it.

Anyway, let me know how it goes, either it works or not.

Good luck,

cf

Dan wrote:
> cf,
> 
> Thanks for the info.  I'm not entirely clear as to what you mean by
> stacking overlapping data sets horizontally.  I believed that a
> requirement of this method was that the data had to be contiguous so
> that any adjacent columns of the augmented data matrix would be exactly
> one step out of phase with each other.  Are you suggesting taking two
> full data sets, creating the augmented matrix for each of them, and
> then concatenating them horizontally to make a double size augmented
> data matrix?  I would say that I clearly do not fully understand how
> the math is working here yet.
> 
> I am taking a look at the paper that you referenced.  It may be a while
> before I have the time to work all the way through it, but it does look
> promising.  Thanks,
> 
> Dan
> 
> cf wrote:
>> Hi Dan,
>>
>> 1. I don't think making the augmented data matrix (the matrix you
>> perform SVD on) as square as possible will make your estimates more
>> accurate.  To make your estimates more accurate, 1) you should collect
>> as much data as you can to construct augmented data matrix.  However, as
>> I recall, you only have a fixed number of samples to work with, so this
>> might not help you.  2) You can try "stacking" (stacking horizontally)
>> more than two overlapping data sets (instead of just one (vertically) as
>> described by many textbooks).  This should make your noise subspace
>> bigger, thereby "absorbing" more of the noise from your signal (I don't
>> have a good technically explanation for this, but I think this will
>> work).  I have done this before for a related problem so you can try it
>> here as well.  Pls. let me know how it goes.
>>
>> 2. Yeah, doing SVD incurs O(n^3) computational complexity (for a square
>> n X n matrix), so it's pretty slow, esp. for big matrices.
>>
>> 3 and 4. I haven't seen any esprit based algorithm for amplitude and
>> phase estimation either, but I have a feeling this can be done.  Perhaps
>> I will look into it :) .  For amplitude estimation, I have seen a least
>> squares based algorithm.  If you have access to IEEE papers, you can
>> look at:
>>
>> P. Stoica, H. Li and J. Li, "Amplitude Estimation of Sinusoidal Signals:
>> Survey, New Results, and an Application", vol. 48(2), pp. 338-352, Feb.
>> 2000.
>>
>> Hope this helps,
>>
>> cf
>>

cf:
How is this any different from the augmented data matrix that we have
to make for ESPRIT?  It looks to me like it is exactly the same
structure....  As for the factorization and later comments:  I'm going
to have to go dig out some references so that I can brush up on my
background... I've read it three times and I'm still not sure what
you're suggesting.  I'll try to take a look at it over the weekend.

Dan

cf wrote:
> Hi Dan,
>
> What I mean by stacking the columns horizontally is that giving the data
> matrix X, you may want to try something like:
>      [   |         |        |        ]
>      [   |         |        |        ]
> X = [  x(1:5)   x(2:6)   x(3:7) ....]
>      [   |         |        |        ]
>      [   |         |        |        ]
>
> where x(i:j) is a column vector (containing your data) with elements
> from index 1 to 5 (using matlab notation here).  Then the next column is
> a shift version of the previous one.  Such a matrix is called a Hankel
> matrix, where the elements on the antidiagonal is the same.  If you have
> something like this, this will increase the noise subspace.  After doing
> this, you might try to factorize this matrix, perhaps something like
> X=ABC, such that the signal subspace vectors of X spans the same space
> as the columns in matrix A (that space being the signal subspace or the
> column space of X).  Then you might be able to obtain the matrix pencil
> relation you got before in the original ESPRIT algorithm.  It sounds
> very complicated, and I don't have all the math worked out yet, but you
> can take a stab at it if you want.
>
> About the paper, I am glad you actually found it because I just realize
> I forgot to tell you which journal it is published in (IEEE Trans on
> Signal Processing).  But I guess you probably googled it.
>
> Anyway, let me know how it goes, either it works or not.
>
> Good luck,
>
> cf
>
> Dan wrote:
> > cf,
> >
> > Thanks for the info.  I'm not entirely clear as to what you mean by
> > stacking overlapping data sets horizontally.  I believed that a
> > requirement of this method was that the data had to be contiguous so
> > that any adjacent columns of the augmented data matrix would be exactly
> > one step out of phase with each other.  Are you suggesting taking two
> > full data sets, creating the augmented matrix for each of them, and
> > then concatenating them horizontally to make a double size augmented
> > data matrix?  I would say that I clearly do not fully understand how
> > the math is working here yet.
> >
> > I am taking a look at the paper that you referenced.  It may be a while
> > before I have the time to work all the way through it, but it does look
> > promising.  Thanks,
> >
> > Dan
> >
> > cf wrote:
> >> Hi Dan,
> >>
> >> 1. I don't think making the augmented data matrix (the matrix you
> >> perform SVD on) as square as possible will make your estimates more
> >> accurate.  To make your estimates more accurate, 1) you should collect
> >> as much data as you can to construct augmented data matrix.  However, as
> >> I recall, you only have a fixed number of samples to work with, so this
> >> might not help you.  2) You can try "stacking" (stacking horizontally)
> >> more than two overlapping data sets (instead of just one (vertically) as
> >> described by many textbooks).  This should make your noise subspace
> >> bigger, thereby "absorbing" more of the noise from your signal (I don't
> >> have a good technically explanation for this, but I think this will
> >> work).  I have done this before for a related problem so you can try it
> >> here as well.  Pls. let me know how it goes.
> >>
> >> 2. Yeah, doing SVD incurs O(n^3) computational complexity (for a square
> >> n X n matrix), so it's pretty slow, esp. for big matrices.
> >>
> >> 3 and 4. I haven't seen any esprit based algorithm for amplitude and
> >> phase estimation either, but I have a feeling this can be done.  Perhaps
> >> I will look into it :) .  For amplitude estimation, I have seen a least
> >> squares based algorithm.  If you have access to IEEE papers, you can
> >> look at:
> >>
> >> P. Stoica, H. Li and J. Li, "Amplitude Estimation of Sinusoidal Signals:
> >> Survey, New Results, and an Application", vol. 48(2), pp. 338-352, Feb.
> >> 2000.
> >>
> >> Hope this helps,
> >>
> >> cf
> >>

On Mon, 16 Oct 2006 15:16:54 -0500, "Jeff Caunter"
<jeffcaunter@sparkysworld.co.uk> wrote:

>If the two frequencies are so close together that they cannot be separately
>resolved by Fourier Transforms, then the hint I suggested has no chance of
>working.
>
>Can you give more details - Sample rate, range of frequency over which the
>two signals span, and an indication of the frequency difference between the
>two signals?
>
>Jeff

Hi Jeff,
   I haven't carefully read all the posts 
in this thread, but maybe the original poster
should search the Internet for:

"Iterative Weighted Phase Averager"

as a method for estimating the center frequency 
of narrowand spectral components.
Jus a suggestion.

[-Rick-]

Rick,
I ran the search that you suggested, and came up with just one useful
hit.  It was a reference to paper by Santamaria et al.  I've gotten the
paper and read through it, but I'm not finding much that is of
significant help, as much of the discussion of how to implement iwpa is
vague.  Do you know of additional resources on this topic?

Thanks,
Dan

Rick Lyons wrote:
> On Mon, 16 Oct 2006 15:16:54 -0500, "Jeff Caunter"
> <jeffcaunter@sparkysworld.co.uk> wrote:
>
> >If the two frequencies are so close together that they cannot be separately
> >resolved by Fourier Transforms, then the hint I suggested has no chance of
> >working.
> >
> >Can you give more details - Sample rate, range of frequency over which the
> >two signals span, and an indication of the frequency difference between the
> >two signals?
> >
> >Jeff
>
> Hi Jeff,
>    I haven't carefully read all the posts
> in this thread, but maybe the original poster
> should search the Internet for:
>
> "Iterative Weighted Phase Averager"
>
> as a method for estimating the center frequency
> of narrowand spectral components.
> Jus a suggestion.
> 
> [-Rick-]