How to do simple extrapolation with Fourier transformation?

Question

I have 1024 sample points, and I would like to do really simple extrapolation using Fourier transformation. First I apply Fast fourier transformation on the data. My first intuition was that I just calculate the inverse fourier transformation on a larger interval. The result of the FFT is a collection of sinusoidal waves, so I expected that it's result is continuous at the end of the original interval.

I calculate the result with this formula:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{\frac{i 2 \pi n k}{N}} $$

The extrapolation at the end of the interval looks like this:

I would like to avoid this "break" in the result.

The python source code is:

import numpy as np
import matplotlib.pyplot as plt

from scipy.fftpack import fft

# generate samples
data = np.zeros(1024)
for k in range(0, 1024):
    data[k] = np.sin(k/20 * np.pi*2)


# apply FFT on samples
waves = fft(data)

# calculate the new data points from the FFT result
res = np.zeros(2048)
for k in range(0, 2048):
    val = 0.0
    for n in range(0, len(waves)):
        # https://dsp.stackexchange.com/a/510/25943
        val += waves[n]*np.exp(1.j * 2*np.pi * n * k / len(waves)) / len(waves)
    res[k] = val.real

#plot the result
plt.plot(res)
plt.show()

that "break" must be due to circular wrap-around or something, because there an extra spacing of the peaks and valleys of the sinuosoid. it's not just a simple noisey click to fix. in the pitch-shifting biz, we sometimes call a delay mismatch like that a "glitch". to mend that glitch, you need to slide the whole signal on the right over a bit. dunno exactly how Fourier helps unless it's Short-Time Fourier (the STFT). — robert bristow-johnson, Oct 29 '17 at 18:20

Fat32 · Answer 1 · 2021-10-28T17:57:37.220

Your computation is right, but your expectation is wrong.

The method that you use to extend the periodic sequence $\tilde{x}[n]$, of period $N$, outside of its base period, $0 \leq n <N$, into a range, say, $0 \leq n <2N$ will simply copy the base period $x[n]$ next to itself. Unless the last and the first elements, $x[N-1]$ and $x[0]$, are the same there will be a mismatch (discontinuity), that you interpret as a glitch, at the period boundaries.

The inverse DFT that you implement is:

$$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} k n}$$

And when evaluated for computing the extrapolated samples, $x[n_0 + N]$, with $n_0$ in $0 \leq n < N$, you will have:

$$x[n_0 + N] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} k (n_0 + N)} = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} kn_0} e^{j \frac{2 \pi}{N} N} = x[n_0] $$

as expected, due to the periodicity of the twiddle (complex exponential) in this case.

How to do simple extrapolation with Fourier transformation?

1 Answers1