spectral estimation algorithm

dani_sb wrote:


> Firstly I tried the quadratic approach whici led to quite some bad
> estimates so I also tried the Gaussian ratio method
> Before I found the Newton stuff and this website I had an idea based on the
> wavelet transform but the problems encountered were about the same as with
> the FFT:

Yes, the problem is exactly the same as with the FFT: you don't have a 
clue and just keep trying various crap googled from the Internet. No 
wonder it doesn't work.

> I even thought of building up some sort of
> mother wavelet that would translate the frequency bands resulting in
> multiple distributions what would allow by comparison (or some algorithm)
> to estimate the tones parameters.

Wow. How about reading a basic textbook? It is so boring and not nearly 
as c00l as wavelets and fuzzy-wuzzy algorithms, but it works in the long 
prospective.


VLV

>For instance, I have some code that uses the same 3 tones used by the
>author:  1*cos(2pi*(10.5 Hz)*t) + .01*cos(2pi*(16 Hz)*t + .
>1*cos(2pi*(25.21 Hz)*t)  (the first 2 are the same as the ones in the
>well-known Harris windows paper).
>Using a 64 point FFT and a sample rate of 64, I compute the Gaussian
>window, apply it to the data, take the FFT, find (via 3 'if'
>statements) the largest peaks in the amplitude squared results, and
>estimate the frequency of the strongest tone.  Then, using that
>estimated frequency, I compute the DFT for that specific frequency and
>obtain amplitude and phase.  Then I subtract the estimated tone from
>the raw input data.  The process is then repeated for the second and
>third tones.  When computing the DFTs, I use a Hann window because it
>gives me better phase results than when using a Gaussian.
>My results are:
>est. frequency     amplitude       phase (radians)
>
>10.500000685    1.000000034     0.000000982
>
>25.210002870    0.099999972     -0.000007120
>
>15.999980676    0.009999999     0.000060729
>
>The above can be tweaked to squeeze out some extra digits of accuracy
>by reducing the computational noise (I used a recursive FFT for the
>above, and the DFT could be computed differently to reduce its error
>contribution).
>
>As I mentioned, I've also re-derived the equations so that I can use
>arbitrary bin spacing.  This is useful if, say, I want to compute the
>estimates for a single tone at, oh ...10.5 Hz.  I just use 2 DFTs: one
>at a point below 10.5 Hz, and another at a point higher than 10.5 Hz
>(how far apart depends on how wide I make the beamwidth of the window
>in the frequency domain).  Then Gaussian window the data, compute the
>2 DFTs, estimate frequency, and then compute one more DFT at the
>estimated frequency to compute amplitude/phase.  And the windowing can
>actually be done for free by using a modified look-up table for the
>DFT sine/cosine parts.  So the total processing is just 3 look-up
>based DFTs plus the natural log calculation when doing the frequency
>estimate.

Correlating the peaks indexes from the squared magnitude FFT spectrum with
the 32 points F[i] and A[i] arrays, I found : 

index          F           A
11          10.4998      1.0001       (or 10    10.5    1.00001)     
16          15.9855      0.0100352
25          25.2093      0.100011
 

In order to obtain the phase information I should compute DFT for the 10.5
Hz, 16 Hz and 25.21 Hz specific frequencies. For the integer one it should
work directly and the non-integer ones should work by the up-mentioned
algorithm 'one at a point below 10.5 Hz, and another at a point higher than
10.5 Hz'. This last method is not yet clear for me along with computing a
DFT at a specific frequency which I believe it's not : 
f=10.5 Hz
         N-1
Re=(2/N)*sum(Data(i)*cos(2*f*i*pi/N));
         i=0

         N-1
Im=(2/N)*sum(Data(i)*sin(2*f*i*pi/N));
         i=0

Mag=sqrt(Re^2+Im^2);

"Of all the idiots, the most obnoxious are the ones that are not 
completely devoid of mind"

Francous De La Rochefocauld





dani_sb wrote:

>>For instance, I have some code that uses the same 3 tones used by the
>>author:  1*cos(2pi*(10.5 Hz)*t) + .01*cos(2pi*(16 Hz)*t + .
>>1*cos(2pi*(25.21 Hz)*t)  (the first 2 are the same as the ones in the
>>well-known Harris windows paper).
>>Using a 64 point FFT and a sample rate of 64, I compute the Gaussian
>>window, apply it to the data, take the FFT, find (via 3 'if'
>>statements) the largest peaks in the amplitude squared results, and
>>estimate the frequency of the strongest tone.  Then, using that
>>estimated frequency, I compute the DFT for that specific frequency and
>>obtain amplitude and phase.  Then I subtract the estimated tone from
>>the raw input data.  The process is then repeated for the second and
>>third tones.  When computing the DFTs, I use a Hann window because it
>>gives me better phase results than when using a Gaussian.
>>My results are:
>>est. frequency     amplitude       phase (radians)
>>
>>10.500000685    1.000000034     0.000000982
>>
>>25.210002870    0.099999972     -0.000007120
>>
>>15.999980676    0.009999999     0.000060729
>>
>>The above can be tweaked to squeeze out some extra digits of accuracy
>>by reducing the computational noise (I used a recursive FFT for the
>>above, and the DFT could be computed differently to reduce its error
>>contribution).
>>
>>As I mentioned, I've also re-derived the equations so that I can use
>>arbitrary bin spacing.  This is useful if, say, I want to compute the
>>estimates for a single tone at, oh ...10.5 Hz.  I just use 2 DFTs: one
>>at a point below 10.5 Hz, and another at a point higher than 10.5 Hz
>>(how far apart depends on how wide I make the beamwidth of the window
>>in the frequency domain).  Then Gaussian window the data, compute the
>>2 DFTs, estimate frequency, and then compute one more DFT at the
>>estimated frequency to compute amplitude/phase.  And the windowing can
>>actually be done for free by using a modified look-up table for the
>>DFT sine/cosine parts.  So the total processing is just 3 look-up
>>based DFTs plus the natural log calculation when doing the frequency
>>estimate.
> 
> 
> Correlating the peaks indexes from the squared magnitude FFT spectrum with
> the 32 points F[i] and A[i] arrays, I found : 
> 
> index          F           A
> 11          10.4998      1.0001       (or 10    10.5    1.00001)     
> 16          15.9855      0.0100352
> 25          25.2093      0.100011
>  
> 
> In order to obtain the phase information I should compute DFT for the 10.5
> Hz, 16 Hz and 25.21 Hz specific frequencies. For the integer one it should
> work directly and the non-integer ones should work by the up-mentioned
> algorithm 'one at a point below 10.5 Hz, and another at a point higher than
> 10.5 Hz'. This last method is not yet clear for me along with computing a
> DFT at a specific frequency which I believe it's not : 
> f=10.5 Hz
>          N-1
> Re=(2/N)*sum(Data(i)*cos(2*f*i*pi/N));
>          i=0
> 
>          N-1
> Im=(2/N)*sum(Data(i)*sin(2*f*i*pi/N));
>          i=0
> 
> Mag=sqrt(Re^2+Im^2);
> 
> 
> 
> 

>
>Yes, the problem is exactly the same as with the FFT: you don't have a 
>clue and just keep trying various crap googled from the Internet. No 
>wonder it doesn't work.
>
I am not an expert in DSP - that's true ...but DSP is not my basic activity
(I am an Electrical Engineer). I need to accomplish this harmonic
identification problem (among many others, more or less DSP related)for 
power quality issues. However, why do you think that 'google crap' isn't
good? 

The 'Newton-Raphson nonsense' that I found in labview examples works fine
and I could have used it in my project. I didn't do it just because I
wanted to understand what I'm doing.

>
>Wow. How about reading a basic textbook? It is so boring and not nearly 
>as c00l as wavelets and fuzzy-wuzzy algorithms, but it works in the long 
>prospective.
>
>
>VLV
>

A basic textbook always helps - I agree . It's not boring but in my case
it's a bit - time consuming. I prefer to read the theory more focused on my
issue and not the whole science. 

Daniel

>But let's burn one bridge at a time.  Get that 'pi' term defined and
>move the power spectrum connection in LabView to before the mag
>squared part.  The McEachern results (outputs 8 to 11) show:
>

I took the pi out from the code because the C node from labview already has
that defined. However,the error in my code was sorted after I figured out
that 'log' in C++ = 'ln' in Labview.

http://img522.imageshack.us/img522/7945/gaussianratio.png

The power spectrum box is a simple plot named "power spectrum" by me. It
plots the array elements against their index. 


>For instance, I have some code that uses the same 3 tones used by the
>author:  1*cos(2pi*(10.5 Hz)*t) + .01*cos(2pi*(16 Hz)*t + .
>1*cos(2pi*(25.21 Hz)*t)  (the first 2 are the same as the ones in the
>well-known Harris windows paper).
>Using a 64 point FFT and a sample rate of 64, I compute the Gaussian
>window, apply it to the data, take the FFT, find (via 3 'if'
>statements) the largest peaks in the amplitude squared results, and
>estimate the frequency of the strongest tone.  Then, using that
>estimated frequency, I compute the DFT for that specific frequency and
>obtain amplitude and phase.  Then I subtract the estimated tone from
>the raw input data.  The process is then repeated for the second and
>third tones.  When computing the DFTs, I use a Hann window because it
>gives me better phase results than when using a Gaussian.

>My results are:
>est. frequency     amplitude       phase (radians)
>
>10.500000685    1.000000034     0.000000982
>
>25.210002870    0.099999972     -0.000007120
>
>15.999980676    0.009999999     0.000060729
>
>The above can be tweaked to squeeze out some extra digits of accuracy
>by reducing the computational noise (I used a recursive FFT for the
>above, and the DFT could be computed differently to reduce its error
>contribution).
>
>As I mentioned, I've also re-derived the equations so that I can use
>arbitrary bin spacing.  This is useful if, say, I want to compute the
>estimates for a single tone at, oh ...10.5 Hz.  I just use 2 DFTs: one
>at a point below 10.5 Hz, and another at a point higher than 10.5 Hz
>(how far apart depends on how wide I make the beamwidth of the window
>in the frequency domain).  Then Gaussian window the data, compute the
>2 DFTs, estimate frequency, and then compute one more DFT at the
>estimated frequency to compute amplitude/phase.  And the windowing can
>actually be done for free by using a modified look-up table for the
>DFT sine/cosine parts.  So the total processing is just 3 look-up
>based DFTs plus the natural log calculation when doing the frequency
>estimate.

Considering this example, I should compute DFT for the 10.5 Hz, 16 Hz and
25.21 Hz specific frequencies. For the integer one it should work directly
and the non-integer ones should work by the up-mentioned algorithm 'one
at a point below 10.5 Hz, and another at a point higher than 10.5 Hz'



Correlating the peaks indexes from the squared magnitude FFT spectrum with
the 32 points F[i] and A[i] arrays, I found : 

index          F           A
11          10.4998      1.0001       (or 10    10.5    1.00001)     
16          15.9855      0.0100352
25          25.2093      0.100011
 

In order to obtain the phase information I should compute DFT for the 10.5
Hz, 16 Hz and 25.21 Hz specific frequencies. For the integer one it should
work directly and the non-integer ones should work by the up-mentioned
algorithm 'one at a point below 10.5 Hz, and another at a point higher than
10.5 Hz'. This last method is not yet clear for me along with computing a
DFT at a specific frequency which I believe it's not : 
f=10.5 Hz
         N-1
Re=(2/N)*sum(Data(i)*cos(2*f*i*pi/N));
         i=0

         N-1
Im=(2/N)*sum(Data(i)*sin(2*f*i*pi/N));
         i=0

Mag=sqrt(Re^2+Im^2);

Thanks
Daniel

P.S. Please ignore my previous message (on this topic) I pressed send by
mistake.

>
>"Of all the idiots, the most obnoxious are the ones that are not 
>completely devoid of mind"
>
>Francous De La Rochefocauld
>

"Idiots and lunatics see only their own wit."

same french guy 

Daniel

(wit is someone skilled in making witty remarks)

dani_sb wrote:

>However,the error in my code was sorted after I figured out
>that 'log' in C++ = 'ln' in Labview

Aha!  When I saw the log( ), I thought it was the same as in C++
(i.e.: natural log).
Your results are much better (correct results are always to be
preferable to the alternative).

>The power spectrum box is a simple plot named "power spectrum" by me. It
>plots the array elements against their index.

That still doesn't look right. When I compute a normalized power
spectrum with that data, I get flat line of -83 dB that starts growing
at the #6 point until it gets to the 11 bin, where its value is 0 dB.
Then it decreases until it gets to the 15/16 bins, where it shows a
weak visible peak of -40 dB.  Then decreases to -80 dB at the 20 bin,
and increases to -20 dB  at the 25/26 bins.  Usually you look at the
power spectrum to find the location of the peaks, and that tells you
that there's something at those points.  LabView doesn't seem to be
calculating what I would expect to see.

>Correlating the peaks indexes from the squared magnitude FFT spectrum with
>the 32 points F[i] and A[i] arrays, I found :
>[...]

You should always use the index that corresponds to calculation
corresponding to the local peak in your power spectrum. Output 10
corresponds to the natural log frequency calculation that uses FFT
output data from bins10 and 11 (in my code, that's where the first
peak is).  So output 10 after the frequency calculation is the one you
should be using.

With regard to your other comments/questions, I'm sorry to say that
I'll have to take a pass at explaining them.  I have no problem
whatsoever helping someone out in replicating an open literature
article, but everything else is new stuff that I developed over the
years, at a not inconsiderable amount of time and effort (e.g.: phase
calculations are not part of the original article; neither is using 2
DFTs or the modified estimation equations that must be used when doing
things that way, nor are different windows discussed, etc. etc.).  As
I mentioned, I have some 30+ pages in my book regarding this
technique, and I'm loathe to publicly disclose it.  The same goes for
some of the new DFT material as well.  Not to mention the software.

However, I can tell you this: 1) your DFT sums are wrong (remember,
the output of the FFT is complex, so your DFT should be one sum based
on complex valued inputs and complex twiddle factors, and 2) the FFT
computes at integer output points - you could get the same thing doing
a single frequency DFT for each integer output point (a slow way of
doing it, but it works).  The important thing is that the DFT need not
be restricted to integer bin numbers - you can compute it for any
frequency.  There's plenty of literature available for computing the
DFT, so I don't think you'll have any problem figuring that part of it
out.

Kevin McGee

On Oct 30, 12:29&#4294967295;am, kevin <kevinjmc...@netscape.net> wrote:
> dani_sb wrote:
> >However,the error in my code was sorted after I figured out
> >that 'log' in C++ = 'ln' in Labview
>
> Aha! &#4294967295;When I saw the log( ), I thought it was the same as in C++
> (i.e.: natural log).
> Your results are much better (correct results are always to be
> preferable to the alternative).
>
> >The power spectrum box is a simple plot named "power spectrum" by me. It
> >plots the array elements against their index.
>
> That still doesn't look right. When I compute a normalized power
> spectrum with that data, I get flat line of -83 dB that starts growing
> at the #6 point until it gets to the 11 bin, where its value is 0 dB.
> Then it decreases until it gets to the 15/16 bins, where it shows a
> weak visible peak of -40 dB. &#4294967295;Then decreases to -80 dB at the 20 bin,
> and increases to -20 dB &#4294967295;at the 25/26 bins. &#4294967295;Usually you look at the
> power spectrum to find the location of the peaks, and that tells you
> that there's something at those points. &#4294967295;LabView doesn't seem to be
> calculating what I would expect to see.
>
> >Correlating the peaks indexes from the squared magnitude FFT spectrum with
> >the 32 points F[i] and A[i] arrays, I found :
> >[...]
>
> You should always use the index that corresponds to calculation
> corresponding to the local peak in your power spectrum. Output 10
> corresponds to the natural log frequency calculation that uses FFT
> output data from bins10 and 11 (in my code, that's where the first
> peak is). &#4294967295;So output 10 after the frequency calculation is the one you
> should be using.
>
> With regard to your other comments/questions, I'm sorry to say that
> I'll have to take a pass at explaining them. &#4294967295;I have no problem
> whatsoever helping someone out in replicating an open literature
> article, but everything else is new stuff that I developed over the
> years, at a not inconsiderable amount of time and effort (e.g.: phase
> calculations are not part of the original article; neither is using 2
> DFTs or the modified estimation equations that must be used when doing
> things that way, nor are different windows discussed, etc. etc.). &#4294967295;As
> I mentioned, I have some 30+ pages in my book regarding this
> technique, and I'm loathe to publicly disclose it. &#4294967295;The same goes for
> some of the new DFT material as well. &#4294967295;Not to mention the software.
>
> However, I can tell you this: 1) your DFT sums are wrong (remember,
> the output of the FFT is complex, so your DFT should be one sum based
> on complex valued inputs and complex twiddle factors, and 2) the FFT
> computes at integer output points - you could get the same thing doing
> a single frequency DFT for each integer output point (a slow way of
> doing it, but it works). &#4294967295;The important thing is that the DFT need not
> be restricted to integer bin numbers - you can compute it for any
> frequency. &#4294967295;There's plenty of literature available for computing the
> DFT, so I don't think you'll have any problem figuring that part of it
> out.
>
> Kevin McGee


Actually, the DFT for the phase estimation is a real input only one.
You know the raw data, you know the frequency estimate, and now you
need to compute the DFT at the estimated frequency using those real
input numbers.  Preferably, you'd use a window on the data to get a
better phase estimate.

Kevin McGee

On Oct 30, 1:21&#4294967295;am, kevin <kevinjmc...@netscape.net> wrote:
> On Oct 30, 12:29&#4294967295;am, kevin <kevinjmc...@netscape.net> wrote:
>
>
>
>
>
> > dani_sb wrote:
> > >However,the error in my code was sorted after I figured out
> > >that 'log' in C++ = 'ln' in Labview
>
> > Aha! &#4294967295;When I saw the log( ), I thought it was the same as in C++
> > (i.e.: natural log).
> > Your results are much better (correct results are always to be
> > preferable to the alternative).
>
> > >The power spectrum box is a simple plot named "power spectrum" by me. It
> > >plots the array elements against their index.
>
> > That still doesn't look right. When I compute a normalized power
> > spectrum with that data, I get flat line of -83 dB that starts growing
> > at the #6 point until it gets to the 11 bin, where its value is 0 dB.
> > Then it decreases until it gets to the 15/16 bins, where it shows a
> > weak visible peak of -40 dB. &#4294967295;Then decreases to -80 dB at the 20 bin,
> > and increases to -20 dB &#4294967295;at the 25/26 bins. &#4294967295;Usually you look at the
> > power spectrum to find the location of the peaks, and that tells you
> > that there's something at those points. &#4294967295;LabView doesn't seem to be
> > calculating what I would expect to see.
>
> > >Correlating the peaks indexes from the squared magnitude FFT spectrum with
> > >the 32 points F[i] and A[i] arrays, I found :
> > >[...]
>
> > You should always use the index that corresponds to calculation
> > corresponding to the local peak in your power spectrum. Output 10
> > corresponds to the natural log frequency calculation that uses FFT
> > output data from bins10 and 11 (in my code, that's where the first
> > peak is). &#4294967295;So output 10 after the frequency calculation is the one you
> > should be using.
>
> > With regard to your other comments/questions, I'm sorry to say that
> > I'll have to take a pass at explaining them. &#4294967295;I have no problem
> > whatsoever helping someone out in replicating an open literature
> > article, but everything else is new stuff that I developed over the
> > years, at a not inconsiderable amount of time and effort (e.g.: phase
> > calculations are not part of the original article; neither is using 2
> > DFTs or the modified estimation equations that must be used when doing
> > things that way, nor are different windows discussed, etc. etc.). &#4294967295;As
> > I mentioned, I have some 30+ pages in my book regarding this
> > technique, and I'm loathe to publicly disclose it. &#4294967295;The same goes for
> > some of the new DFT material as well. &#4294967295;Not to mention the software.
>
> > However, I can tell you this: 1) your DFT sums are wrong (remember,
> > the output of the FFT is complex, so your DFT should be one sum based
> > on complex valued inputs and complex twiddle factors, and 2) the FFT
> > computes at integer output points - you could get the same thing doing
> > a single frequency DFT for each integer output point (a slow way of
> > doing it, but it works). &#4294967295;The important thing is that the DFT need not
> > be restricted to integer bin numbers - you can compute it for any
> > frequency. &#4294967295;There's plenty of literature available for computing the
> > DFT, so I don't think you'll have any problem figuring that part of it
> > out.
>
> > Kevin McGee
>
> Actually, the DFT for the phase estimation is a real input only one.
> You know the raw data, you know the frequency estimate, and now you
> need to compute the DFT at the estimated frequency using those real
> input numbers. &#4294967295;Preferably, you'd use a window on the data to get a
> better phase estimate.
>
> Kevin McGee- Hide quoted text -
>
> - Show quoted text -

dani_sb wrote:

>The power spectrum box is a simple plot named "power spectrum" by me. It >plots the array elements against their index.

OK. I get it now.  I should have read it more carefully.  I thought
that the power spectrum box was a LabView function that converted its
inputs to a power spectrum.  No problem - simple fix - just take
10*log10 (amplitude squared) before putting the data into the graph
(that's 10*log base 10).  This will be a non-normalized power
spectrum.  You should be able to see the tones in the power spectrum
much more clearly after that.

Kevin McGee