Hi I'm using fftw3 on ubuntu linux in a convolution of 1-dimensional data.
for ( int k=0 ; k<=2*N+1 ; k++ ){
fyz[k][0] = (fy[k][0] * fz[k][0]) - (fy[k][1] * fz[k][1]) ;
fyz[k][1] = (fy[k][0] * fz[k][1]) + (fy[k][1] * fz[k][0]) ;
}
a sample of the last output is
fyz[131068] 8419794.642330 3060741.589590
fyz[131069] -8637597.560357 -3493977.608001
fyz[131070] 8298749.249280 3911285.476558
fyz[131071] -7704225.570165 -3227891.546761
fyz[131072] 8505818.548396 2791756.773350
I run this on Matlab and I get
Column 131068
8436336.58969844 + 3059158.22027252i
Column 131069
-8631164.17521021 - 3505892.14995256i
Column 131070
8278270.82666789 + 3902196.48863977i
Column 131071
-7705102.14197655 - 3203236.6220612i
Column 131072
8524950.31393116 + 2790538.23412558i
spot the obvious slight differences ioj results.
the data on ubuntu and matlab matches going into the convolution, but not
on the other side. so must be roundoff error? how can I curtail this to
make results the same?
also eps is the same on my ubuntu machine as it is on my matlab machine.
i need this accurate because the inverse is subsequently way different
between the two machines.
thanks
Ian
fftw3 roundoff error in convolution c++
Started by ●September 17, 2009
Reply by ●September 17, 20092009-09-17
"Kwisatch" <ibellinfantie@yahoo.co.uk> writes:> Hi I'm using fftw3 on ubuntu linux in a convolution of 1-dimensional data. > > for ( int k=0 ; k<=2*N+1 ; k++ ){ > fyz[k][0] = (fy[k][0] * fz[k][0]) - (fy[k][1] * fz[k][1]) ; > fyz[k][1] = (fy[k][0] * fz[k][1]) + (fy[k][1] * fz[k][0]) ; > } > > a sample of the last output is > > fyz[131068] 8419794.642330 3060741.589590 > fyz[131069] -8637597.560357 -3493977.608001 > fyz[131070] 8298749.249280 3911285.476558 > fyz[131071] -7704225.570165 -3227891.546761 > fyz[131072] 8505818.548396 2791756.773350 > > > I run this on Matlab and I get > > Column 131068 > > 8436336.58969844 + 3059158.22027252i > > Column 131069 > > -8631164.17521021 - 3505892.14995256i > > Column 131070 > > 8278270.82666789 + 3902196.48863977i > > Column 131071 > > -7705102.14197655 - 3203236.6220612i > > Column 131072 > > 8524950.31393116 + 2790538.23412558i > > spot the obvious slight differences ioj results. > > the data on ubuntu and matlab matches going into the convolution, but not > on the other side. so must be roundoff error? how can I curtail this to > make results the same? > also eps is the same on my ubuntu machine as it is on my matlab machine. > > i need this accurate because the inverse is subsequently way different > between the two machines. > > thanks > > IanIan, how big is your FFT? (What is "N"?) This seems way off - not just roundoff error - but I haven't analyzed it closely. I know Matlab uses 8-byte floats (double floats), and I think you can configure FFTW to use either single- or double-precision floats - perhaps that's you're discrepancy? -- Randy Yates % "Maybe one day I'll feel her cold embrace, Digital Signal Labs % and kiss her interface, mailto://yates@ieee.org % til then, I'll leave her alone." http://www.digitalsignallabs.com % 'Yours Truly, 2095', *Time*, ELO
Reply by ●September 18, 20092009-09-18
>"Kwisatch" <ibellinfantie@yahoo.co.uk> writes: > >> Hi I'm using fftw3 on ubuntu linux in a convolution of 1-dimensionaldata.>> >> for ( int k=0 ; k<=2*N+1 ; k++ ){ >> fyz[k][0] = (fy[k][0] * fz[k][0]) - (fy[k][1] * fz[k][1]) ; >> fyz[k][1] = (fy[k][0] * fz[k][1]) + (fy[k][1] * fz[k][0]) ; >> } >> >> a sample of the last output is >> >> fyz[131068] 8419794.642330 3060741.589590 >> fyz[131069] -8637597.560357 -3493977.608001 >> fyz[131070] 8298749.249280 3911285.476558 >> fyz[131071] -7704225.570165 -3227891.546761 >> fyz[131072] 8505818.548396 2791756.773350 >> >> >> I run this on Matlab and I get >> >> Column 131068 >> >> 8436336.58969844 + 3059158.22027252i >> >> Column 131069 >> >> -8631164.17521021 - 3505892.14995256i >> >> Column 131070 >> >> 8278270.82666789 + 3902196.48863977i >> >> Column 131071 >> >> -7705102.14197655 - 3203236.6220612i >> >> Column 131072 >> >> 8524950.31393116 + 2790538.23412558i >> >> spot the obvious slight differences ioj results. >> >> the data on ubuntu and matlab matches going into the convolution, butnot>> on the other side. so must be roundoff error? how can I curtail thisto>> make results the same? >> also eps is the same on my ubuntu machine as it is on my matlabmachine.>> >> i need this accurate because the inverse is subsequently way different >> between the two machines. >> >> thanks >> >> Ian > >Ian, how big is your FFT? (What is "N"?) > >This seems way off - not just roundoff error - but I haven't analyzed it >closely. > >I know Matlab uses 8-byte floats (double floats), and I think you can >configure FFTW to use either single- or double-precision floats - >perhaps that's you're discrepancy? >-- >Randy Yates % "Maybe one day I'll feel her coldembrace,>Digital Signal Labs % and kiss herinterface,>mailto://yates@ieee.org % til then, I'll leave heralone.">http://www.digitalsignallabs.com % 'Yours Truly, 2095', *Time*,ELO>Hi Randy for the transform N=2*(2^16) i checked this again, the loss of accuracy occurs from the complex arithmetic BEFORE the fftw. the examples below are at the end of the 131072 samples. Ubuntu or C++ rounds the real part to 1 zz[131070] (1,-5.23869e-06) zz[131071] (1,-2.32831e-06) zz[131072] (1,-5.82077e-07) whereas in matlab they are Column 131070 0.999999999986278 - 5.2386894821881e-006i Column 131071 0.99999999999729 - 2.32830643653659e-006i Column 131072 0.999999999999831 - 5.82076609134641e-007i after fftw3 accuracy is out, and this loss progresses through to the inverse fftw. ubuntu fz[131070] 1675.330280 -1647.737086 fz[131071] 1633.353689 -1616.370818 fz[131072] 1662.342630 -1660.870246 matlab Column 131070 1674.30929408701 - 1647.99273838256i Column 131071 1632.17793684431 - 1616.34051114058i Column 131072 1661.40285267096 - 1661.02707675479i so the problem must be in ubuntu or C++ rounding? unless you have any other suggestions, I guess I should go back to looking at ubuntu and c++. Thanks Ian
Reply by ●September 18, 20092009-09-18
"Kwisatch" <ibellinfantie@yahoo.co.uk> writes:> [...] > Hi Randy > > for the transform > N=2*(2^16)That would be 2^17. That is a HUGE FFT! 128K?> i checked this again, the loss of accuracy occurs from the complex > arithmetic BEFORE the fftw.Then it has nothing to do with the FFT. Unless you show us these operations, I'm not sure how we can help you. -- Randy Yates % "Bird, on the wing, Digital Signal Labs % goes floating by mailto://yates@ieee.org % but there's a teardrop in his eye..." http://www.digitalsignallabs.com % 'One Summer Dream', *Face The Music*, ELO
Reply by ●September 18, 20092009-09-18
>>"Kwisatch" <ibellinfantie@yahoo.co.uk> writes: >> >>> Hi I'm using fftw3 on ubuntu linux in a convolution of 1-dimensional >data. >>> >>> for ( int k=0 ; k<=2*N+1 ; k++ ){ >>> fyz[k][0] = (fy[k][0] * fz[k][0]) - (fy[k][1] * fz[k][1]) ; >>> fyz[k][1] = (fy[k][0] * fz[k][1]) + (fy[k][1] * fz[k][0]) ; >>> } >>> >>> a sample of the last output is >>> >>> fyz[131068] 8419794.642330 3060741.589590 >>> fyz[131069] -8637597.560357 -3493977.608001 >>> fyz[131070] 8298749.249280 3911285.476558 >>> fyz[131071] -7704225.570165 -3227891.546761 >>> fyz[131072] 8505818.548396 2791756.773350 >>> >>> >>> I run this on Matlab and I get >>> >>> Column 131068 >>> >>> 8436336.58969844 + 3059158.22027252i >>> >>> Column 131069 >>> >>> -8631164.17521021 - 3505892.14995256i >>> >>> Column 131070 >>> >>> 8278270.82666789 + 3902196.48863977i >>> >>> Column 131071 >>> >>> -7705102.14197655 - 3203236.6220612i >>> >>> Column 131072 >>> >>> 8524950.31393116 + 2790538.23412558i >>> >>> spot the obvious slight differences ioj results. >>> >>> the data on ubuntu and matlab matches going into the convolution, but >not >>> on the other side. so must be roundoff error? how can I curtail this >to >>> make results the same? >>> also eps is the same on my ubuntu machine as it is on my matlab >machine. >>> >>> i need this accurate because the inverse is subsequently waydifferent>>> between the two machines. >>> >>> thanks >>> >>> Ian >> >>Ian, how big is your FFT? (What is "N"?) >> >>This seems way off - not just roundoff error - but I haven't analyzedit>>closely. >> >>I know Matlab uses 8-byte floats (double floats), and I think you can >>configure FFTW to use either single- or double-precision floats - >>perhaps that's you're discrepancy? >>-- >>Randy Yates % "Maybe one day I'll feel her cold >embrace, >>Digital Signal Labs % and kiss her >interface, >>mailto://yates@ieee.org % til then, I'll leave her >alone." >>http://www.digitalsignallabs.com % 'Yours Truly, 2095', *Time*, >ELO >> >Hi Randy > >for the transform >N=2*(2^16) > >i checked this again, the loss of accuracy occurs from the complex >arithmetic BEFORE the fftw. the examples below are at the end of the131072>samples. Ubuntu or C++ rounds the real part to 1 > >zz[131070] (1,-5.23869e-06) >zz[131071] (1,-2.32831e-06) >zz[131072] (1,-5.82077e-07) > >whereas in matlab they are > > Column 131070 > > 0.999999999986278 - 5.2386894821881e-006i > > Column 131071 > > 0.99999999999729 - 2.32830643653659e-006i > > Column 131072 > > 0.999999999999831 - 5.82076609134641e-007i > >after fftw3 accuracy is out, and this loss progresses through to the >inverse fftw. > >ubuntu >fz[131070] 1675.330280 -1647.737086 >fz[131071] 1633.353689 -1616.370818 >fz[131072] 1662.342630 -1660.870246 > >matlab > Column 131070 > > 1674.30929408701 - 1647.99273838256i > > Column 131071 > > 1632.17793684431 - 1616.34051114058i > > Column 131072 > > 1661.40285267096 - 1661.02707675479i > >so the problem must be in ubuntu or C++ rounding? unless you have any >other suggestions, I guess I should go back to looking at ubuntu andc++.> >Thanks > >Ian >That is, I need greater accuracy in C++/ubuntu so thay my values do not go to 1 too soon. zz[131026] (0.999999,-0.00128581) zz[131027] (0.999999,-0.00123167) zz[131028] (0.999999,-0.0011787) zz[131029] (0.999999,-0.0011269) zz[131030] (0.999999,-0.00107626) zz[131031] (0.999999,-0.00102678) zz[131032] (1,-0.000978471) zz[131033] (1,-0.000931322) zz[131034] (1,-0.000885338) zz[131035] (1,-0.000840519) zz[131036] (1,-0.000796863)
Reply by ●September 18, 20092009-09-18
Randy Yates wrote:> "Kwisatch" <ibellinfantie@yahoo.co.uk> writes: >> [...] >> Hi Randy >> >> for the transform >> N=2*(2^16) > > That would be 2^17. That is a HUGE FFT! 128K? >A little while ago I knocked up a simple FFT convolver program, initially for experimenting with diy FIR impulse responses, but inevitably also as a first foray into convolutional reverb. For offline processes the otherwise extreme latency does not matter. Nothing fancy such as partitioned convolution (the next project). One of our CDP users asked if he could beta-test it for a project he was working on. He took a long track (I forget whether it was the Dam Busters march or the Thunderbirds theme), but about 30secs anyway at sr=44100, and convolved it with itself. He reported himself very pleased with the results. So I guess my code worked - clean sound - no obvious distortion or other nasties. You can get some amazing complex but coherent quasi-granular textures that way. Of course the FFT has to be twice the IR length, so that was an FFT the next power-of-two above 2.646MB. There are even more left-field computer musicians who like to experiment with FFTing whole songs/tracks/CDs and then hacking the bins in various semi-dsp-ish ways that would make the hair of some folk here turn white. Years ago Notam provided a tool "mammut" for just this job: http://archive.notam02.no/arkiv/doc/mammut All of which is a long-winded way of saying that for some computer musicians, a 128K FFT is small change. Do need double-precision though. :-) Richard Dobson
Reply by ●September 18, 20092009-09-18
Richard Dobson wrote:> Randy Yates wrote: > >> "Kwisatch" <ibellinfantie@yahoo.co.uk> writes: >> >>> [...] >>> Hi Randy >>> >>> for the transform >>> N=2*(2^16) >> >> >> That would be 2^17. That is a HUGE FFT! 128K? > > All of which is a long-winded way of saying that for some computer > musicians, a 128K FFT is small change.I don't care what those shameful queers do, but there are the real problems where you have to deal with the events of the very different time scales. An example of such problem is an analysis of a Class D amplifier. The FFT of the size of 128M samples is not unusual there.> Do need double-precision though.Certainly. Vladimir Vassilevsky DSP and Mixed Signal Design Consultant http://www.abvolt.com
Reply by ●September 18, 20092009-09-18
>>>"Kwisatch" <ibellinfantie@yahoo.co.uk> writes: >>> >>>> Hi I'm using fftw3 on ubuntu linux in a convolution of 1-dimensional >>data. >>>> >>>> for ( int k=0 ; k<=2*N+1 ; k++ ){ >>>> fyz[k][0] = (fy[k][0] * fz[k][0]) - (fy[k][1] * fz[k][1]) ; >>>> fyz[k][1] = (fy[k][0] * fz[k][1]) + (fy[k][1] * fz[k][0]) ; >>>> } >>>> >>>> a sample of the last output is >>>> >>>> fyz[131068] 8419794.642330 3060741.589590 >>>> fyz[131069] -8637597.560357 -3493977.608001 >>>> fyz[131070] 8298749.249280 3911285.476558 >>>> fyz[131071] -7704225.570165 -3227891.546761 >>>> fyz[131072] 8505818.548396 2791756.773350 >>>> >>>> >>>> I run this on Matlab and I get >>>> >>>> Column 131068 >>>> >>>> 8436336.58969844 + 3059158.22027252i >>>> >>>> Column 131069 >>>> >>>> -8631164.17521021 - 3505892.14995256i >>>> >>>> Column 131070 >>>> >>>> 8278270.82666789 + 3902196.48863977i >>>> >>>> Column 131071 >>>> >>>> -7705102.14197655 - 3203236.6220612i >>>> >>>> Column 131072 >>>> >>>> 8524950.31393116 + 2790538.23412558i >>>> >>>> spot the obvious slight differences ioj results. >>>> >>>> the data on ubuntu and matlab matches going into the convolution,but>>not >>>> on the other side. so must be roundoff error? how can I curtail this >>to >>>> make results the same? >>>> also eps is the same on my ubuntu machine as it is on my matlab >>machine. >>>> >>>> i need this accurate because the inverse is subsequently way >different >>>> between the two machines. >>>> >>>> thanks >>>> >>>> Ian >>> >>>Ian, how big is your FFT? (What is "N"?) >>> >>>This seems way off - not just roundoff error - but I haven't analyzed >it >>>closely. >>> >>>I know Matlab uses 8-byte floats (double floats), and I think you can >>>configure FFTW to use either single- or double-precision floats - >>>perhaps that's you're discrepancy? >>>-- >>>Randy Yates % "Maybe one day I'll feel her cold >>embrace, >>>Digital Signal Labs % and kiss her >>interface, >>>mailto://yates@ieee.org % til then, I'll leave her >>alone." >>>http://www.digitalsignallabs.com % 'Yours Truly, 2095', *Time*, >>ELO >>> >>Hi Randy >> >>for the transform >>N=2*(2^16) >> >>i checked this again, the loss of accuracy occurs from the complex >>arithmetic BEFORE the fftw. the examples below are at the end of the >131072 >>samples. Ubuntu or C++ rounds the real part to 1 >> >>zz[131070] (1,-5.23869e-06) >>zz[131071] (1,-2.32831e-06) >>zz[131072] (1,-5.82077e-07) >> >>whereas in matlab they are >> >> Column 131070 >> >> 0.999999999986278 - 5.2386894821881e-006i >> >> Column 131071 >> >> 0.99999999999729 - 2.32830643653659e-006i >> >> Column 131072 >> >> 0.999999999999831 - 5.82076609134641e-007i >> >>after fftw3 accuracy is out, and this loss progresses through to the >>inverse fftw. >> >>ubuntu >>fz[131070] 1675.330280 -1647.737086 >>fz[131071] 1633.353689 -1616.370818 >>fz[131072] 1662.342630 -1660.870246 >> >>matlab >> Column 131070 >> >> 1674.30929408701 - 1647.99273838256i >> >> Column 131071 >> >> 1632.17793684431 - 1616.34051114058i >> >> Column 131072 >> >> 1661.40285267096 - 1661.02707675479i >> >>so the problem must be in ubuntu or C++ rounding? unless you have any >>other suggestions, I guess I should go back to looking at ubuntu and >c++. >> >>Thanks >> >>Ian >> > >That is, I need greater accuracy in C++/ubuntu so thay my values do notgo>to 1 too soon. > >zz[131026] (0.999999,-0.00128581) >zz[131027] (0.999999,-0.00123167) >zz[131028] (0.999999,-0.0011787) >zz[131029] (0.999999,-0.0011269) >zz[131030] (0.999999,-0.00107626) >zz[131031] (0.999999,-0.00102678) >zz[131032] (1,-0.000978471) >zz[131033] (1,-0.000931322) >zz[131034] (1,-0.000885338) >zz[131035] (1,-0.000840519) >zz[131036] (1,-0.000796863) > >Identified that cout setprecision (15) shows me that the input values to the fft are correct zz[131064] (0.999999998888523,-4.71482053224406e-05) zz[131065] (0.999999999306111,-3.72529029760027e-05) zz[131066] (0.999999999593255,-2.8521753843732e-05) zz[131067] (0.999999999780449,-2.09547579273147e-05) zz[131068] (0.999999999894121,-1.45519152278533e-05) zz[131069] (0.999999999956632,-9.31322574602015e-06) zz[131070] (0.999999999986278,-5.2386894821881e-06) zz[131071] (0.99999999999729,-2.32830643653659e-06) zz[131072] (0.999999999999831,-5.82076609134641e-07) but when it comes out of the fft it strays, compared to matlab as said initially fz[131065] 1727.50831421296,-1515.52568733716 fz[131066] 1738.43586424091,-1580.81341761358 fz[131067] 1684.34619863036,-1563.24876927625 fz[131068] 1700.90551563079,-1621.25191184388 fz[131069] 1652.6897324293,-1596.62142662382 fz[131070] 1675.33028045191,-1647.73708634193 fz[131071] 1633.35368929937,-1616.37081795208 fz[131072] 1662.34262952489,-1660.87024640258 Column 131065 1726.12467509267 - 1515.86633403907i Column 131066 1737.27588448741 - 1581.31505759045i Column 131067 1683.02379469184 - 1563.44617159844i Column 131068 1699.81142853965 - 1621.62038834973i Column 131069 1651.44065609826 - 1596.69358159125i Column 131070 1674.30929408701 - 1647.99273838256i Column 131071 1632.17793684431 - 1616.34051114058i Column 131072 1661.40285267096 - 1661.02707675479i Found the problem is with the first value of the zz input zz[0] (0.999999999999831,-5.82076609134641e-07) << should be zero zz[1] (1,-0) zz[2] (0.999999999999831,-5.82076609134641e-07) zz[3] (0.99999999999729,-2.32830643653659e-06) with matlab 1 0.999999999999831 - 5.82076609134641e-007i 0.99999999999729 - 2.32830643653659e-006i this improves things fz[131070] 1674.33028045191,-1647.73708575986 fz[131071] 1632.35368929937,-1616.37081737 fz[131072] 1661.34262952489,-1660.8702458205 and suggests some other small error in the input data. Thanks
