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
On Sep 17, 3:10 pm, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote:> 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 > > IanWhy are you comparing the elements of a 131072 element array with the elements of a 131073 element array? Dale B. Dalrymple
Reply by ●September 18, 20092009-09-18
>On Sep 17, 3:10 pm, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote: >> 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 > >Why are you comparing the elements of a 131072 element array with the >elements of a 131073 element array? > >Dale B. Dalrymple >if you are refering to the fact that matlab starts with 1 and c++ starts with 0, I have already accounted for that so entries are in line with each other. See this followup thread http://www.dsprelated.com/showmessage/118116/1.php
Reply by ●September 18, 20092009-09-18
On Sep 18, 6:12 am, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote:> >On Sep 17, 3:10 pm, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote: > >> 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 > > >Why are you comparing the elements of a 131072 element array with the > >elements of a 131073 element array? > > >Dale B. Dalrymple > > if you are refering to the fact that matlab starts with 1 and c++ starts > with 0, I have already accounted for that so entries are in line with each > other. > See this followup thread > > http://www.dsprelated.com/showmessage/118116/1.phpA process that generates 131073 values isn't the same process as one that generates 131072 values. Setting the zero indexed value to zero and only looking at the last 5 values don't guarantee a fix. As Randy said: "Unless you show us these operations, I'm not sure how we can help you." Dale B. Dalrymple
Reply by ●September 18, 20092009-09-18
>On Sep 18, 6:12 am, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote: >> >On Sep 17, 3:10 pm, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote: >> >> Hi I'm using fftw3 on ubuntu linux in a convolution of1-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 curtailthis>> 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 >> >> >Why are you comparing the elements of a 131072 element array with the >> >elements of a 131073 element array? >> >> >Dale B. Dalrymple >> >> if you are refering to the fact that matlab starts with 1 and c++starts>> with 0, I have already accounted for that so entries are in line witheach>> other. >> See this followup thread >> >> http://www.dsprelated.com/showmessage/118116/1.php > >A process that generates 131073 values isn't the same process as one >that generates 131072 values. Setting the zero indexed value to zero >and only looking at the last 5 values don't guarantee a fix. > >As Randy said: > "Unless you show us these >operations, I'm not sure how we can help you." > >Dale B. Dalrymple >Hi Dale thanks , I'm checking the two arrays going into the fft. zz is ok now but yy needs to be shifted, and it has incorrect imaginary values. so I have more work to do.
Reply by ●September 21, 20092009-09-21
>>On Sep 18, 6:12 am, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote: >>> >On Sep 17, 3:10 pm, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote: >>> >> 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 >>> >>> >Why are you comparing the elements of a 131072 element array withthe>>> >elements of a 131073 element array? >>> >>> >Dale B. Dalrymple >>> >>> if you are refering to the fact that matlab starts with 1 and c++ >starts >>> with 0, I have already accounted for that so entries are in line with >each >>> other. >>> See this followup thread >>> >>> http://www.dsprelated.com/showmessage/118116/1.php >> >>A process that generates 131073 values isn't the same process as one >>that generates 131072 values. Setting the zero indexed value to zero >>and only looking at the last 5 values don't guarantee a fix. >> >>As Randy said: >> "Unless you show us these >>operations, I'm not sure how we can help you." >> >>Dale B. Dalrymple >> > >Hi Dale > >thanks , I'm checking the two arrays going into the fft. zz is ok nowbut>yy needs to be shifted, and it has incorrect imaginary values. so I have >more work to do. > > >Hi Dale you're so right about the sensitivity of the process values, and the differences in arrays. I've fixed this, so results come out OK, until the second time I run it through an inverse. I think its something to do with reuse of plans. Can you offer any advice on destroying plans or efficiency of plans. My code has a function which is reused in a loop and within that function creates and destroys 2 x 1d forward plans and 1 X 1d inverse plan. I just use fftw_destroy_plan(name) 3 times within the function. could anything be left behind distorting the second call to the inverse ? Thanks
Reply by ●September 22, 20092009-09-22
On Sep 21, 6:07�pm, "Kwisatch" <ibellinfan...@yahoo.co.uk> wrote:> you're so right about the sensitivity of the process values, and the > differences in arrays. I've fixed this, so results come out OK, until the > second time I run it through an inverse. I think its something to do with > reuse of plans. Can you offer any advice on destroying plans or efficiency > of plans. My code has a function which is reused in a loop and within that > function creates and destroys 2 x 1d forward plans and 1 X 1d inverse plan. > I just use fftw_destroy_plan(name) 3 times within the function. could > anything be left behind distorting the second call to the inverse ?No, you can create and destroy plans as many times as you want; you have some other bug in your code. (Note, however, that FFTW is optimized for the case where you create a plan once and re-use it for many transforms, so if you find yourself continually re-creating plans for FFTs of the same size then use are not using FFTW optimally. This is a question of efficiency, not correctness, however.) Regards, Steven G. Johnson
