numerical precision of FFTs when interpolating

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.

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

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