> > >I just did some tests with complex signals. >In that case you can get the exact value if "a" is chosen somewhere >between 0.089 and 0.129. >In the complex case "a" seems to depend on N in a different way. > >MichaelRemember, you are dealing with the noiseless case (as far as your numerical precision allows). These "a" values correspond to much narrower bell curves. I recommend that your read the following section in Julien's comparison paper for some perspective on what happens in the presence of noise. "4.3.5 Choice of the gaussian window shape factor" http://www.tsdconseil.fr/log/scriptscilab/festim/index-en.html Ced --------------------------------------- Posted through http://www.DSPRelated.com
Question about Gaussian window and parabolic interpolation
Started by ●March 9, 2017
Reply by ●March 10, 20172017-03-10
Reply by ●March 10, 20172017-03-10
Michael Modern general purpose dsp libraries support direct computation of FFTs in sizes that are products of powers of factors from 2 to 7 (as well as larger prime and composite factors). Examples are MATLAB and the GNU scientific library. Modern personal computers and cell phones can readily calculate DFTs in the hundreds if high real-time sample frequencies are not required. On Friday, March 10, 2017 at 3:19:23 PM UTC-8, Cedron wrote:> ... > Sigh > ... > A FFT cannot have 100 points. A FFT is a DFT with a size of either a > power of two and there is a different version involving prime numbers. > > Your friendly neighborhood stick puppet, > > CedThis isn't one of those cases of engineers and mathematicians thinking differently. This is someone incompetent at signal processing posing as a sock puppet to troll newsgroups. There have been good reasons for readers of comp.dsp to place no confidence in the statements of this troll about the values of algorithms for signal processing. Dale B. Dalrymple
Reply by ●March 11, 20172017-03-11
On Friday, March 10, 2017 at 3:53:46 PM UTC-5, Michael Plet wrote:> On Fri, 10 Mar 2017 11:44:25 -0800 (PST), robert bristow-johnson > <rbj@audioimagination.com> wrote: > > >On Friday, March 10, 2017 at 5:10:26 AM UTC-5, Michael Plet wrote: > >> 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, what maths are you using to get the peak estimate from the log-magnitude FFT data? are you just interpolating with the discrete peak and the two neighboring FFT points to the left and right of the discrete peak? if you do as suggested in the Mohonk paper, you want to do a discrete-derivative, which will turn your parabola into points that appear like a straight line. and then fit a straight line through those points. you might do even better in fitting to the data. > > > >r b-j > > > I am calculating the log of three consecutive bin magnitudes and then > calculating the peak location like this > > x="Index of max bin" - (y2-y0)/(2*(y2-2*y1+y0)) >yeah, that's the max point and its two neighboring points. it's okay as far as it goes, but you might get even better results fitting the quadratic to more points than three. (or, as in the 2001 paper fitting a line to the finite differences.) r b-j
Reply by ●March 11, 20172017-03-11
>> >> I am calculating the log of three consecutive bin magnitudes and then >> calculating the peak location like this >> >> x="Index of max bin" - (y2-y0)/(2*(y2-2*y1+y0)) >> > >yeah, that's the max point and its two neighboring points. it's okay as >far as it goes, but you might get even better results fitting thequadratic>to more points than three. (or, as in the 2001 paper fitting a line tothe>finite differences.) > > >r b-jI don't think he will, but I wanted to do a little coding to make sure before I said so. From my other work on frequency estimators I have found that when there is noise present bringing in more bins ends up lowering your overall effective SNR. This is because the peak bin(s) have the most signal and the further out you go the less signal you get with the same (nominally) amount of noise. So the bin SNR drops rapidly as you move away from the peak. In the noiseless case, ideally you could fit a line to two difference points, or as many as you wanted, and get the same results as Michael's method. However, because this is the finite Discrete case, your difference points aren't going to be exactly linear because there is some distortion in the parabola found in the log of the DFT bin values. The further from the peak you go, the more distortion there is. Therefore doing your calculation as tight to the peak bin as you can produces the most distortion results. The distortion gets worse with too large or too small "a" values. The point I made earlier about using the log of the bin value vs using the log of the bin magnitude turns out to be irrelevant. The first difference seems to have a constant "twist" and the second difference washes this out. Therefore Michael's quotient would have a little bit of an imaginary part to it, which you discard anyway, leaving the exact same result as using the magnitude. I've run this program for many different parameters and have selected one run that shows what I am trying to say the best. There is an inverse relationship between the width of the input parabola (in log space) and the width of output parabola similar to Equation (1) in the Mohonk paper. A_time * A_freq = (Pi/N)^2 This relationship is fairly accurate for reasonable values of "a", meaning about .10 to .25. Ced This is long, please snip if you reply. ================================ Calculation of Differences L[k] = log( Z[k] ) d = First Diff = ( L[k+1] - L[k-1] ) / 2 s = Second Diff = L[k+1] - 2 * L[k] + L[k-1] =========================================== Calculation of Parabola Parameters Time domain: theCenter = ( N - 1 ) / 2 A_time = -1 / ( 2 * a^2 * theCenter^2 ); Frequency domain: y = Ax^2 + Bx + C d = ( y(x+1) - y(x-1) ) / 2 = ( ( Ax^2 + 2Ax + A + Bx + B + C ) - ( Ax^2 - 2Ax + A + Bx - B + C ) ) / 2 = 2Ax + B s = y(x+1) - 2y(x) + y(x-1) = ( Ax^2 + 2Ax + A + Bx + B + C ) - 2 ( Ax^2 + Bx + C ) + ( Ax^2 - 2Ax + A + Bx - B + C ) = 2A A_freq = s / 2 B_freq = d - 2A_freq * n C_freq = y - A_freq * n^2 - B_freq * n Freq = -B_freq / ( 2 * A_freq ) Combined: parabola product = N * sqrt( A_time * A_freq ) =========================================== Sample Run Frequency = 6.1 N = 16 Michael's a = 0.2000 Sigma Range is -5.0000 to 5.0000 Julien's Sigma = N / 1.6000 A_time = -0.2222 n Window Signal Windowed Signal -- ------- ----------------- ----------------- 0 0.0000 0.5653 0.8249 0.0000 0.0000 1 0.0001 -0.9750 -0.2220 -0.0001 -0.0000 2 0.0012 0.8667 -0.4988 0.0010 -0.0006 3 0.0111 -0.2978 0.9546 -0.0033 0.0106 4 0.0657 -0.4293 -0.9032 -0.0282 -0.0594 5 0.2494 0.9283 0.3718 0.2315 0.0927 6 0.6065 -0.9341 0.3571 -0.5665 0.2166 7 0.9460 0.4435 -0.8963 0.4195 -0.8478 8 0.9460 0.2827 0.9592 0.2674 0.9074 9 0.6065 -0.8587 -0.5124 -0.5208 -0.3108 10 0.2494 0.9784 -0.2066 0.2440 -0.0515 11 0.0657 -0.5782 0.8159 -0.0380 0.0536 12 0.0111 -0.1292 -0.9916 -0.0014 -0.0110 13 0.0012 0.7680 0.6405 0.0009 0.0008 14 0.0001 -0.9987 0.0510 -0.0001 0.0000 15 0.0000 0.6988 -0.7154 0.0000 -0.0000 k Windowed DFT Log of WDFT First Diff Second Diff -- ---------------- ------------- ------------- ------------- 0 0.0004 0.0000 -7.90 0.09 2.12 -2.95 -0.35 -0.00 1 -0.0025 -0.0007 -5.96 -2.86 1.77 0.20 -0.35 6.28 2 0.0113 0.0059 -4.36 0.48 1.42 0.20 -0.35 -6.28 3 -0.0346 -0.0277 -3.12 -2.47 1.08 0.20 -0.35 6.28 4 0.0704 0.0837 -2.21 0.87 0.73 0.20 -0.35 -6.28 5 -0.0918 -0.1669 -1.66 -2.07 0.38 0.20 -0.35 6.28 6 0.0707 0.2237 -1.45 1.26 0.03 0.20 -0.35 -6.28 7 -0.0224 -0.2030 -1.59 -1.68 -0.31 0.20 -0.35 6.28 8 -0.0108 0.1252 -2.07 1.66 -0.66 0.20 -0.35 -6.28 9 0.0152 -0.0525 -2.91 -1.29 -1.01 0.20 -0.35 6.28 10 -0.0077 0.0149 -4.09 2.05 -1.35 0.20 -0.35 -6.28 11 0.0023 -0.0028 -5.61 -0.90 -1.70 0.20 -0.35 6.28 12 -0.0004 0.0004 -7.49 2.44 -2.05 0.20 -0.35 -6.28 13 0.0001 -0.0000 -9.71 -0.50 -2.82 0.20 -1.20 6.28 14 -0.0000 0.0000 -13.13 2.84 -0.25 1.77 6.35 -3.14 15 -0.0000 0.0000 -10.20 3.03 2.61 -1.37 -0.63 -3.14 k A B C Freq DFT Mag -- -------- ------- ------- -------- ------- 0 -0.1769 2.1200 -7.9037 5.992783 0.0004 1 -0.1735 2.1166 -7.9037 6.099591 0.0026 2 -0.1735 2.1166 -7.9037 6.100024 0.0127 3 -0.1735 2.1166 -7.9037 6.099993 0.0444 4 -0.1735 2.1166 -7.9037 6.100002 0.1093 5 -0.1735 2.1166 -7.9037 6.099999 0.1905 6 -0.1735 2.1166 -7.9037 6.100000 0.2346 7 -0.1735 2.1166 -7.9037 6.100001 0.2042 8 -0.1735 2.1166 -7.9037 6.099997 0.1256 9 -0.1735 2.1166 -7.9038 6.100010 0.0546 10 -0.1735 2.1165 -7.9035 6.099953 0.0168 11 -0.1735 2.1169 -7.9054 6.100395 0.0036 12 -0.1745 2.1404 -8.0401 6.131973 0.0006 13 -0.5991 12.7544 -74.2717 10.644875 0.0001 14 3.1757 -89.1646 612.7377 14.038612 0.0000 15 -0.3168 12.1174 -120.6839 19.125181 0.0000 parabola product = 3.1416 --------------------------------------- Posted through http://www.DSPRelated.com
Reply by ●March 11, 20172017-03-11
Here is another run with a broader window function and thus a narrower peak: Michael's a = 0.3000 Sigma Range is -3.3333 to 3.3333 k A B C Freq DFT Mag -- -------- ------- ------- -------- ------- 0 1.5163 0.1921 -12.3566 -0.063356 0.0000 1 0.6392 1.0692 -12.3566 -0.836309 0.0000 2 -0.0570 3.1579 -13.7490 27.710585 0.0005 3 -0.4251 4.9985 -15.9578 5.879191 0.0083 4 -0.3860 4.7251 -15.4891 6.119904 0.0630 5 -0.3915 4.7744 -15.5986 6.097238 0.2199 6 -0.3896 4.7531 -15.5405 6.100209 0.3510 7 -0.3916 4.7794 -15.6255 6.102276 0.2571 8 -0.3852 4.6829 -15.2653 6.078934 0.0860 9 -0.4437 5.6778 -19.4787 6.398249 0.0133 10 0.2803 -8.0779 45.6798 14.410109 0.0008 11 0.7801 -18.5735 100.6567 11.904927 0.0001 12 0.1817 -4.8107 21.6699 13.238517 0.0001 13 -0.0387 0.6994 -12.7128 9.033657 0.0000 14 -0.0890 2.0581 -21.8717 11.558128 0.0000 15 -0.4002 11.0811 -87.2105 13.845464 0.0000 parabola product = 3.1385 --------------------------------------- Posted through http://www.DSPRelated.com
Reply by ●March 12, 20172017-03-12
The window function as Michael used it is "symmetric", that is centered at (N-1)/2. With regard to the fact that there are two versions of the Von Hann window, I thought I would try the "periodic" version of the Gaussian, that is centered at N/2. The results were much better. So I did some web searching and it seems that the only the "symmetric" version is used anywhere. This includes MatLab (I am not a user) which offers both versions of the Von Hann, but only the symmetric version of the Gaussian. Here are the results. Notice that the imaginary part of the First Diff column is now all zeroes in the region of interest. This will make the use of this window much more accurate when its use gets expanded to include the frequency sweep and amplitude ramp functionality as described in r b-j's Mohonk paper. Currently the A, B, C values are calculated solely on the real part of the Difference values. Ced ============================================================= N = 16 Michael's a = 0.2000 Sigma Range is -5.0000 to 4.6875 Julien's Sigma = N / 1.6000 The parabola A = -0.1953 n Window Signal Windowed Signal -- ------- ----------------- ----------------- 0 0.0000 0.5653 0.8249 0.0000 0.0000 1 0.0001 -0.9750 -0.2220 -0.0001 -0.0000 2 0.0009 0.8667 -0.4988 0.0008 -0.0004 3 0.0076 -0.2978 0.9546 -0.0023 0.0072 4 0.0439 -0.4293 -0.9032 -0.0189 -0.0397 5 0.1724 0.9283 0.3718 0.1601 0.0641 6 0.4578 -0.9341 0.3571 -0.4276 0.1635 7 0.8226 0.4435 -0.8963 0.3648 -0.7373 8 1.0000 0.2827 0.9592 0.2827 0.9592 9 0.8226 -0.8587 -0.5124 -0.7064 -0.4215 10 0.4578 0.9784 -0.2066 0.4480 -0.0946 11 0.1724 -0.5782 0.8159 -0.0997 0.1407 12 0.0439 -0.1292 -0.9916 -0.0057 -0.0436 13 0.0076 0.7680 0.6405 0.0058 0.0049 14 0.0009 -0.9987 0.0510 -0.0009 0.0000 15 0.0001 0.6988 -0.7154 0.0000 -0.0000 k Windowed DFT Log of WDFT First Diff Second Diff -- ---------------- ------------- ------------- ------------- 0 0.0000 0.0002 -8.73 1.28 2.40 -0.00 -0.37 -6.28 1 -0.0004 -0.0014 -6.52 -1.86 2.01 0.00 -0.40 6.28 2 0.0026 0.0087 -4.70 1.28 1.62 -0.00 -0.39 -6.28 3 -0.0106 -0.0361 -3.28 -1.86 1.22 0.00 -0.39 6.28 4 0.0297 0.1007 -2.25 1.28 0.83 -0.00 -0.39 -6.28 5 -0.0558 -0.1894 -1.62 -1.86 0.43 0.00 -0.39 6.28 6 0.0707 0.2400 -1.39 1.28 0.04 -0.00 -0.39 -6.28 7 -0.0604 -0.2049 -1.54 -1.86 -0.36 -0.00 -0.39 6.28 8 0.0348 0.1179 -2.10 1.28 -0.75 0.00 -0.39 -6.28 9 -0.0135 -0.0457 -3.04 -1.86 -1.14 -0.00 -0.39 6.28 10 0.0035 0.0119 -4.39 1.28 -1.54 0.00 -0.39 -6.28 11 -0.0006 -0.0021 -6.12 -1.86 -1.93 -0.00 -0.40 6.28 12 0.0001 0.0002 -8.26 1.28 -2.32 0.00 -0.38 -6.28 13 -0.0000 -0.0000 -10.77 -1.85 -2.57 -0.02 -0.12 6.23 14 0.0000 0.0000 -13.40 1.24 -0.27 0.00 4.72 -6.18 15 -0.0000 -0.0000 -11.31 -1.85 2.34 0.02 0.49 6.22 k A B C Freq DFT Mag -- -------- ------- ------- -------- ------- 0 -0.1859 2.3981 -8.7299 6.450833 0.0002 1 -0.1982 2.4104 -8.7299 6.081274 0.0015 2 -0.1973 2.4077 -8.7281 6.101858 0.0091 3 -0.1974 2.4083 -8.7288 6.099731 0.0376 4 -0.1974 2.4081 -8.7285 6.100055 0.1050 5 -0.1974 2.4082 -8.7287 6.099986 0.1974 6 -0.1974 2.4082 -8.7285 6.100001 0.2502 7 -0.1974 2.4082 -8.7287 6.100010 0.2136 8 -0.1974 2.4081 -8.7283 6.099957 0.1229 9 -0.1974 2.4084 -8.7297 6.100199 0.0477 10 -0.1973 2.4068 -8.7220 6.098719 0.0125 11 -0.1979 2.4192 -8.7871 6.111886 0.0022 12 -0.1910 2.2609 -7.8785 5.917735 0.0003 13 -0.0591 -1.0383 12.7084 -8.789554 0.0000 14 2.3613 -66.3876 453.2111 14.057524 0.0000 15 0.2472 -5.0781 9.2462 10.272795 0.0000 parabola product = 3.1416 --------------------------------------- Posted through http://www.DSPRelated.com
Reply by ●March 12, 20172017-03-12
>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 N2, 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. > >MichaelI finally have an answer for you, sort of. It seems that the best value for "a" is the one where the size of the Gaussian in the time domain matches the size in the frequency domain. I have found the using "N" instead of "N-1" in the window function produces better results. This answer is also based on a complex signal, not a real one. Why you got such consistent results with different Ns on a real signal is still an open question. The formula is a = sqrt( 2 / (N*Pi) ) For N=32 the value is 0.141047396... Ced --------------------------------------- Posted through http://www.DSPRelated.com
Reply by ●March 12, 20172017-03-12
On Sun, 12 Mar 2017 08:28:59 -0500, "Cedron" <103185@DSPRelated> 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 N2, 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 finally have an answer for you, sort of. > >It seems that the best value for "a" is the one where the size of the >Gaussian in the time domain matches the size in the frequency domain. > >I have found the using "N" instead of "N-1" in the window function >produces better results. > >This answer is also based on a complex signal, not a real one. Why you >got such consistent results with different Ns on a real signal is still an >open question. > >The formula is a = sqrt( 2 / (N*Pi) ) > >For N=32 the value is 0.141047396... > >Ced >--------------------------------------- >Posted through http://www.DSPRelated.comThanks. ( 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
Reply by ●March 12, 20172017-03-12
On Sun, 12 Mar 2017 15:08:33 +0100, Michael Plet <me@home.com> wrote:>On Sun, 12 Mar 2017 08:28:59 -0500, "Cedron" <103185@DSPRelated> >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 N2, 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 finally have an answer for you, sort of. >> >>It seems that the best value for "a" is the one where the size of the >>Gaussian in the time domain matches the size in the frequency domain. >> >>I have found the using "N" instead of "N-1" in the window function >>produces better results. >> >>This answer is also based on a complex signal, not a real one. Why you >>got such consistent results with different Ns on a real signal is still an >>open question. >> >>The formula is a = sqrt( 2 / (N*Pi) ) >> >>For N=32 the value is 0.141047396... >> >>Ced >>--------------------------------------- >>Posted through http://www.DSPRelated.com > >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. > >MichaelI 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
Reply by ●March 12, 20172017-03-12
> >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. > >MichaelYou'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
