3
$\begingroup$

I've got a situation where I'd like to use an FFT to do interpolation in time on some complex data (I need to go to the frequency domain anyways to window my data).

The notional way of doing this is to take an FFT, which, in the case of something like FFTW will put the frequency components "in-order" which can be viewed as having positive frequencies first, followed by negative.

$\mathcal{F}([\begin{array}{ccccccc} a & b & c & d & e & f & g & h \end{array}]) = \begin{array}{cccccccc} [A & B & C & D & E & F & G & H] \end{array}$

Where ABCD are positive frequency components, and EFGH are negative frequency components.

You can split the frequency response in two and zero-pad between D and E, but you have to be careful if your transform is an odd length, as well as possibly splitting the nyquist bin, and it's more cumbersome from a buffer-juggling perspective too. It'd be better if I could pre-multiply by Fs/2, which will rotate the spectrum:

$\mathcal{F}(exp(j2\pi fs/2t) * [\begin{array}{ccccccc} a & b & c & d & e & f & g & h \end{array}]) = \begin{array}{cccccccc} [E & F & G & H & A & B & C & D] \end{array}$

Then it's just a matter of zero padding on the end, and getting back to the time domain. This might be as simple as just transforming into a sufficiently large zeroed buffer.

But, it's the getting back to the time domain that's got me wrapped around the axle. You can certainly inverse-FFT, and then multiply by -Fs/4 to shift things back to the way they should be, but I'd prefer to avoid the post-multiplication if I can (performance reasons).

My gut says that the shifted spectrum is now in the proper order (again in FFTW land) to just do a forward-FFT on it, but that might time-reverse the data, so perhaps a conjugate before transforming would be in order, but I'm not sure. Can anyone help me puzzle this out?

Improve this question
$\endgroup$
5
  • $\begingroup$ Where you suggest "zero padding on the end", do you mean to zero pad on both ends? Seems like if you had, for example, [E F G H A B C D 0 0 0 0] and then inverse transformed, you'd get a mess. $\endgroup$
    – The Photon
    Dec 1, 2012 at 6:35
  • $\begingroup$ A good point, I've had it in my head that I could zero pad on the end and then forward FFT, but I'll have to think on that. $\endgroup$
    – gct
    Dec 1, 2012 at 16:18
  • $\begingroup$ I think that this thread will be helpful to you. dsp.stackexchange.com/questions/2962/… $\endgroup$
    – Jim Clay
    Dec 2, 2012 at 2:28
  • $\begingroup$ @The Photon: After thinking on it more, I think it's OK, tuning by Fs/2 centers the spectrum, then you can zero pad on the end, and you're in-order to use an FFT for the inverse transform (conj(fft(conj(x)))). Then you have to tune by Fs/2/interp to re-center the zero padded spectrum. Since there's two tunes in there, maybe not a win vs buffer twiddling to zero pad. $\endgroup$
    – gct
    Dec 2, 2012 at 15:52
  • $\begingroup$ You may have a look at dsp.stackexchange.com/questions/79356 for analytic derivation. $\endgroup$
    – Eric Johnson
    Jan 22, 2022 at 13:33

2 Answers 2

Reset to default
5
$\begingroup$

It's important to remember that, for an even N, you have both F(0) and F(N/2) points, neither of which has a counterpart. So, for instance, with N = 8, you have:

F(0), F(1), F(2), F(3), F(4), F(-3), F(-2), F(-1)

(some people consider the N/2 point as a negative frequency ( ie: F(-4) ), but we'll consider it to be positive)

So to interpolate an 8 point FFT by a factor of 2, you'd do a 16 point inverse FFT using:

F(0), F(1), F(2), F(3), F(4)/2, 0, 0, 0, 0, 0, 0, 0, F(4)/2, F(-3), F(-2), F(-1)

So you split the N/2 point of your forward FFT, and add the appropriate number of zeroes in the center, and do the inverse.

Now, for an odd number N, you DON'T have an N/2 point (but you still have a unique 0 point), and the splitting business doesn't happen (but you still have to add the appropriate zeroes in the center).

I suggest you experiment with the FFT/IFFT method by using: a real input time domain waveform; a complex input time domain waveform; even and odd N; an integer number of cycles in the input; a non-integer number of cycles in the input. Then compare/contrast your results from the FFT/IFFT method with to those generated by a strictly time domain based interpolation (eg: use A*cos(omega*t) for the FFT/IFFT method and compare to A*cos(omega*t/2) generated in the time domain). Experience can be enlightening.

As for your multiplying method – I don't think it will work for an even N. And the FFT/IFFT timing, unless you're using very small FFT's, will probably dominate the total run time (of, course, you should always profile your programs to make sure).

Improve this answer
$\endgroup$
10
  • $\begingroup$ I've hacked together a small program using a complex chirp to play around with. The best I've come up with is: multiply by Fs/2, zero-pad at the end, use conj(FFT(conj(x))) to inverse-transform (because data is in-order for FFT and shifted for IFFT in fftw), then tune by Fs/2/interp to re-center the spectrum. All-told probably not a win over buffer twiddling. $\endgroup$
    – gct
    Dec 2, 2012 at 15:54
  • $\begingroup$ Two quick things: When you split the Fs/2 bin, should you conjugate it as well as dividing by two? And using this: pastebin.com/aLw6Jmtn code for splitting the spectrum, it seems to work well, except in the odd N case, I get pretty strong ringing at the beginning/end of my signal on the real channel (but not imaginary). Thoughts? $\endgroup$
    – gct
    Dec 2, 2012 at 18:40
  • $\begingroup$ If your input was real only, then you want to make sure that your FFT displays conjugate symmetry. As for your code - consider comparing your results with interpolated numbers generated strictly in the time domain. You will find that, when doing FFT/IFFT interpolation, half your output will be the original signal, and the in-between points will be the interpolated ones. $\endgroup$
    – Kevin McGee
    Dec 2, 2012 at 20:40
  • $\begingroup$ Boy I'm puzzled, as far as I can tell I'm doing the right thing. I've explicitly split my code for even/odd now to help me get a better grip on it: pastebin.com/nngQaSTB You can see some images of what I'm getting here: imgur.com/zjn1Q,hq7as,dJ1Fr,smyGU (pink is results from even-sized transform, green is odd sized). $\endgroup$
    – gct
    Dec 2, 2012 at 23:19
  • $\begingroup$ That graph of the chirp doesn't look like a chirp. Perhaps try a simpler, known waveform to start out with - say a sinusoid. You know what the waveform should be in the time domain, and what the interpolated waveform should be - it'll be easier to plot and test. That way you can be more sure of your code. $\endgroup$
    – Kevin McGee
    Dec 3, 2012 at 4:50
1
$\begingroup$

This can be a tricky problem for complex-valued data samples. If you have MATLAB, the interpft() function, whose code is viewable with edit interpft, shows how to interpolate. Be careful though, it assumes your data are periodic, just as the DFT/FFT do.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.