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:

extrapolation error

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
  • 2
    $\begingroup$ 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). $\endgroup$ Oct 29 '17 at 18:20

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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