In the standard definition of Discrete Fourier transform, as in Wikipedia, only positive frequencies exist. I want zero frequency to be at the center of the spectrum. In such a case, how should I change formula for Discrete Fourier transform?
2 Answers
Please don't change the formula, people have given it some thought. If you look at it
$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$
you can see that
$$X[k]=X[k+N]\tag{2}$$
holds (because $e^{-j2\pi nk/N}$ is $N$-periodic in $k$ and $n$). Consequently, values corresponding to negative indices $k$ can simply be derived from positive values by increasing the negative index by the DFT length $N$. E.g.,
$$X[-1]=X[N-1]\tag{3}$$
et cetera.
Negatives exist also; $k=0$ to $N - 1$ spans positive frequencies up to $N/2$, then negatives.
To center the zero bin, shift indexing into the complex sinusoid:
$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j2\pi n(k + N/2)/N} $$
equivalent to fftshift
. Python code below. (For "positives only", there's real DFT.)
import numpy as np
from numpy.fft import fft, fftshift
def dft(x, center_zero=False):
N = len(x)
offset = np.ceil(N / 2) if center_zero else 0
out = np.zeros(N, dtype='complex128')
for k in range(N):
for n in range(N):
out[k] += x[n] * np.exp(-2j*np.pi * n * (k + offset) / N)
return out
for N in (128, 129): # check even & odd case
x = np.random.randn(N)
xf0 = fftshift(fft(x))
xf1 = dft(x, center_zero=True)
assert np.allclose(xf0, xf1)