"extended" IFFT

Question

If I have a waveform x such as

x = [math.sin(W*t + Ph) for t in range(16)]

with arbitrary W and Ph, and I calculate its (Real) FFT f with

f = numpy.fft.rfft(x)

I can get the original x with

numpy.fft.irfft(f)

Now, what if I need to extend the range of the recovered waveform a number of samples to the left and to the right? I.e. a waveform y such that len(y) == 48, y[16:32] == x and y[0:16], y[32:48] are the periodic extensions of the original waveform.

In other words, if the FFT assumes its input is an infinite function f(t) sampled over t = 0, 1, ... N-1, how can I recover the values of f(t) for t<0 and t>=N?

Note: I used a perfect sine wave as an example, but in practice x could be anything: arbitrary signals such as x = range(16) or x = np.random.rand(16), or a segment of any length taken from a random .wav file.

Carree

You might want to ask this question in dsp.stackexchange.com too, just in case. – Licit

Don't ask a question on two forums. If you realized you asked in the wrong place, delete and ask somewhere else, or ask a moderator to move it. meta.stackexchange.com/questions/64068/… – Stoichiometry

Answer 1

Now, what if I need to extend the range of the recovered waveform a number of samples to the left and to the right? I.e. a waveform y such that len(y) == 48, y[16:32] == x and y[0:16], y[32:48] are the periodic extensions of the original waveform.

The periodic extension are also just x because it's the periodic extension.

In other words, if the FFT assumes its input is an infinite function f(t) sampled over t = 0, 1, ... N-1, how can I recover the values of f(t) for t<0 and t>=N?

The "N-point FFT assumes" that your signal is periodic with a periodicity of N. That's because all the harmonic base functions your block is decomposed into are periodic in the way that the previous N and succeding N samples are just a copy of the main N samples.

If you allow any value for W your input sinusoid won't be periodic with periodicity of N. But that does not stop the FFT function from decomposing it into a sum of many periodic sinusiods. And the sum of periodic sinusoids with periodicity of N will also have a periodicity of N.

Clearly, you have to rethink the problem.

Maybe you could make use of linear prediction. Compute a couple of linear prediction coefficients based on your fragment's windowed auto-correlation and the Levinson-Durbin recursion and extrapolate using those prediction coefficients. However, for a stable prediction filter, the prediction will converge to zero and the speed of convergence depends on what kind of signal you have. The perfect linear prediction coefficients for white noise, for example, are all zero. In that case you would "extrapolate" zeros to the left and the right. But there's not much you can do about it. If you have white noise, there is just no information in your fragment about surrounding samples because all the samples are independent (that's what white noise is about).

This kind of linear prediction is actually able to predict sinusoid samples perfectly. So, if your input is sin(W*t+p) for arbitrary W and p you will only need linear prediction with order two. For more complex signals I suggest an order of 10 or 16.

Wrangle

It looks like linear predictors are the kind of thing I was looking for. My question was not clearly formulated and my ignorance misled me to thing FFT was the tool for the job, but you somehow read through all that noise and found out what my problem was. Thank you! – Carree

Answer 2

The following examples should give you a good idea of how to go about it:

>>> x1 = np.random.rand(4)
>>> x2 = np.concatenate((x1, x1))
>>> x3 = np.concatenate((x1, x1, x1))
>>> np.fft.rfft(x1)
array([ 2.30410617+0.j        , -0.89574460-0.26838271j, -0.26468792+0.j        ])
>>> np.fft.rfft(x2)
array([ 4.60821233+0.j        ,  0.00000000+0.j        ,
       -1.79148921-0.53676542j,  0.00000000+0.j        , -0.52937585+0.j        ])
>>> np.fft.rfft(x3)
array([ 6.91231850+0.j        ,  0.00000000+0.j        ,
        0.00000000+0.j        , -2.68723381-0.80514813j,
        0.00000000+0.j        ,  0.00000000+0.j        , -0.79406377+0.j        ])

Of course the easiest way to get three periods is to concatenate 3 copies of the inverse FFT in the time domain:

np.concatenate((np.fft.irfft(f),) * 3)

But if you want or have to do this in the frequency domain, you can do the following:

>>> a = np.arange(4)
>>> f = np.fft.rfft(a)
>>> n = 3
>>> ext_f = np.zeros(((len(f) - 1) * n + 1,), dtype=f.dtype)
>>> ext_f[::n] = f * n
>>> np.fft.irfft(ext_f)
array([ 0.,  1.,  2.,  3.,  0.,  1.,  2.,  3.,  0.,  1.,  2.,  3.])

Lance

Thank you. Will this solution still work if the original waveform is not phase-aligned? Say, x = [math.sin(W*i + Ph) for i in range(16)] with arbitrary W and Ph. Time-domain concatenation certainly won't. – Carree

What about arbitrary waveforms such as x = range(16) or x = np.random.rand(16)? – Carree

@MartinBlech I think sellibitze's answer covers your new question in this comment. FFTs are great for periodic functions, but if you want to feed it [0, 1, 2, 3] and get back [-4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7] then you are going to have to look somewhere else. And the first thing you need to do is to define your question in an unambiguous way. FFT does the periodic extension by concatenating the original data, which is a very reasonable way of going about it. If that is not what you want, how exactly do you want it done? – Lance

I was unable to formulate my question correctly due to my own ignorance. From sellibitze's answer, it looks like linear predictors are the thing I was looking for but didn't know the name of. – Carree

Answer 3

For stationary waveforms that are periodic in the FFT aperture or length, you can just cyclicly repeat the waveform, or the IFFT(FFT()) re-synthesized equivalent waveform, to extend them in the time domain. For waveforms which are widowed in time from sources that are not periodic in the FFT aperture or length, the FFT result will be the spectrum convolved with a Sinc function. So some sort of equivalent to a de-convolution will be required to recover the original un-windowed spectral content. Since this deconvolution is difficult or impossible, most commonly an analysis/re-synthesis method is used instead, such as a phase-vocoder process or other frequency estimators. Then those estimated frequencies, which may be different from those in the bins of a single raw FFT result, can be fed to a bank of sinusoidal synthesizers, a mix of phase-modified IFFTs, or other re-synthesis methods, to create a longer waveform with approximately the same spectral content.

Wivestad