Thanks for the helpful comments everyone.  As correctly predicted by
Gordon, my issue was dealing with the Nyquist frequency incorrectly.  I was
neglecting to copy the nyquist component to both negative and positive
frequencies (and dividing by sqrt 2) in the padded spectrum, since I had
even FFTs.  This resulted in my spectrum being non-Hermitian, and then when
I took a complex-to-real IFFT I lost some of the data, hence my norm wasn't
quite 1.0.  Dealing with the Nyquist frequencies gets a little tricky with
3D data, but I managed to get it coded and it's all working beautifully
now.  I can interpolate and undo the interpolation and the vectors match
within 1.0E-15.  

As for the Gibb's phenomenon stuff, that's expected, and I don't want to
filter.  Essentially what I'm doing this for is to do a linear convolution
in reciprocal space via the usual Pad, FFT, diagonal matrix multiply,
IFFT.

Thanks for all the help folks!
-- 
Will

I looked at that post I wrote at 3:51 AM, and I have no idea why I
added that comment about adding zeros between the frequency data.  I
think it had to do with something I once did when I was trying to
accomplish some entirely different task.

So today I ran a simple C++ test program to see what I could find.
This is how it works:  I generate 8 points of a sawtooth wave (sum of
4 sine waves of different amplitudes and frequencies 0f to 3f).  Then
I generate 16 points over the same time span to see what I should get
when I use the FFT to interpolate the 8 point data.  Then I did an 8
point FFT on the 8 point data.  I get the proper amplitudes and phases
at points 0, 1, 2, 3 and 5, 6, 7.   The point at 4 (fs/2) is zero
because my data is exactly periodic within the 8 points.  I zero pad
by adding 8 zeros to the center of the data, so points 0, 1, 2, and 3
remain the same, and points 5, 6, and 7 become points 13, 14, and 15.
Then I inverse FFT with 16 points.

Everything works out perfectly - the 16 points obtained via FFT
interpolation are exactly the same as the points I expected to get
(the first 15 decimal digits of both data sets are identical, and my
imaginary parts via the FFT are all zero up to 15 digits as well). So
far, so good (now the fun begins).  I know that things won't work out
as well when I have input data that isn't periodic in the sampling
interval, so I changed the sawtooth input to get a non-periodic data
set.

As expected, I got a non-zero FFT output at point 4.  So I have 8
outputs and they all have non-zero real and imaginary parts (except
for point 4, which has a zero imaginary part, as it should, since my
input data is real). To zero pad this, I add 8 zeros after the 4
point.  I also split the real part of point 4 between point 4 and
point 12 (this is what Gordon's post was referring to).  So now I have
my zero padded data in points 0, 1, 2, 3, and 4(half of what it was),
and points 12 (half of the former 4 point), 13, 14, and 15. All the
other points are zero.  After inverse transforming, my imaginary
results are all zero up to 15 decimal digits, which leads me to
believe that the algorithm is correct.  The even numbered outputs are
excellent (they differ by, at most, the number in decimal digit 16).
But the odd numbered points are an entirely different story. One of
them starts to differ from the expected result in the first point to
the right of the decimal point.  More than likely, this is a Gibbs
type 'rippling' effect (which has already been pointed out by other
posters).

I've seen it before when using certain types of image processing
filters, but I never expected to see so much of it when doing
interpolation.

So it would seem that if your data is not periodic within the sampling
interval, then the interpolation process is subject to the Gibbs
effect, which will cause output errors.  I suppose you could window
your 3D data before the FFTs - that would help to minimize the errors.

I�ll check my code again just to make sure (and I might try some
windowing).

I was just about to post when I saw yours.

This was what I was going to say:  Gordon's post brings up an
excellent point that may be the source of your problem.  I presumed
that when you zero pad your frequency data that you're adding zeros
between your output points.  Perhaps you're doing it differently by
appending zeros to your frequency data and then inverse transforming?
That could very well be the problem.  How  exactly are you zero
padding your frequency data?

OK.  I think your post answered it - you're appending zeros.  That's
fine if you're zero padding a time domain sequence, but the frequency
domain is different.  By adding zeros between your points, you're
keeping all the frequency domain points in the same relative positions
when you use a doubled N.

So try adding zeros between your frequency domain points, and see if
that solves your problem.  I'll be back later tonight when I rest up
from the (too much) food I had earlier.

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

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?

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.

> 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?

>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?

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.
>
> --
> Will

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

>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.

-- 
Will

On Nov 26, 11:19&#4294967295;pm, rajesh <getrajes...@gmail.com> wrote:
> On Nov 27, 9:16&#4294967295;am, rajesh <getrajes...@gmail.com> wrote:
>
>
>
>
>
> > On Nov 27, 6:51&#4294967295;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. &#4294967295;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. &#4294967295;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. &#4294967295;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