spectral estimation algorithm

Hi,

I am a newbie to DSP and I hope to get some help related to the following
issue:

Suppose I have a sampled voltage signal (sampled at fs=6400Hz):

x(t)=sin(2*pi*f*t)+0.5*sin(11.4*2*pi*f*t) ... f=50Hz, t=0:0.02/128:0.02

If I take a FFT over 10 periods I will be able to find the frequency
(570Hz), amplitude (0.5V) and phase (0rad) of my harmonic. The problem is
to find the same components from only one period.

The FFT spectrum produces leakage due to the noninteger harmonic
(interharmonic) and the method that I found (which I do not understand
completely unfortunately) is based on the Newton-Raphson root finding
algorithm. 
The signal is first windowed in the time domain with a 128 samples hanning
window f(t)=x(t).w(t) and F=FFT(f(t)) 
The algorithm takes 3 complex points from the spectrum : F[n-1] (at 500
Hz), F[n] (550 Hz) and F[n+1] (600Hz) where F[n] is the peak and the
following formula is given :

http://img21.imageshack.us/img21/4129/formulab.jpg

fn is and initial guess frequency w is the hanning complex window. w*
denotes the complex conjugate

Could someone please help me understand how was that formula built. 
I actually found it in a Labview example (Complex Multi-tone
Measurements.vi)

After having G(f), the problem is easy : a simple loop that would update
the guess frequency by:

f[n+1]=f[n]-G(f[n])/G'(f[n])

Any help is greatly appreciated.  Thanks.

Daniel

On Oct 25, 10:29=A0am, "dani_sb" <dan1_sb@n_o_s_p_a_m.yahoo.com> wrote:
> Hi,
>
> I am a newbie to DSP and I hope to get some help related to the following
> issue:
>
> Suppose I have a sampled voltage signal (sampled at fs=3D6400Hz):
>
> x(t)=3Dsin(2*pi*f*t)+0.5*sin(11.4*2*pi*f*t) ... f=3D50Hz, t=3D0:0.02/128:=
0.02
>
> If I take a FFT over 10 periods I will be able to find the frequency
> (570Hz), amplitude (0.5V) and phase (0rad) of my harmonic. The problem is
> to find the same components from only one period.
>
> The FFT spectrum produces leakage due to the noninteger harmonic
> (interharmonic) and the method that I found (which I do not understand
> completely unfortunately) is based on the Newton-Raphson root finding
> algorithm.
> The signal is first windowed in the time domain with a 128 samples hannin=
g
> window f(t)=3Dx(t).w(t) and F=3DFFT(f(t))
> The algorithm takes 3 complex points from the spectrum : F[n-1] (at 500
> Hz), F[n] (550 Hz) and F[n+1] (600Hz) where F[n] is the peak and the
> following formula is given :
>
> http://img21.imageshack.us/img21/4129/formulab.jpg
>
> fn is and initial guess frequency w is the hanning complex window. w*
> denotes the complex conjugate
>
> Could someone please help me understand how was that formula built.
> I actually found it in a Labview example (Complex Multi-tone
> Measurements.vi)
>
> After having G(f), the problem is easy : a simple loop that would update
> the guess frequency by:
>
> f[n+1]=3Df[n]-G(f[n])/G'(f[n])
>
> Any help is greatly appreciated. =A0Thanks.
>
> Daniel


I Googled it a little bit, but couldn't find a description of the
exact algorithm they're using.  I think it involves a quadratic fit,
based on the info I found below:

http://zone.ni.com/devzone/cda/tut/p/id/3770

http://lavag.org/topic/11240-waveform-peak-detection-vi-bogus-numbers/

I don't know offhand how the equations provided in your link relate to
what they're doing because I don't know how all the terms are
defined.  But it's easy to see that the equation for G(f) is just a
reordering of the ratios above it.  For the ratios shown, just
multiply both sides with the product of the two denominators; you end
up with G(f).  Now if you let k =3D 0 and M =3D 3, you can see that the
indexing works out.  I think the Gprime term is a derivative, but I'm
not sure.  And I don't know why they're using the conjugate of the
window (cross-correlation perhaps?).  And the denominator terms at the
top of the page certainly looks like a scaling factor.  I think the
'iterate on the guess' part of their algorithm is due to the Newton-
Raphson method.

At any rate, there are a lot of ways of interpolating between FFT bins
to estimate the location of a spectral peaks. You might, for example,
consider using:

https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks=
.html

At least it comes with and explanation.

There is a significantly better method using a Gaussian ratio:

http://www.edn.com/archives/1994/030394/05df1.htm

It's several orders of magnitude more accurate than other methods.
I've used it for many years and developed a lot of new things with it
(e.g.: subtracting estimates from the raw data to 'unmask' the
presence of weaker tones, using a more accurate Gaussian window, re-
deriving the equations to allow for arbitrary bin spacing, estimating
the rate of change of non-stationary signals, etc.).  I devote some
30+ pages to it in a book I wrote, and that represents a lot of my
time and effort.  So I'm kind of reluctant to post any of that here.

But I think you'll find that if you follow the example and code in the
above, you could probably get it right.  Maybe you'd try the Quadratic
approach first, after which you can try the Gaussian.

Kevin McGee

dani_sb wrote:

> Hi,
> 
> I am a newbie to DSP and I hope to get some help related to the following
> issue:
> 
> Suppose I have a sampled voltage signal (sampled at fs=6400Hz):
> 
> x(t)=sin(2*pi*f*t)+0.5*sin(11.4*2*pi*f*t) ... f=50Hz, t=0:0.02/128:0.02
> 
> If I take a FFT over 10 periods I will be able to find the frequency
> (570Hz), amplitude (0.5V) and phase (0rad) of my harmonic. The problem is
> to find the same components from only one period.
> 
> The FFT spectrum produces leakage due to the noninteger harmonic
> (interharmonic) and the method that I found (which I do not understand
> completely unfortunately) is based on the Newton-Raphson

[...]

Don't take FFT and don't do that Newton-Raphson nonsense, especially if 
you don't understand what it does. There is approximately a zillion of 
ways to accomplish the result depending on what exactly the signal is 
and what do you need.

> Any help is greatly appreciated.  Thanks.

Does your great appreciation imply any real content, or is it purely 
imaginary?

VLV

>I don't know offhand how the equations provided in your link relate to
>what they're doing because I don't know how all the terms are
>defined.  But it's easy to see that the equation for G(f) is just a
>reordering of the ratios above it.  For the ratios shown, just
>multiply both sides with the product of the two denominators; you end
>up with G(f).  Now if you let k =3D 0 and M =3D 3, you can see that the
>indexing works out.  I think the Gprime term is a derivative, but I'm
>not sure.  And I don't know why they're using the conjugate of the
>window (cross-correlation perhaps?).  And the denominator terms at the
>top of the page certainly looks like a scaling factor.  I think the
>'iterate on the guess' part of their algorithm is due to the Newton-
>Raphson method.

The G prime term is indeed a derivative (it comes from the Newton-Raphson
algorithm - which I do understand how it works). The cross correlation
might have a sense because I found out that the FFT spectrum is previously
shifted.

I did some further reading on this topic and the issue equation (the one
with the ratios) might be related to Parseval's equality which i found in
another topic here :

http://www.dsprelated.com/showmessage/34956/1.php

>Parseval's equality tells us
>the sum of the squared time sequence is related to
>the sum of the squared (magnitude) of the fourier sequence
>by a constant factor
>
>i.e.
>
>sum( abs(x).^2 ) == K * sum( abs(fft(x)).^2 )
>
>where K depends on the fft algorithm, and the length of the dft.
>
>-- Mark

abs(fft(x)).^2 would be the squared magnitude of the FFT spectrum and so if
my harmonic component "energy" is distributed in 3 fft bins (mag1) it
should have an equivalent distribution in 2 bins (mag2). The K might be the
denominator inverse since it's a constant value.

K1 * sum( mag1.^2 )=K2 * sum( mag2.^2 )    (looks a bit like the first
equation - but I'm not sure)



>
>At any rate, there are a lot of ways of interpolating between FFT bins
>to estimate the location of a spectral peaks. You might, for example,
>consider using:
>
>https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks=
>.html
>
>At least it comes with and explanation.
>
>There is a significantly better method using a Gaussian ratio:
>
>http://www.edn.com/archives/1994/030394/05df1.htm
>
>It's several orders of magnitude more accurate than other methods.
>I've used it for many years and developed a lot of new things with it
>(e.g.: subtracting estimates from the raw data to 'unmask' the
>presence of weaker tones, using a more accurate Gaussian window, re-
>deriving the equations to allow for arbitrary bin spacing, estimating
>the rate of change of non-stationary signals, etc.).  I devote some
>30+ pages to it in a book I wrote, and that represents a lot of my
>time and effort.  So I'm kind of reluctant to post any of that here.
>
>But I think you'll find that if you follow the example and code in the
>above, you could probably get it right.  Maybe you'd try the Quadratic
>approach first, after which you can try the Gaussian.
>
>Kevin McGee
>

I will first try to implement the Quadratic approach which looks quite
good.

Further on this topic I have a new question:
If my initial sampled signal would have two close harmonics, i.e. 

x(t)=sin(2*pi*f*t)+0.5*sin(11.4*2*pi*f*t)+0.3*sin(11*2*pi*f*t),

is there any method to identify both harmonics from one cycle? (to me it
looks a bit unlikely since both harmonics are spread in the same frequency
bins)

Thanks a lot !

Daniel

>
>Don't take FFT and don't do that Newton-Raphson nonsense, especially if 
>you don't understand what it does. There is approximately a zillion of 
>ways to accomplish the result depending on what exactly the signal is 
>and what do you need.
>

I will try other methods ... Thanks

>> Any help is greatly appreciated.  Thanks.
>
>Does your great appreciation imply any real content, or is it purely 
>imaginary?
>

The content is Complex. (So it might include some real part too.)

Daniel

dani_sb wrote:

> 
>>>Any help is greatly appreciated.  Thanks.
>>
>>Does your great appreciation imply any real content, or is it purely 
>>imaginary?
> 
> The content is Complex. (So it might include some real part too.)

My contact is at the web site in the signature. Please come up with the 
detailed description of the problem and what should be accomplished.


Vladimir Vassilevsky
DSP and Mixed Signal Design Consultant
http://www.abvolt.com

On Oct 25, 7:29&#4294967295;am, "dani_sb" <dan1_sb@n_o_s_p_a_m.yahoo.com> wrote:
> Hi,
>
> I am a newbie to DSP and I hope to get some help related to the following
> issue:
>
> Suppose I have a sampled voltage signal (sampled at fs=6400Hz):
>
> x(t)=sin(2*pi*f*t)+0.5*sin(11.4*2*pi*f*t) ... f=50Hz, t=0:0.02/128:0.02
>
> If I take a FFT over 10 periods I will be able to find the frequency
> (570Hz), amplitude (0.5V) and phase (0rad) of my harmonic. The problem is
> to find the same components from only one period.
>
> The FFT spectrum produces leakage due to the noninteger harmonic
> (interharmonic) and the method that I found (which I do not understand
> completely unfortunately) is based on the Newton-Raphson root finding
> algorithm.
> The signal is first windowed in the time domain with a 128 samples hanning
> window f(t)=x(t).w(t) and F=FFT(f(t))
> The algorithm takes 3 complex points from the spectrum : F[n-1] (at 500
> Hz), F[n] (550 Hz) and F[n+1] (600Hz) where F[n] is the peak and the
> following formula is given :
>
> http://img21.imageshack.us/img21/4129/formulab.jpg
>
> fn is and initial guess frequency w is the hanning complex window. w*
> denotes the complex conjugate
>
> Could someone please help me understand how was that formula built.
> I actually found it in a Labview example (Complex Multi-tone
> Measurements.vi)
>
> After having G(f), the problem is easy : a simple loop that would update
> the guess frequency by:
>
> f[n+1]=f[n]-G(f[n])/G'(f[n])
>
> Any help is greatly appreciated. &#4294967295;Thanks.

This appears to be a peak interpolation method.  Note that
you can get a similar kind of spectral interpolation by
simply zero padding your data (1 "cycle" worth) by large
enough factor (say 10X or 100X), computing the longer
FFT, and looking for the local maxima there.  Center the
data in the zero-padded FFT aperture if you want the phase
calculation to be a bit more intuitive.

The result would also be equivalent to interpolation using
the transform of your window function as the interpolation
kernel.  This interpolation is usually more accurate than
parabolic interpolation, unless the transform of the window
function you use is "close enough" to a parabola.


IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
 http://www.nicholson.com/rhn/dsp.html

dani_sb wrote:
[...]
Well, I agree with Vladimir that using the Newton-Raphson is not a
very good idea.  And there are other ways of getting much more
accurate results.
And just to avoid any confusion for anyone attempting to replicate the
results in the Gaussian ratio paper, just replace the
'GetRealPowerSpectrum code with a 64 point FFT and then compute and
plot the amplitudes squared - the plot will be similar to the author's
graph, but not normalized or in dB.
As I mentioned, I've used the technique ever since I read about it in
Mar. 1994, and have been very impressed with how well it works.  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.

>Further on this topic I have a new question:
>If my initial sampled signal would have two close harmonics, i.e.
>x(t)=sin(2*pi*f*t)+0.5*sin(11.4*2*pi*f*t)+0.3*sin(11*2*pi*f*t),
>is there any method to identify both harmonics from one cycle? (to me it
>looks a bit unlikely since both harmonics are spread in the same frequency
>bins)

Well, if your 'f' is still 50, then you'd have tones at 550 and 570
Hz.  Your 'bin width' and spacing is sample_rate/N, so you've got
6400/128 = 50, which means that one tone (550) is bin-centered, and
the other one (570) lies between it and the next bin up.  If you use
the minimal sample rate and allow N to be large enough, you can always
get enough separation between them to make it work.  But the moment
you throw in that 'identify both harmonics from one cycle'
requirement, it becomes more difficult.  There are in fact things you
can do (e.g: I've been able to use it for 1/2 cycle or 3/4 cycle
inputs), but it becomes more sensitive to noise.  And there wasn't
another very close-by interfering tone.

There may be some parametric model that you can use. But even then,
you've got that time-bandwith problem to deal with.  The more specific
you make your problem statement, the better you can select a
processing model to deal with it.

Kevin McGee

>And just to avoid any confusion for anyone attempting to replicate the
>results in the Gaussian ratio paper, just replace the
>'GetRealPowerSpectrum code with a 64 point FFT and then compute and
>plot the amplitudes squared - the plot will be similar to the author's
>graph, but not normalized or in dB.

Firstly I tried the quadratic approach whici led to quite some bad
estimates so I also tried the Gaussian ratio method but I am not sure if my
results are good because the link to the author's graph is broken (or the
problem is at my end). My implementation and results would be the following
: 

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


>As I mentioned, I've used the technique ever since I read about it in
>Mar. 1994, and have been very impressed with how well it works.  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.
>
>>Further on this topic I have a new question:
>>If my initial sampled signal would have two close harmonics, i.e.
>>x(t)=sin(2*pi*f*t)+0.5*sin(11.4*2*pi*f*t)+0.3*sin(11*2*pi*f*t),
>>is there any method to identify both harmonics from one cycle? (to me it
>>looks a bit unlikely since both harmonics are spread in the same
frequency
>>bins)
>
>Well, if your 'f' is still 50, then you'd have tones at 550 and 570
>Hz.  Your 'bin width' and spacing is sample_rate/N, so you've got
>6400/128 = 50, which means that one tone (550) is bin-centered, and
>the other one (570) lies between it and the next bin up.  If you use
>the minimal sample rate and allow N to be large enough, you can always
>get enough separation between them to make it work.  But the moment
>you throw in that 'identify both harmonics from one cycle'
>requirement, it becomes more difficult.  There are in fact things you
>can do (e.g: I've been able to use it for 1/2 cycle or 3/4 cycle
>inputs), but it becomes more sensitive to noise.  And there wasn't
>another very close-by interfering tone.
>
>There may be some parametric model that you can use. But even then,
>you've got that time-bandwith problem to deal with.  The more specific
>you make your problem statement, the better you can select a
>processing model to deal with it.
>
>Kevin McGee
>

The 'one cycle thing' is more like the ideal case. My system input is a
voltage or current waveform measured to the input of a nonlinear load (i.e.
an arc furnance) so the shorter the analysis window, the higher the chances
of analyzing a stationary piece of signal. The fundamental frequency is 50
Hz (with some small variations) and the best I can do in terms of
acquisition is 25600 S/s (512 samples/cycle).

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: one tone was placed in more than one frequency band and closed up
tones could not be separated. 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.

Daniel

dani_sb wrote:

>Firstly I tried the quadratic approach whici led to quite some bad
>estimates so I also tried the Gaussian ratio method but I am not sure if my
>results are good because the link to the author's graph is broken.

Yes, the link does appear to be broken.  And your results are off, but
I'm sure we can fix it.  First: what happened to pi?  The McEachern
code has a statement: pi = 3.1415926;  Your code is missing that, and
if  it's not defined and the compiler sets it to zero, then all your
window points are all = 1.  Your power spectrum results look very much
like an unwindowed 10.5 Hz tone.  Perhaps you should plot the window
to make sure that it's being computed correctly (it should show small
numbers at first, then increase towards the middle (i.e: 32 point of
the 64 total points), and then decrease towards the end point at 63).

Second, you probably want to move that connection to your power
spectrum box to be before the magnitude squared part.  My guess is
that LabView requires inputs to a power spectrum to be real/imaginary,
not mag/phase.  As I said, your 'power spectrum' looks very much like
an unwindowed tone at 10.5, and something very small scale around the
23-26 points, and nothing at 16 (but I'm not surprised about the
latter, since that tone is 40 dB down).

Also check to see if Labview scales their FFT results or not.  that
will affect amplitude estimates.

Your XY graph actually doesn't look too bad, but the amplitudes tell
me that someting is not being scaled properly.

>The 'one cycle thing' is more like the ideal case. My system input is a
>voltage or current waveform measured to the input of a nonlinear load (i.e.
>an arc furnance) so the shorter the analysis window, the higher the chances
>of analyzing a stationary piece of signal.

The 'one cycle' requirement coupled with the '2 close-by tones'  is a
very severe problem.  And non-stationarity signals can be difficult at
times, but there are methods for dealing with them.

>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: one tone was placed in more than one frequency band and closed up
>tones could not be separated.

Yes, that time-bandwidth problem is a tough one to deal with, no
matter what transform type you're dealing with.
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:

      frequency    amplitude
8    10.49958     0.999533
9    10.50005     1.00001
10  10.50005     1.000012
11  10.49977     1.000104

Kevin McGee