Question about Gaussian window and parabolic interpolation

>
>I forgot to say that I knew there was an inverse relationship between
>the width of the window and the width of it's DFT.
>
>So intuitively I thought the best "a" would correspond to a certain
>ratio of those two widths.
>
>Michael

Yep.  My playing around seems to indicate that the best ratio is 1:1.  I
was not rigorous about testing this.

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

On Sun, 12 Mar 2017 10:11:19 -0500, "Cedron" <103185@DSPRelated>
wrote:

>>
>>Thanks. ( N-1)/2 is the midpoint of the sequence, so I am surprised
>>that N/2 is better. E.g. if you have 8 numbers 0,...,7 the midpoint is
>>3.5.
>>So how does your window function look exactly?
>>
>>With the complex signals I testet I got exact results for a range of
>>"a" values.
>>But I am still curious about the real case.
>>
>>If I use a = 0.2 (instead of my constant 0.265...) I get these results
>>for a real frequency:
>>
>>3	3,00853376764318
>>3,1	3,09972285705032
>>3,2	3,19535649624848
>>3,3	3,29522039607512
>>3,4	3,39729649299645
>>3,5	3,49933200116124
>>3,6	3,59993967651388
>>3,7	3,70027774300613
>>3,8	3,80035465470425
>>
>>So even a small change in "a" make the results a lot worse.
>>
>>Michael
>
>You're welcome.  For the frequency determination the (N-1) vs N isn't that
>important.  I am using a steady signal with a constant frequency, so
>according to r b-j's equation I expect that the imaginary part of the
>first and second differences should be zero.
>
>With all the quibbling we've had around here about the meaning of the word
>"exact" it would be nicer if you said "numerically correct".  Since the
>formula is based on the assumption of a single complex tone your results
>are what is expected.  I haven't even started to test in the presence of
>noise yet which may have a big impact on what "a" should be.
>
>The real case is a special case of two complex tones.  Perhaps solving the
>general case of two tones will give you a better answer for the specific
>case.  Anyway, I ran the 3.3 real signal case with your same other
>parameters (other than the center) and got the results below.  Changing
>the center back to (N-1)/2 didn't change my frequency results much.
>
>It's strange that your results don't match mine.  What phase value did you
>use?
>
>Notice that the frequency determination at bin 6 is much more accurate
>than the one at bin 3.  This is due to the underlying parabola not being
>as affected by the other parabola at the -3.3 (28.7) frequency.
>
>I also ran the test with a 8.3 frequency real signal, which being near N/4
>is about as far as you can separate the two tones, and the interference
>becomes symmetric.  The results with the best frequency read was at bin 8
>as expected.
>
>Ced
>
>================================
>Code for the Window Definition
>
>
>   theCenter = N / 2.0;
>   theParabolaA = -0.5 / ( a * a * theCenter * theCenter );
>
>   x = n - theCenter;
>   theWindow[n] = exp( theParabolaA * x * x );
>
>============================
>Excerpted sample run
>
>Real   3.3 * 2Pi/N + Pi/4
>
>   N = 32
>
> This run's a = 0.2000
>
> Sigma Range is -5.0000 to 4.8438
>
>parabola A_time =  -0.0488
>
>
> k    Windowed DFT      Log of WDFT    First Diff    Second Diff
>--  ----------------   -------------  -------------  -------------
> 0   0.0046   0.0000   -5.39   0.00    0.00   1.75    4.39   0.00
> 1  -0.0074   0.0404   -3.19   1.75    1.49  -0.71   -1.42  -4.92
> 2   0.0141  -0.0882   -2.42  -1.41    0.55  -0.01   -0.46   6.30
> 3  -0.0193   0.1216   -2.09   1.73    0.12  -0.00   -0.40  -6.28
> 4   0.0178  -0.1124   -2.17  -1.41   -0.28  -0.00   -0.40   6.28
> 5  -0.0111   0.0700   -2.65   1.73   -0.67  -0.00   -0.39  -6.28
> 6   0.0047  -0.0294   -3.52  -1.41   -1.07  -0.00   -0.39   6.28
> 7  -0.0013   0.0083   -4.78   1.73   -1.46   0.00   -0.39  -6.28
> 8   0.0003  -0.0016   -6.44  -1.41   -1.86  -0.00   -0.40   6.28
>..
>22   0.0000   0.0000  -10.93   1.41    2.67  -0.05   -0.44  -6.17
>23  -0.0000  -0.0002   -8.49  -1.73    2.25   0.00   -0.39   6.28
>24   0.0003   0.0016   -6.44   1.41    1.86  -0.00   -0.40  -6.28
>25  -0.0013  -0.0083   -4.78  -1.73    1.46   0.00   -0.39   6.28
>26   0.0047   0.0294   -3.52   1.41    1.07  -0.00   -0.39  -6.28
>27  -0.0111  -0.0700   -2.65  -1.73    0.67  -0.00   -0.39   6.28
>28   0.0178   0.1124   -2.17   1.41    0.28  -0.00   -0.40  -6.28
>29  -0.0193  -0.1216   -2.09  -1.73   -0.12  -0.00   -0.40   6.28
>30   0.0141   0.0882   -2.42   1.41   -0.55  -0.01   -0.46  -6.30
>31  -0.0074  -0.0404   -3.19  -1.75   -1.49  -0.71   -1.42   4.92
>
>
> k         A          B         C        Freq       DFT Mag
>--     --------   -------    -------    --------    -------
> 0      2.1948     0.0000    -5.3883   -0.000000    0.0046
> 1     -0.7084     2.9032    -5.3883    2.049062    0.0410
> 2     -0.2287     1.4639    -4.4288    3.201164    0.0893
> 3     -0.1996     1.3188    -4.2546    3.303126    0.1231
> 4     -0.1976     1.3043    -4.2298    3.301068    0.1138
> 5     -0.1974     1.3029    -4.2267    3.300124    0.0708
> 6     -0.1974     1.3029    -4.2266    3.300092    0.0297
> 7     -0.1974     1.3024    -4.2252    3.299528    0.0084
> 8     -0.1976     1.3056    -4.2370    3.304009    0.0016
>..
>22     -0.2219    12.4282  -176.9734   28.008722    0.0000
>23     -0.1956    11.2456  -163.6754   28.748935    0.0002
>24     -0.1976    11.3392  -164.7754   28.695991    0.0016
>25     -0.1974    11.3289  -164.6490   28.700472    0.0084
>26     -0.1974    11.3307  -164.6716   28.699908    0.0297
>27     -0.1974    11.3308  -164.6736   28.699876    0.0708
>28     -0.1976    11.3394  -164.7911   28.698932    0.1138
>29     -0.1996    11.4571  -166.4677   28.696874    0.1231
>30     -0.2287    13.1698  -191.7227   28.798836    0.0893
>31     -0.7084    42.4364  -637.9190   29.950938    0.0410
>
>  parabola product =  3.1593
>
>---------------------------------------
>Posted through http://www.DSPRelated.com


My parameters:
    Amp = 1
    Phase = PI / 4

My signal:
    Factor = 2 * PI * Freq / N
    Real(k) = Amp * Sin(Factor * k + Phase)

>Thanks. ( N-1)/2 is the midpoint of the sequence, so I am surprised
>that N/2 is better. E.g. if you have 8 numbers 0,...,7 the midpoint is
>3.5.
..
>
>Michael


I am also surprised by the alternating nature of the bin values.  I am
going to try to do some theoretical work to figure out why.  The
difference equations I used compensate for this implicitly.  A little
quirk that is beneficial.

ForwardDifference  = f[k+1] - f[k]
BackwardDifference = f[k] - f[k-1]

FirstDifference = ( ForwardDifference + BackwardDifference ) / 2

SecondDifference = ForwardDifference - BackwardDifference

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

>
>
>My parameters:
>    Amp = 1
>    Phase = PI / 4
>
>My signal:
>    Factor = 2 * PI * Freq / N
>    Real(k) = Amp * Sin(Factor * k + Phase)

Okay, I feel better about it now.  I used Cos instead of Sin.

Here is a run with Sin:

k         A          B         C        Freq       DFT Mag
--     --------   -------    -------    --------    -------
 0      0.4373    -0.0000    -3.4467    0.000000    0.0318
 1      0.0788     0.3585    -3.4467   -2.275141    0.0493
 2     -0.1525     1.0523    -3.9093    3.450714    0.0894
 3     -0.1823     1.2014    -4.0882    3.295220    0.1195
 4     -0.1850     1.2203    -4.1205    3.298204    0.1109
 5     -0.1852     1.2224    -4.1252    3.299720    0.0711
 6     -0.1852     1.2226    -4.1258    3.299981    0.0315
 7     -0.1852     1.2226    -4.1257    3.299912    0.0096
 8     -0.1853     1.2230    -4.1274    3.300610    0.0020
 9     -0.1850     1.2190    -4.1101    3.293867    0.0003
10     -0.1878     1.2708    -4.3559    3.384139    0.0000

The values match.  Phew.

I reset my center to N/2 and created a formula to find the expected A for
the frequency parabola.  It turns out that the frequency calculations on
the bins where the calculated A matches the expected A give much better
results.  Perhaps this should be the criteria rather than the peak bin (in
the noiseless case).

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

>
>
>I am also surprised by the alternating nature of the bin values.

Not any more.  The Fourier transform works the same way forward and
backward.  So, if you multiply a signal with a complex tone you get a
shift in the frequency domain.  The reverse is if you multiply the
frequency domain with a complex "tone" you get a shift in the time domain,
aka a phase shift.

Alternating positive and negative ones corresponds to a complex tone with
a frequency of N/2.  So here is another reason that the window should be
centered at N/2.

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

On Sunday, March 12, 2017 at 2:59:12 PM UTC-4, Cedron wrote:
> >
> >
> >I am also surprised by the alternating nature of the bin values.
> 
> Not any more.  The Fourier transform works the same way forward and
> backward.  So, if you multiply a signal with a complex tone you get a
> shift in the frequency domain.  The reverse is if you multiply the
> frequency domain with a complex "tone" you get a shift in the time domain,
> aka a phase shift.
> 
> Alternating positive and negative ones corresponds to a complex tone with
> a frequency of N/2.  So here is another reason that the window should be
> centered at N/2.

it's centering the window in the time domain at N/2 that causes the alternating of the signs after the FFT.  this is why MATLAB has this function called fftshift().

r b-j

>On Sunday, March 12, 2017 at 2:59:12 PM UTC-4, Cedron wrote:
>> >
>> >
>> >I am also surprised by the alternating nature of the bin values.
>> 
>> Not any more.  The Fourier transform works the same way forward and
>> backward.  So, if you multiply a signal with a complex tone you get a
>> shift in the frequency domain.  The reverse is if you multiply the
>> frequency domain with a complex "tone" you get a shift in the time
>domain,
>> aka a phase shift.
>> 
>> Alternating positive and negative ones corresponds to a complex tone
with
>> a frequency of N/2.  So here is another reason that the window should
be
>> centered at N/2.
>
>it's centering the window in the time domain at N/2 that causes the
>alternating of the signs after the FFT.  this is why MATLAB has this
>function called fftshift().
>
>r b-j

Thanks, that's the "duh moment" I had yesterday, aka "No $#!^ Sherlock".

I don't use MATLAB (or SciLab) so I had to look up the documentation:

https://www.mathworks.com/help/matlab/ref/fftshift.html

This is not the same thing.  fftshift does a shift in the frequency
domain, it has nothing to do with alternating the signs which corresponds
to a shift in the time domain.

You can find the same concept in Vicanek's two bin frequency paper where
he calls it a "Centered DFT".  When I coded it, I called it "Rotating the
Bins" in the commentary.

For the case of N/2, the rotations turn out to be ones with alternating
signs making for a short cut implementation.  For the case of (N-1)/2 you
have to do a complex multiplication with the correct "tone" value for that
bin.  The rotations for two subsequent bins are almost 180 degrees apart,
i.e. their ratio is nearly -1.

What this means is that I have to alter my previous advice about centering
the window.  For even values of N the center should be at N/2 and for odd
values the center should be at (N-1)/2.  If you don't, after you apply the
centering fix you get a 180 degree phase discontinuity in the circularity
of your DFT.  If you are only using the log of the magnitudes, this
doesn't matter, but to me the whole point of using this type of window is
for the frequency sweep and amplitude ramp estimation capability which
depend on the imaginary part of the Log of the DFT.

Ced

   N = 16

parabola Middle A = -0.1963

         Middle a = 0.1995

 This run's a = 0.1995

 Sigma Range is -5.0133 to 4.6999

parabola A_time =  -0.1963

expected A_freq =  -0.1963

 k    Rotated WDFT      Log of RWDFT    First Diff    Second Diff
--  ----------------   -------------  -------------  -------------
 0   0.0046  -0.0291   -3.52  -1.41    1.30   0.00   -0.39  -0.00
 1   0.0138  -0.0874   -2.42  -1.41    0.90  -0.00   -0.39   0.00
 2   0.0281  -0.1772   -1.72  -1.41    0.51   0.00   -0.39  -0.00
 3   0.0384  -0.2426   -1.40  -1.41    0.12  -0.00   -0.39   0.00
 4   0.0355  -0.2243   -1.48  -1.41   -0.27  -0.00   -0.39  -0.00
 5   0.0222  -0.1400   -1.95  -1.41   -0.67   0.00   -0.39   0.00
 6   0.0093  -0.0590   -2.82  -1.41   -1.06  -0.00   -0.39  -0.00
 7   0.0027  -0.0168   -4.07  -1.41   -1.45   0.00   -0.39   0.00
 8   0.0005  -0.0032   -5.72  -1.41   -1.85  -0.00   -0.39  -0.00
 9   0.0001  -0.0004   -7.77  -1.41   -2.24   0.00   -0.39   0.01
10   0.0000  -0.0000  -10.20  -1.41   -2.58  -0.04   -0.31  -0.08
11   0.0000  -0.0000  -12.93  -1.49   -0.81   0.01    3.85   0.17
12   0.0000  -0.0000  -11.82  -1.39    1.88   0.04    1.53  -0.12
13   0.0000  -0.0001   -9.18  -1.42    2.46  -0.01   -0.36   0.03
14   0.0002  -0.0010   -6.90  -1.41    2.08   0.00   -0.39  -0.00
15   0.0010  -0.0065   -5.02  -1.41    1.69  -0.00   -0.39   0.00


 k         A          B         C        Freq       WDFT Mag
--     --------   -------    -------    --------    -------
 0     -0.1964     1.2959    -3.5245    3.299796    0.0295
 1     -0.1963     1.2959    -3.5245    3.300040    0.0885
 2     -0.1964     1.2959    -3.5245    3.299990    0.1794
 3     -0.1963     1.2959    -3.5245    3.300002    0.2456
 4     -0.1964     1.2959    -3.5246    3.300004    0.2271
 5     -0.1963     1.2959    -3.5245    3.299980    0.1417
 6     -0.1964     1.2960    -3.5248    3.300088    0.0597
 7     -0.1963     1.2955    -3.5232    3.299475    0.0170
 8     -0.1966     1.2990    -3.5362    3.304434    0.0033
 9     -0.1942     1.2586    -3.3652    3.240889    0.0004
10     -0.1526     0.4686     0.3770    1.535379    0.0000
11      1.9237   -43.1328   228.7651   11.211127    0.0000
12      0.7647   -16.4769    75.7836   10.773327    0.0000
13     -0.1810     7.1653   -71.7438   19.795857    0.0001
14     -0.1970     7.5981   -74.6611   19.283604    0.0010
15     -0.1963     7.5768   -74.5068   19.301514    0.0066

  parabola product =  3.1416

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

On Sun, 12 Mar 2017 22:43:13 -0700 (PDT), robert bristow-johnson
<rbj@audioimagination.com> wrote:

>On Sunday, March 12, 2017 at 2:59:12 PM UTC-4, Cedron wrote:
>> >
>> >
>> >I am also surprised by the alternating nature of the bin values.
>> 
>> Not any more.  The Fourier transform works the same way forward and
>> backward.  So, if you multiply a signal with a complex tone you get a
>> shift in the frequency domain.  The reverse is if you multiply the
>> frequency domain with a complex "tone" you get a shift in the time domain,
>> aka a phase shift.
>> 
>> Alternating positive and negative ones corresponds to a complex tone with
>> a frequency of N/2.  So here is another reason that the window should be
>> centered at N/2.
>
>it's centering the window in the time domain at N/2 that causes the alternating of the signs after the FFT.  this is why MATLAB has this function called fftshift().
>
>r b-j

fftshift just swaps the output halves of the FFT depending on whether
you want DC in the middle or not.   Many people used to call this and
how it is used in algorithms (like fast convolution) "folding" the
output or not.

Naturally you can exploit some of the symmetry properties with this,
too, like some alternating sign conditions with even or odd symmetry
signals, etc.



---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus

On Thu, 09 Mar 2017 23:58:56 +0100, Michael Plet <me@home.com> 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

Michael,

I got your email but can't respond due to "relaying not allowed",
which appears to have something to do with your mail server.

FWIW, though, yes, you have it right.   I was going to add some info
on tradeoffs, but they're not a big deal.

Cheers,

Eric

---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus

On Wed, 15 Mar 2017 16:51:16 GMT, eric.jacobsen@ieee.org wrote:

>On Thu, 09 Mar 2017 23:58:56 +0100, Michael Plet <me@home.com> 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
>
>Michael,
>
>I got your email but can't respond due to "relaying not allowed",
>which appears to have something to do with your mail server.
>
>FWIW, though, yes, you have it right.   I was going to add some info
>on tradeoffs, but they're not a big deal.
>
>Cheers,
>
>Eric
>
>---
>This email has been checked for viruses by Avast antivirus software.
>https://www.avast.com/antivirus

Thank you for your quick answer Eric.
I'm glad I got it right.
"Relaying not allowed". That is strange.

Thanks,
Michael