Hi all, I'm running into some problems when using FFTs to interpolate data from coarse grids to finer grids. The operation that I'm doing to achieve this is: * Forward 3D FFT on coarse grid * Copy data to fine reciprocal space grid vector with zero padding * Backward 3D FFT The backward FFT is therefore of larger size than the forward FFT. Despite using double precision FFTs, I find that I lose some numerical precision when doing this, for example, starting from a normalized input vector, if I interpolate, then undo the interpolation, the resulting vector's norm is different from 1.0 by about 1E-10. This numerical precision is a lot worse than what results from a forward and backward FFT of the same size, which is more like 1E-15. For my particular application, I need high numerical precision. I have a suspicion that the issue is due to me using FFTW, which is probably selecting different FFT algorithms for the forward and backward FFT when interpolating, and perhaps using two different FFT algorithms sequentially increases numerical noise. Does anyone else interpolate in a similar manner to me, if so, do they see this kind of precision loss and is it avoidable? Many thanks, -- Will Glover UCLA Department of Chemistry & Biochemistry
numerical precision of FFTs when interpolating
Started by ●November 26, 2008
Reply by ●November 27, 20082008-11-27
On Nov 27, 6:51�am, "WillGlover" <will_glo...@yahoo.com> wrote:> Hi all, > I'm running into some problems when using FFTs to interpolate data from > coarse grids to finer grids. �The operation that I'm doing to achieve this > is: > * Forward 3D FFT on coarse grid > * Copy data to fine reciprocal space grid vector with zero padding > * Backward 3D FFT > > The backward FFT is therefore of larger size than the forward FFT. > > Despite using double precision FFTs, I find that I lose some numerical > precision when doing this, for example, starting from a normalized input > vector, if I interpolate, then undo the interpolation, the resulting > vector's norm is different from 1.0 by about 1E-10. �This numerical > precision is a lot worse than what results from a forward and backward FFT > of the same size, which is more like 1E-15. �For my particular application, > I need high numerical precision. > > I have a suspicion that the issue is due to me using FFTW, which is > probably selecting different FFT algorithms for the forward and backward > FFT when interpolating, and perhaps using two different FFT algorithms > sequentially increases numerical noise. > > Does anyone else interpolate in a similar manner to me, if so, do they see > this kind of precision loss and is it avoidable? > > Many thanks, > -- > Will Glover > UCLA Department of Chemistry & Biochemistryi dont think simply taking a FFT and zero padding and IFFT gives interpolation. Interpolation means zero padding in time domain and convolution with any kernel like 'sinc' for example.You can do the convolution process faster by using FFT/IFFT and overlap add method.
Reply by ●November 27, 20082008-11-27
On Nov 27, 9:16�am, rajesh <getrajes...@gmail.com> wrote:> On Nov 27, 6:51�am, "WillGlover" <will_glo...@yahoo.com> wrote: > > > > > > > Hi all, > > I'm running into some problems when using FFTs to interpolate data from > > coarse grids to finer grids. �The operation that I'm doing to achieve this > > is: > > * Forward 3D FFT on coarse grid > > * Copy data to fine reciprocal space grid vector with zero padding > > * Backward 3D FFT > > > The backward FFT is therefore of larger size than the forward FFT. > > > Despite using double precision FFTs, I find that I lose some numerical > > precision when doing this, for example, starting from a normalized input > > vector, if I interpolate, then undo the interpolation, the resulting > > vector's norm is different from 1.0 by about 1E-10. �This numerical > > precision is a lot worse than what results from a forward and backward FFT > > of the same size, which is more like 1E-15. �For my particular application, > > I need high numerical precision. > > > I have a suspicion that the issue is due to me using FFTW, which is > > probably selecting different FFT algorithms for the forward and backward > > FFT when interpolating, and perhaps using two different FFT algorithms > > sequentially increases numerical noise. > > > Does anyone else interpolate in a similar manner to me, if so, do they see > > this kind of precision loss and is it avoidable? > > > Many thanks, > > -- > > Will Glover > > UCLA Department of Chemistry & Biochemistry > > i dont think simply taking a FFT and zero padding and IFFT gives > interpolation. > > Interpolation means zero padding in time domain and convolution with > any > kernel like 'sinc' for example.You can do the convolution process > faster by > using FFT/IFFT and overlap add method.- Hide quoted text - > > - Show quoted text -by zero padding in time domain i mean inserting appropriate number of zeros in between samples.
Reply by ●November 27, 20082008-11-27
On Nov 26, 11:19�pm, rajesh <getrajes...@gmail.com> wrote:> On Nov 27, 9:16�am, rajesh <getrajes...@gmail.com> wrote: > > > > > > > On Nov 27, 6:51�am, "WillGlover" <will_glo...@yahoo.com> wrote: > > > > Hi all, > > > I'm running into some problems when using FFTs to interpolate data from > > > coarse grids to finer grids. �The operation that I'm doing to achieve this > > > is: > > > * Forward 3D FFT on coarse grid > > > * Copy data to fine reciprocal space grid vector with zero padding > > > * Backward 3D FFT > > > > The backward FFT is therefore of larger size than the forward FFT. > > > > Despite using double precision FFTs, I find that I lose some numerical > > > precision when doing this, for example, starting from a normalized input > > > vector, if I interpolate, then undo the interpolation, the resulting > > > vector's norm is different from 1.0 by about 1E-10. �This numerical > > > precision is a lot worse than what results from a forward and backward FFT > > > of the same size, which is more like 1E-15. �For my particular application, > > > I need high numerical precision. > > > > I have a suspicion that the issue is due to me using FFTW, which is > > > probably selecting different FFT algorithms for the forward and backward > > > FFT when interpolating, and perhaps using two different FFT algorithms > > > sequentially increases numerical noise. > > > > Does anyone else interpolate in a similar manner to me, if so, do they see > > > this kind of precision loss and is it avoidable? > > > > Many thanks, > > > -- > > > Will Glover > > > UCLA Department of Chemistry & Biochemistry > > > i dont think simply taking a FFT and zero padding and IFFT gives > > interpolation. > > > Interpolation means zero padding in time domain and convolution with > > any > > kernel like 'sinc' for example.You can do the convolution process > > faster by > > using FFT/IFFT and overlap add method.- Hide quoted text - > > > - Show quoted text - > > by zero padding in time domain i mean inserting appropriate number of > zeros in between samples.- Hide quoted text - > > - Show quoted text -(to rajesh): Well, ...., actually, you can zero pad in either domain (e.g.: for 1D forward/inverse, zero padding in frequency gives you interpolated time points). (to Will): This might be a Steven Johnson question (he's the FFTW expert), and he seems to find posts related to FFTW, so perhaps he'll find yours. If not, you might also try posting at the MATLAB group, since MATLAB also uses FFTW. As for the interpolation part, I do hope you realize that you don't really get more 'resolution' when you zero pad. After all, if you have 1000 data points, you could compute the DFT at a million different frequency points (e.g.: 'k' = .000001, .000002, etc.) to get a million points in the frequency domain. But your 'resolution' is still determined by the fact that you have a data record that spans only 1000 sample times (and, actually, resolution depends on the SNR, too). Sorry I can't be of more help with regard to the loss of precision, but I haven't seen that kind of deterioration in zero padding 2D FFTs. Kevin
Reply by ●November 27, 20082008-11-27
>On Nov 26, 11:19=A0pm, rajesh <getrajes...@gmail.com> wrote: >> On Nov 27, 9:16=A0am, rajesh <getrajes...@gmail.com> wrote: >> >> >> >> >> >> > On Nov 27, 6:51=A0am, "WillGlover" <will_glo...@yahoo.com> wrote: >> >> > > Hi all, >> > > I'm running into some problems when using FFTs to interpolate datafr=>om >> > > coarse grids to finer grids. =A0The operation that I'm doing toachie=>ve this >> > > is: >> > > * Forward 3D FFT on coarse grid >> > > * Copy data to fine reciprocal space grid vector with zero padding >> > > * Backward 3D FFT >> >> > > The backward FFT is therefore of larger size than the forward FFT. >> >> > > Despite using double precision FFTs, I find that I lose somenumerica=>l >> > > precision when doing this, for example, starting from a normalizedin=>put >> > > vector, if I interpolate, then undo the interpolation, theresulting>> > > vector's norm is different from 1.0 by about 1E-10. =A0Thisnumerical>> > > precision is a lot worse than what results from a forward andbackwar=>d FFT >> > > of the same size, which is more like 1E-15. =A0For my particularappl=>ication, >> > > I need high numerical precision. >> >> > > I have a suspicion that the issue is due to me using FFTW, whichis>> > > probably selecting different FFT algorithms for the forward andbackw=>ard >> > > FFT when interpolating, and perhaps using two different FFTalgorithm=>s >> > > sequentially increases numerical noise. >> >> > > Does anyone else interpolate in a similar manner to me, if so, dothe=>y see >> > > this kind of precision loss and is it avoidable? >> >> > > Many thanks, >> > > -- >> > > Will Glover >> > > UCLA Department of Chemistry & Biochemistry >> >> > i dont think simply taking a FFT and zero padding and IFFT gives >> > interpolation. >> >> > Interpolation means zero padding in time domain and convolution with >> > any >> > kernel like 'sinc' for example.You can do the convolution process >> > faster by >> > using FFT/IFFT and overlap add method.- Hide quoted text - >> >> > - Show quoted text - >> >> by zero padding in time domain i mean inserting appropriate number of >> zeros in between samples.- Hide quoted text - >> >> - Show quoted text - > >(to rajesh): Well, ...., actually, you can zero pad in either domain >(e.g.: for 1D forward/inverse, zero padding in frequency gives you >interpolated time points). > >(to Will): This might be a Steven Johnson question (he's the FFTW >expert), and he seems to find posts related to FFTW, so perhaps he'll >find yours. If not, you might also try posting at the MATLAB group, >since MATLAB also uses FFTW. > >As for the interpolation part, I do hope you realize that you don't >really get more 'resolution' when you zero pad. After all, if you >have 1000 data points, you could compute the DFT at a million >different frequency points (e.g.: 'k' =3D .000001, .000002, etc.) to get >a million points in the frequency domain. But your 'resolution' is >still determined by the fact that you have a data record that spans >only 1000 sample times (and, actually, resolution depends on the SNR, >too). > >Sorry I can't be of more help with regard to the loss of precision, >but I haven't seen that kind of deterioration in zero padding 2D FFTs. > >Kevin >Hi Kevin, Thanks for the reply. I know I won't get more resolution, I just want to interpolate the data in between the original grid points. Essentially, the data in frequency space is adequately sampled, but I need to linearly convolve it with a function that has higher frequency components. I've already emailed the fftw developers, but I thought I would also post here in case I don't get a response. -- Will
Reply by ●November 27, 20082008-11-27
On Nov 27, 2:13 am, "WillGlover" <will_glo...@yahoo.com> wrote:> >On Nov 26, 11:19=A0pm, rajesh <getrajes...@gmail.com> wrote: > >> On Nov 27, 9:16=A0am, rajesh <getrajes...@gmail.com> wrote: > > >> > On Nov 27, 6:51=A0am, "WillGlover" <will_glo...@yahoo.com> wrote: > > >> > > Hi all, > >> > > I'm running into some problems when using FFTs to interpolate data > fr= > >om > >> > > coarse grids to finer grids. =A0The operation that I'm doing to > achie= > >ve this > >> > > is: > >> > > * Forward 3D FFT on coarse grid > >> > > * Copy data to fine reciprocal space grid vector with zero padding > >> > > * Backward 3D FFT > > >> > > The backward FFT is therefore of larger size than the forward FFT. > > >> > > Despite using double precision FFTs, I find that I lose some > numerica= > >l > >> > > precision when doing this, for example, starting from a normalized > in= > >put > >> > > vector, if I interpolate, then undo the interpolation, the > resulting > >> > > vector's norm is different from 1.0 by about 1E-10. =A0This > numerical > >> > > precision is a lot worse than what results from a forward and > backwar= > >d FFT > >> > > of the same size, which is more like 1E-15. =A0For my particular > appl= > >ication, > >> > > I need high numerical precision. > > >> > > I have a suspicion that the issue is due to me using FFTW, which > is > >> > > probably selecting different FFT algorithms for the forward and > backw= > >ard > >> > > FFT when interpolating, and perhaps using two different FFT > algorithm= > >s > >> > > sequentially increases numerical noise. > > >> > > Does anyone else interpolate in a similar manner to me, if so, do > the= > >y see > >> > > this kind of precision loss and is it avoidable? > > >> > > Many thanks, > >> > > -- > >> > > Will Glover > >> > > UCLA Department of Chemistry & Biochemistry > > >> > i dont think simply taking a FFT and zero padding and IFFT gives > >> > interpolation. > > >> > Interpolation means zero padding in time domain and convolution with > >> > any > >> > kernel like 'sinc' for example.You can do the convolution process > >> > faster by > >> > using FFT/IFFT and overlap add method.- Hide quoted text - > > >> > - Show quoted text - > > >> by zero padding in time domain i mean inserting appropriate number of > >> zeros in between samples.- Hide quoted text - > > >> - Show quoted text - > > >(to rajesh): Well, ...., actually, you can zero pad in either domain > >(e.g.: for 1D forward/inverse, zero padding in frequency gives you > >interpolated time points). > > >(to Will): This might be a Steven Johnson question (he's the FFTW > >expert), and he seems to find posts related to FFTW, so perhaps he'll > >find yours. If not, you might also try posting at the MATLAB group, > >since MATLAB also uses FFTW. > > >As for the interpolation part, I do hope you realize that you don't > >really get more 'resolution' when you zero pad. After all, if you > >have 1000 data points, you could compute the DFT at a million > >different frequency points (e.g.: 'k' =3D .000001, .000002, etc.) to get > >a million points in the frequency domain. But your 'resolution' is > >still determined by the fact that you have a data record that spans > >only 1000 sample times (and, actually, resolution depends on the SNR, > >too). > > >Sorry I can't be of more help with regard to the loss of precision, > >but I haven't seen that kind of deterioration in zero padding 2D FFTs. > > >Kevin > > Hi Kevin, > Thanks for the reply. I know I won't get more resolution, I just want to > interpolate the data in between the original grid points. Essentially, the > data in frequency space is adequately sampled, but I need to linearly > convolve it with a function that has higher frequency components. > I've already emailed the fftw developers, but I thought I would also post > here in case I don't get a response. > > -- > WillWill, Zero padding in the frequency domain without applying any type of filter results in ripples in your time/space domain due to wrap around of sinc() function. This is more likely the cause of the artifacts than the FFT algorithm. Cheers, David
Reply by ●November 27, 20082008-11-27
>Will, Zero padding in the frequency domain without applying any type >of filter results in ripples in your time/space domain due to wrap >around of sinc() function. This is more likely the cause of the >artifacts than the FFT algorithm. > >Cheers, >David >Yep I know that, but those artifacts should not be causing what I see (at least in principle). Let me explain again and I'll try to modify my language to be more DSP-like (I'm not in the field of DSP as you can probably tell ;) ), even though I'm not dealing with signals... Let A be my initial signal represented on N points. Let it be normalized, so that if treated as a vector, |A|=1.0 I apply an FFT to get the spectrum of A: B=[FFT]A, represented by N points in the frequency domain (I know, FFTW stores the output vector on just N/2+1 frequencies in real-to-complex transforms). I zero-pad B out to 2N frequencies, I'll call this new spectrum C. Then I IFFT C to get D on 2N points and finally scale D by 1/(sqrt(N)*sqrt(2N)) since my FFTs were not normalized. Now I'm seeing that abs( |D|-1.0 ) = 1.0E-10 Compared to if I did G=IFFT FFT A, abs( |G|-1.0 ) = 1.0E-15, without padding and with FFTs of the same size. This is what I mean by loss of precision. I am reluctant to say it's round-off error, because I'm not sure it's that, but it could well be, especially as the Norm changes with different input vectors. To remind, I'm using double precision throughout. I've also come across a DSP buzzword, "quantization", could this be related?
Reply by ●November 27, 20082008-11-27
On 2008-11-27 13:17:38 -0400, "WillGlover" <will_glover@yahoo.com> said:> >> Will, Zero padding in the frequency domain without applying any type >> of filter results in ripples in your time/space domain due to wrap >> around of sinc() function. This is more likely the cause of the >> artifacts than the FFT algorithm. >> >> Cheers, >> David >> > > Yep I know that, but those artifacts should not be causing what I see (at > least in principle). > Let me explain again and I'll try to modify my language to be more > DSP-like (I'm not in the field of DSP as you can probably tell ;) ), even > though I'm not dealing with signals... > > Let A be my initial signal represented on N points. Let it be normalized, > so that if treated as a vector, |A|=1.0Did you actually check that |A| - 1.0 is small? (Not just printing |A| to 5 decimal places.)> I apply an FFT to get the spectrum of A: B=[FFT]A, represented by N points > in the frequency domain (I know, FFTW stores the output vector on just > N/2+1 frequencies in real-to-complex transforms).How close to 1.0 is |B|? Given that the FFT is unitary it should be 1.0 also.> I zero-pad B out to 2N frequencies, I'll call this new spectrum C. > > Then I IFFT C to get D on 2N points and finally scale D by > 1/(sqrt(N)*sqrt(2N)) since my FFTs were not normalized. > > Now I'm seeing that abs( |D|-1.0 ) = 1.0E-10How close is the interpolation on the original points (after correcting the normailization)? What is the norm of this decimate? What is N? It may be that some roundoff should be expected but 1.0e4 certainly seems excessive.> Compared to if I did G=IFFT FFT A, abs( |G|-1.0 ) = 1.0E-15, without > padding and with FFTs of the same size. > > This is what I mean by loss of precision. I am reluctant to say it's > round-off error, because I'm not sure it's that, but it could well be, > especially as the Norm changes with different input vectors. To remind, > I'm using double precision throughout. > I've also come across a DSP buzzword, "quantization", could this be > related?
Reply by ●November 27, 20082008-11-27
On 2008-11-27 13:54:49 -0400, Gordon Sande <g.sande@worldnet.att.net> said:> On 2008-11-27 13:17:38 -0400, "WillGlover" <will_glover@yahoo.com> said: > >> >>> Will, Zero padding in the frequency domain without applying any type >>> of filter results in ripples in your time/space domain due to wrap >>> around of sinc() function. This is more likely the cause of the >>> artifacts than the FFT algorithm. >>> >>> Cheers, >>> David >>> >> >> Yep I know that, but those artifacts should not be causing what I see (at >> least in principle). >> Let me explain again and I'll try to modify my language to be more >> DSP-like (I'm not in the field of DSP as you can probably tell ;) ), even >> though I'm not dealing with signals... >> >> Let A be my initial signal represented on N points. Let it be normalized, >> so that if treated as a vector, |A|=1.0 > > Did you actually check that |A| - 1.0 is small? > (Not just printing |A| to 5 decimal places.) > >> I apply an FFT to get the spectrum of A: B=[FFT]A, represented by N points >> in the frequency domain (I know, FFTW stores the output vector on just >> N/2+1 frequencies in real-to-complex transforms). > > How close to 1.0 is |B|? Given that the FFT is unitary it should be > 1.0 also. > >> I zero-pad B out to 2N frequencies, I'll call this new spectrum C. >> >> Then I IFFT C to get D on 2N points and finally scale D by >> 1/(sqrt(N)*sqrt(2N)) since my FFTs were not normalized. >> >> Now I'm seeing that abs( |D|-1.0 ) = 1.0E-10 > > How close is the interpolation on the original points (after > correcting the normailization)? What is the norm of this decimate? > > What is N? It may be that some roundoff should be expected but 1.0e4 > certainly seems excessive.Forgot to ask.. How big is the N/2 coefficient and did you remember to spilt it into two pieces when zero padding in reciprocal space?>> Compared to if I did G=IFFT FFT A, abs( |G|-1.0 ) = 1.0E-15, without >> padding and with FFTs of the same size. >> >> This is what I mean by loss of precision. I am reluctant to say it's >> round-off error, because I'm not sure it's that, but it could well be, >> especially as the Norm changes with different input vectors. To remind, >> I'm using double precision throughout. >> I've also come across a DSP buzzword, "quantization", could this be >> related?
Reply by ●November 28, 20082008-11-28
Ok, I've found the problem, thanks to the useful suggestions. It's not the FFT, it's how I am zero padding, or rather stuffing. I've found the following guide on zero-stuffing: http://www.dspguru.org/howto/tech/zeropad.htm It tells me to stuff the zeros after N/2 frequencies, i.e. just before the nyquist frequency of the original sample. Now if I do just this, the resulting stuffed spectrum will no longer have Hermitian symmetry because the value at the original nyquist frequency (at index N/2 +1) has been replaced by zero and the value at the negative of the nyquist frequency (at index N/2+1 + (M-N), where M is the size of the larger FFT)). Now when I IFFT and then just take the real component, I'll lose some information. In the example on that page, the amplitude at the nyquist frequency was zero, but that's not always going to be the case... It "feels" like the correct thing to do when stuffing is to leave the value at the original Nyquist frequency intact and copy its conjugate into the negative Nyquist frequency value. This retains Hermitian symmetry, but will change the norm of the resulting vector, so do I divide these Nyquist frequencies by sqrt(2)? Hmmm, starting to get arbitrary... Most of the information I can find online is for zero-padding time signals (which is a lot easier, since you just add zeros to the end). Any further suggestions? Much appreciated, -- Will Glover
