Question about Gaussian window and parabolic interpolation

I have done some experiments with a Gaussian window applied to real
signals. The signals were pure sinusoids of a fixed frequency and
phase.

I used this window:

w(k) = e^( -(1/2) * (k - (N-1)/2)^2 / (a*(N-1)/2)^2 )

with N=32, k = 0, ..., N-1 and a = 0.265391249974779

I multiplied my signal with the window and did the DFT.
Then I calculated the magnitudes around the peak.
I took the natural logarithm of the magnitudes and did parabolic
interpolation to find the location of the true peak.
The results were almost exact.

Now I have this question:
I found a = 0.265391249974779 by minimizing the MSE between the real
and the estimated frequency.
That value gave the most exact results.
I believe it is supposed to have something to do with std dev, but can
someone explain why this is the correct value to use?
And maybe give an analytical expression for the number.

Michael

On Thursday, March 9, 2017 at 5:59:01 PM UTC-5, Michael Plet wrote:
> I have done some experiments with a Gaussian window applied to real
> signals. The signals were pure sinusoids of a fixed frequency and
> phase.
> 
> I used this window:
> 
> w(k) = e^( -(1/2) * (k - (N-1)/2)^2 / (a*(N-1)/2)^2 )
> 
> with N=32, k = 0, ..., N-1 and a = 0.265391249974779
> 
> I multiplied my signal with the window and did the DFT.
> Then I calculated the magnitudes around the peak.
> I took the natural logarithm of the magnitudes and did parabolic
> interpolation to find the location of the true peak.
> The results were almost exact.
> 
> Now I have this question:
> I found a = 0.265391249974779 by minimizing the MSE between the real
> and the estimated frequency.
> That value gave the most exact results.
> I believe it is supposed to have something to do with std dev, but can
> someone explain why this is the correct value to use?
> And maybe give an analytical expression for the number.
> 
> Michael

i have no idea why 0.265... is better than a smaller number for "a".  the Gaussian window gets down to a very small value (0.0008259) before it is truncated to zero.  but that's hardly 10 bits.  i would think a smaller value for "a" would be closer to an ideal Gaussian.  the closer to "ideal" means smaller sidelobes.  a perfect Gaussian has no side lobes.

and, it appears that you understand that Fourier Transform of a Gaussian is a Gaussian.  and that doing the log magnitude should get you the quadratic function that is in the exponent of a Gaussian.  so quadratic interpolation of the log magnitude to get the peak location should be perfect.

r b-j

Am 10.03.17 um 06:41 schrieb robert bristow-johnson:
> On Thursday, March 9, 2017 at 5:59:01 PM UTC-5, Michael Plet wrote:
>> I have done some experiments with a Gaussian window applied to real
>> signals. The signals were pure sinusoids of a fixed frequency and
>> phase.
>>
>> I used this window:
>>
>> w(k) = e^( -(1/2) * (k - (N-1)/2)^2 / (a*(N-1)/2)^2 )
>>
>> with N=32, k = 0, ..., N-1 and a = 0.265391249974779

>> Now I have this question:
>> I found a = 0.265391249974779 by minimizing the MSE between the real
>> and the estimated frequency.
>> That value gave the most exact results.
>> I believe it is supposed to have something to do with std dev, but can
>> someone explain why this is the correct value to use?
>> And maybe give an analytical expression for the number.

I do no not believe that this is a universal constant. Have you tried 
different sample lengths N? I'd expect that it converges to some 
reasonable value for N->oo, maybe 0.25.

The standard form for a Gaussian is e^(-x^2 / (2*sigma^2)) where sigma 
is the standard deviation. In your usage, the more you limit the 
Gaussian bandwidth, the poorer the frequency resolution will become due 
to the uncertainty principle. OTOH if you choose too small a, then the 
width of the Gaussian window surpasses the finite input, which 
approximates the window with a rectangular window, thereby disturbing 
your interpolation formula which expects a Gaussian peak shape.

Because of that, I believe that the constant you found is only valid for 
this special window length.

	Christian

On Fri, 10 Mar 2017 07:25:04 +0100, Christian Gollwitzer
<auriocus@gmx.de> wrote:

>Am 10.03.17 um 06:41 schrieb robert bristow-johnson:
>> On Thursday, March 9, 2017 at 5:59:01 PM UTC-5, Michael Plet wrote:
>>> I have done some experiments with a Gaussian window applied to real
>>> signals. The signals were pure sinusoids of a fixed frequency and
>>> phase.
>>>
>>> I used this window:
>>>
>>> w(k) = e^( -(1/2) * (k - (N-1)/2)^2 / (a*(N-1)/2)^2 )
>>>
>>> with N=32, k = 0, ..., N-1 and a = 0.265391249974779
>
>>> Now I have this question:
>>> I found a = 0.265391249974779 by minimizing the MSE between the real
>>> and the estimated frequency.
>>> That value gave the most exact results.
>>> I believe it is supposed to have something to do with std dev, but can
>>> someone explain why this is the correct value to use?
>>> And maybe give an analytical expression for the number.
>
>I do no not believe that this is a universal constant. Have you tried 
>different sample lengths N? I'd expect that it converges to some 
>reasonable value for N->oo, maybe 0.25.
>
>The standard form for a Gaussian is e^(-x^2 / (2*sigma^2)) where sigma 
>is the standard deviation. In your usage, the more you limit the 
>Gaussian bandwidth, the poorer the frequency resolution will become due 
>to the uncertainty principle. OTOH if you choose too small a, then the 
>width of the Gaussian window surpasses the finite input, which 
>approximates the window with a rectangular window, thereby disturbing 
>your interpolation formula which expects a Gaussian peak shape.
>
>Because of that, I believe that the constant you found is only valid for 
>this special window length.
>
>	Christian
>


Thanks. I have tried many different values of N from 2^3 to 2^14.
The same value of a is the best in all cases.
I guess that's not surprising since N is in my denominator:

(a*(N-1)/2)^2

I think that a decides the with of the Gaussian. Somehow this
particular value makes the Gaussian in freq domain have the right
width.
I just don't understand why this value?

Michael

On Thu, 9 Mar 2017 21:41:40 -0800 (PST), robert bristow-johnson
<rbj@audioimagination.com> wrote:

>On Thursday, March 9, 2017 at 5:59:01 PM UTC-5, Michael Plet wrote:
>> I have done some experiments with a Gaussian window applied to real
>> signals. The signals were pure sinusoids of a fixed frequency and
>> phase.
>> 
>> I used this window:
>> 
>> w(k) = e^( -(1/2) * (k - (N-1)/2)^2 / (a*(N-1)/2)^2 )
>> 
>> with N=32, k = 0, ..., N-1 and a = 0.265391249974779
>> 
>> I multiplied my signal with the window and did the DFT.
>> Then I calculated the magnitudes around the peak.
>> I took the natural logarithm of the magnitudes and did parabolic
>> interpolation to find the location of the true peak.
>> The results were almost exact.
>> 
>> Now I have this question:
>> I found a = 0.265391249974779 by minimizing the MSE between the real
>> and the estimated frequency.
>> That value gave the most exact results.
>> I believe it is supposed to have something to do with std dev, but can
>> someone explain why this is the correct value to use?
>> And maybe give an analytical expression for the number.
>> 
>> Michael
>
>i have no idea why 0.265... is better than a smaller number for "a".  the Gaussian window gets down to a very small value (0.0008259) before it is truncated to zero.  but that's hardly 10 bits.  i would think a smaller value for "a" would be closer to an ideal Gaussian.  the closer to "ideal" means smaller sidelobes.  a perfect Gaussian has no side lobes.
>
>and, it appears that you understand that Fourier Transform of a Gaussian is a Gaussian.  and that doing the log magnitude should get you the quadratic function that is in the exponent of a Gaussian.  so quadratic interpolation of the log magnitude to get the peak location should be perfect.
>
>r b-j

Yes, 0.265... is better than a smaller (or a bigger) value.
And the results are quite accurate. Here is a test with N=32 and
frequencies 3.0 - 3.9:

Real      Parabolic estimate
3	3,00017992565646
3,1	3,10006318924732
3,2	3,20004028033627
3,3	3,30002322785603
3,4	3,39996269738451
3,5	3,49988325820929
3,6	3,59989685072534
3,7	3,69993366240555
3,8	3,79995380779119

Michael

>  so quadratic interpolation
>of the log magnitude to get the peak location should be perfect.
>>
>>r b-j
>
>Yes, 0.265... is better than a smaller (or a bigger) value.
>And the results are quite accurate. Here is a test with N2 and
>frequencies 3.0 - 3.9:
>
>Real      Parabolic estimate
>3	3,00017992565646
>3,1	3,10006318924732
>3,2	3,20004028033627
>3,3	3,30002322785603
>3,4	3,39996269738451
>3,5	3,49988325820929
>3,6	3,59989685072534
>3,7	3,69993366240555
>3,8	3,79995380779119
>
>Michael

Michael,

The formula in r b-j's paper is based on a complex signal.  I think you
should try the same approach with a complex signal instead.  This may
account for a lot of the inaccuracy you are getting.  You will notice in
the data that the results are generally a little better as you move away
from the DC bin.

This approach is based on the Gaussian being the Eigenfunction of the
Fourier transform in the continuous case.  Thus, you should expect that a
Gaussian signal should be close to an Eigenvector of the DFT matrix.  If
you do a search on this you will find a lot of papers talking about
finding Gaussian like Eigenvectors.

You can reframe the operation in a slightly different manner.  Consider
the Gaussian your signal and the complex signal as a frequency shift. 
Doing the DFT on the straight up Gaussian would be the same as a windowed
constant signal or a frequency shift of zero.  The results are the same,
you should get an approximate Gaussian centered at the DC bin.  Ideally,
all the imaginary values in the DFT should be zero.

Anyway, because of symmetry, a parabolic fit at the DC bin should find a
peak at the DC bin, which is what you would expect.

Any complex signal with an integer frequency should shift the peak by the
frequency value and you should end up centered on a bin and once again get
an exact answer.

When you have a non-integer frequency is the case when the trueness of the
fit to a parabola matters.  As I mentioned in the other thread you have
two other sources of inaccuracy that seem to be counter balancing one
another.  If you make your "a" smaller, then you get more of the Bell
curve and reduce the inaccuracy due to the tails being smaller.  However,
you also make the sides of the bell steeper which increases the inaccuracy
between a summation and an integration.  I would have expected your value
to be dependent on N as well.

Just some thoughts from your friendly neighborhood pendant hobbyist.

Ced
---------------------------------------
Posted through http://www.DSPRelated.com

On Fri, 10 Mar 2017 09:33:17 -0600, "Cedron" <103185@DSPRelated>
wrote:

>>  so quadratic interpolation
>>of the log magnitude to get the peak location should be perfect.
>>>
>>>r b-j
>>
>>Yes, 0.265... is better than a smaller (or a bigger) value.
>>And the results are quite accurate. Here is a test with N2 and
>>frequencies 3.0 - 3.9:
>>
>>Real      Parabolic estimate
>>3	3,00017992565646
>>3,1	3,10006318924732
>>3,2	3,20004028033627
>>3,3	3,30002322785603
>>3,4	3,39996269738451
>>3,5	3,49988325820929
>>3,6	3,59989685072534
>>3,7	3,69993366240555
>>3,8	3,79995380779119
>>
>>Michael
>
>Michael,
>
>The formula in r b-j's paper is based on a complex signal.  I think you
>should try the same approach with a complex signal instead.  This may
>account for a lot of the inaccuracy you are getting.  You will notice in
>the data that the results are generally a little better as you move away
>from the DC bin.
>
>This approach is based on the Gaussian being the Eigenfunction of the
>Fourier transform in the continuous case.  Thus, you should expect that a
>Gaussian signal should be close to an Eigenvector of the DFT matrix.  If
>you do a search on this you will find a lot of papers talking about
>finding Gaussian like Eigenvectors.
>
>You can reframe the operation in a slightly different manner.  Consider
>the Gaussian your signal and the complex signal as a frequency shift. 
>Doing the DFT on the straight up Gaussian would be the same as a windowed
>constant signal or a frequency shift of zero.  The results are the same,
>you should get an approximate Gaussian centered at the DC bin.  Ideally,
>all the imaginary values in the DFT should be zero.
>
>Anyway, because of symmetry, a parabolic fit at the DC bin should find a
>peak at the DC bin, which is what you would expect.
>
>Any complex signal with an integer frequency should shift the peak by the
>frequency value and you should end up centered on a bin and once again get
>an exact answer.
>
>When you have a non-integer frequency is the case when the trueness of the
>fit to a parabola matters.  As I mentioned in the other thread you have
>two other sources of inaccuracy that seem to be counter balancing one
>another.  If you make your "a" smaller, then you get more of the Bell
>curve and reduce the inaccuracy due to the tails being smaller.  However,
>you also make the sides of the bell steeper which increases the inaccuracy
>between a summation and an integration.  I would have expected your value
>to be dependent on N as well.
>
>Just some thoughts from your friendly neighborhood pendant hobbyist.
>
>Ced
>---------------------------------------
>Posted through http://www.DSPRelated.com


Thanks Cedron. I don't think I'm getting inaccuracy. On the contrary I
am impressed with the accuracy. I can' t be perfect because the
Gaussian extends to -infinity and +infinity.

If the standard form of the Gaussian is

e^(-x^2 / (2*sigma^2)), where sigma is std dev for the Gaussian

and my window is

e^( -(1/2) * (k - (N-1)/2)^2 / (a*(N-1)/2)^2 )

that means that my "a" is

a = 2*sigma / (N-1)

So it does depend on N.

I am not interested in using complex signals. Only in explaining why
the particular value for "a" is the best.
Maybe it is, as you (and others) say, a question of balance.
A balance of getting the bell curve wide enough to get accuracy, but
not so wide that to high values are cut off at the ends.
I just can't help thinking that there must be a mathematical way of
finding "a".

Michael

>
>Thanks Cedron. I don't think I'm getting inaccuracy. On the contrary I
>am impressed with the accuracy. I can' t be perfect because the
>Gaussian extends to -infinity and +infinity.
>
>If the standard form of the Gaussian is
>
>e^(-x^2 / (2*sigma^2)), where sigma is std dev for the Gaussian
>
>and my window is
>
>e^( -(1/2) * (k - (N-1)/2)^2 / (a*(N-1)/2)^2 )
>
>that means that my "a" is
>
>a = 2*sigma / (N-1)
>
>So it does depend on N.
>
>I am not interested in using complex signals. Only in explaining why
>the particular value for "a" is the best.
>Maybe it is, as you (and others) say, a question of balance.
>A balance of getting the bell curve wide enough to get accuracy, but
>not so wide that to high values are cut off at the ends.
>I just can't help thinking that there must be a mathematical way of
>finding "a".
>
>Michael

Inaccuracy is a relative term, just meaning not exact.  Another way to say
it is that you have a very good approximation.

Your "a" may depend on N, but the rescaled chunk of the bell curve you are
using remains the same, i.e. from about -3.768 to 3.768 on the standard
curve (mean = 0, std dev = 1).

If the balance as I stated it is true, then I would expect the balance to
shift as the summation vs. integration error shrank with a larger N, but
the size of the tails you chop off remain relatively the same.  It is
complicated by the fact that you are then passing it through the DFT and
using a parabola fit.

The formula is based on a complex signal.  Since a real signal can be
thought of as the sum of two complex signals when you use this formula on
a real signal it is the same as if you had another signal interfering with
it some bin distance away.  So just for testing purposes, using a complex
signal will eliminate one source of inaccuracy.

This is like when Jacobsen applied your three bin magnitude exact
frequency formula for real signals to complex signals and claimed it
wasn't as good as his estimator.  

Finding an analytic formula will likely prove to be difficult.  At least
it will include a finite integration of the bell curve which will have to
be left as is because there is no analytic solution for that.

If the best "a" remains constant for any N in the complex signal case then
you have made a real contribution to the best practices of using this
formula. 

Ced
---------------------------------------
Posted through http://www.DSPRelated.com

On Friday, March 10, 2017 at 8:43:13 AM UTC-8, Michael Plet wrote:
> ...

> 
> I am not interested in using complex signals. Only in explaining why
> the particular value for "a" is the best.
> ...

> 
> Michael

So try putting your signal from 43.0 to 44.0 and see if you get the same result.

Dale B. Dalrymple

>
>If the balance as I stated it is true, then I would expect the balance
to
>shift as the summation vs. integration error shrank with a larger N, but
>the size of the tails you chop off remain relatively the same.  It is
>complicated by the fact that you are then passing it through the DFT and
>using a parabola fit.
>

Thinking this through a little more makes me change my mind a little bit. 
Even though you are chopping off relatively the same size tails you are
also looking a smaller portion of the top of the curve so the effect of
any error due to the distortion of the ideal parabola will also be
reduced.  So it is quite conceivable that the balance should remain at the
same spot as N increases.

I believe the vast portion of your current inaccuracy is due to using a
real signal instead of a complex signal.

Ced
---------------------------------------
Posted through http://www.DSPRelated.com