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
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
On Fri, 10 Mar 2017 10:44:35 -0800 (PST), dbd <d.dalrymple@sbcglobal.net> wrote:>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. DalrympleI don't understand. Why 43-44? With N=32 the mirror of 3 should be at 29. Michael
Reply by ●March 10, 20172017-03-10
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-jI 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)) Michael
Reply by ●March 10, 20172017-03-10
On Fri, 10 Mar 2017 11:40:23 -0600, "Cedron" <103185@DSPRelated> wrote:>> >>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.comI would feel proud if I can contribute with anything. Now I need to find out how I can convert a real signal to a complex signal. Regarding an analytic formula all I can say is that there is not an analytic anti derivative to a Gaussian, but with limits it may be possible to find an analytic expression for the value. Michael
Reply by ●March 10, 20172017-03-10
On Fri, 10 Mar 2017 22:00:22 +0100, Michael Plet <me@home.com> wrote:>On Fri, 10 Mar 2017 11:40:23 -0600, "Cedron" <103185@DSPRelated> >wrote: > >>> >>>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 > > >I would feel proud if I can contribute with anything. > >Now I need to find out how I can convert a real signal to a complex >signal.If there a barrier to just creating it as a complex signal initially? S(x) = e^j(x) = cos(x) + j*sin(x) Elsewise if the signal is spectrally rich, you can mix it down to baseband and filter it, but for just testing with pure tones it's easy enough to just generate the complex tone initially.>Regarding an analytic formula all I can say is that there is not an >analytic anti derivative to a Gaussian, but with limits it may be >possible to find an analytic expression for the value. > >Michael--- This email has been checked for viruses by Avast antivirus software. https://www.avast.com/antivirus
Reply by ●March 10, 20172017-03-10
On Fri, 10 Mar 2017 21:05:43 GMT, eric.jacobsen@ieee.org wrote:>On Fri, 10 Mar 2017 22:00:22 +0100, Michael Plet <me@home.com> wrote: > >>On Fri, 10 Mar 2017 11:40:23 -0600, "Cedron" <103185@DSPRelated> >>wrote: >> >>>> >>>>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 >> >> >>I would feel proud if I can contribute with anything. >> >>Now I need to find out how I can convert a real signal to a complex >>signal. > >If there a barrier to just creating it as a complex signal initially? > >S(x) = e^j(x) = cos(x) + j*sin(x) > >Elsewise if the signal is spectrally rich, you can mix it down to >baseband and filter it, but for just testing with pure tones it's easy >enough to just generate the complex tone initially. > >>Regarding an analytic formula all I can say is that there is not an >>analytic anti derivative to a Gaussian, but with limits it may be >>possible to find an analytic expression for the value. >> >>Michael > > >--- >This email has been checked for viruses by Avast antivirus software. >https://www.avast.com/antivirusThanks. I'll just create a complex signal initially. Btw, I use "i" ;-) Michael
Reply by ●March 10, 20172017-03-10
On Fri, 10 Mar 2017 22:12:34 +0100, Michael Plet <me@home.com> wrote:>On Fri, 10 Mar 2017 21:05:43 GMT, eric.jacobsen@ieee.org wrote: > >>On Fri, 10 Mar 2017 22:00:22 +0100, Michael Plet <me@home.com> wrote: >> >>>On Fri, 10 Mar 2017 11:40:23 -0600, "Cedron" <103185@DSPRelated> >>>wrote: >>> >>>>> >>>>>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 >>> >>> >>>I would feel proud if I can contribute with anything. >>> >>>Now I need to find out how I can convert a real signal to a complex >>>signal. >> >>If there a barrier to just creating it as a complex signal initially? >> >>S(x) = e^j(x) = cos(x) + j*sin(x) >> >>Elsewise if the signal is spectrally rich, you can mix it down to >>baseband and filter it, but for just testing with pure tones it's easy >>enough to just generate the complex tone initially. >> >>>Regarding an analytic formula all I can say is that there is not an >>>analytic anti derivative to a Gaussian, but with limits it may be >>>possible to find an analytic expression for the value. >>> >>>Michael >> >> >>--- >>This email has been checked for viruses by Avast antivirus software. >>https://www.avast.com/antivirus > >Thanks. I'll just create a complex signal initially. > >Btw, I use "i" ;-) > >MichaelI 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. Michael
Reply by ●March 10, 20172017-03-10
On Friday, March 10, 2017 at 12:50:16 PM UTC-8, Michael Plet wrote:> ... > > I don't understand. Why 43-44? With N=32 the mirror of 3 should be at > 29. > > MichaelMichael Here's a problem with inadequate processing descriptions. The only transform size like number I remember from your past posts is something like data size of 100. Wasn't that your fft? (If not, what was it and why did you post that number?) I went to the original post in this thread and found no other value. I suggest you better document whenever (that is, every time) you publish results. Then you can try different values and readers will know the appropriate value. You can't (successfully) expect readers to be as familiar with the context of you test cases as you are. Actually I hoped to suggest that you calculate with images at a different separation, but when I considered how to do that I thought it might also identify this problem with knowing what your processing parameters actually are. How did you expect to get meaningful answers to questions when you fail to accurately describe the situation? Dale B. Dalrymple
Reply by ●March 10, 20172017-03-10
> > >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)) > >MichaelYou may actually be introducing another approximation here. Consider this: ln( a + i*b ) = ln( r * e^( i * theta ) ) = ln( r ) + ( i * theta ) = ln( a^2 + b^2 ) / 2 + i * atan2( b, a ) Now the real part is the same as the log of the magnitude. If the DFT doesn't have all zero imaginary bins then they are still going to be very very small, so the error introduced is also going to be very very small. If all your bins have the same theta, then the peak formula cancels them out nicely and this is not an issue. Ced --------------------------------------- Posted through http://www.DSPRelated.com
Reply by ●March 10, 20172017-03-10
>On Friday, March 10, 2017 at 12:50:16 PM UTC-8, Michael Plet wrote: >> ... >> >> I don't understand. Why 43-44? With N=32 the mirror of 3 should be at >> 29. >> >> Michael > >Michael > >Here's a problem with inadequate processing descriptions. The only >transform size like number I remember from your past posts is somethinglike data>size of 100. Wasn't that your fft? (If not, what was it and why did you >post that number?) I went to the original post in this thread and foundno>other value. I suggest you better document whenever (that is, every time)you>publish results. Then you can try different values and readers will know >the appropriate value. You can't (successfully) expect readers to be as >familiar with the context of you test cases as you are. > >Actually I hoped to suggest that you calculate with images at adifferent>separation, but when I considered how to do that I thought it might also >identify this problem with knowing what your processing parametersactually>are. How did you expect to get meaningful answers to questions when youfail>to accurately describe the situation? > >Dale B. DalrympleSigh. The original post, and several others, clearly state N=32. In one post he mentions testing different N sizes. 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, Ced --------------------------------------- Posted through http://www.DSPRelated.com
