Shifted discrete Fourier transform

Forward transform

We define the shifted discrete Fourier transform as

\[Z_k \equiv \sum_{n = 0}^{N - 1} z_n \twiddle{- 2 \pi}{\left( n + \frac{1}{2} \right) k}{N} \equiv \dft{k}{N}{z_0}{z_1}{z_{N - 1}},\]

where \(\seq{k}{0}{1}{N - 1}\).

Note

Since this relation leads to

\[Z_k = \twiddle{- \pi}{k}{N} \sum_{n = 0}^{N - 1} z_n \twiddle{- 2 \pi}{n k}{N},\]

the shifted transform can be computed using the normal discrete Fourier transform followed by multiplication by a pre-factor.

By utilizing the decimation in time, the forward transform leads to

\[\begin{split}\begin{pmatrix} Z_{k } \\ Z_{k + N / 2} \\ \end{pmatrix} = \begin{pmatrix} \twiddle{\pi}{k}{N} & \twiddle{- \pi}{k}{N} \\ - I \twiddle{\pi}{k}{N} & I \twiddle{- \pi}{k}{N} \\ \end{pmatrix} \begin{pmatrix} \dft{k}{N / 2}{z_0}{z_2}{z_{N - 2}} \\ \dft{k}{N / 2}{z_1}{z_3}{z_{N - 1}} \\ \end{pmatrix}\end{split}\]

for \(\seq{k}{0}{1}{N / 2 - 1}\).

src/sdft.c

static int sdft (
    const size_t nitems,
    const size_t stride,
    const trig_t * const trigs,
    const double complex * const xs,
    double complex * const ys
) {
  if (1 == nitems) {
    ys[0] = xs[0];
  } else if (0 == nitems % 2) {
    const size_t nh = nitems / 2;
    // divide-and-conquer
    sdft(nh, stride * 2, trigs, xs         , ys     );
    sdft(nh, stride * 2, trigs, xs + stride, ys + nh);
    // unify two sub problem results
    for (size_t i = 0; i < nh; i++) {
      const size_t j = i + nh;
      const trig_t * const trig = trigs + stride * i;
      const double complex twiddle[2] = {
        trig->cos + I * trig->sin,
        trig->cos - I * trig->sin,
      };
      double complex * const yi = ys + i;
      double complex * const yj = ys + j;
      const double complex ye = *yi;
      const double complex yo = *yj;
      *yi =       twiddle[0] * ye +     twiddle[1] * yo;
      *yj = - I * twiddle[0] * ye + I * twiddle[1] * yo;
    }
  } else {
    // naive O(N^2) sdft
    for (size_t k = 0; k < nitems; k++) {
      double complex * const y = ys + k;
      *y = 0. + I * 0.;
      for (size_t n = 0; n < nitems; n++) {
        *y += xs[stride * n] * cexp(- 2. * pi * (n + 0.5) * k * I / nitems);
      }
    }
  }
  return 0;
}

Backward transform

The inverse transform is defined as

\[z_n = \frac{1}{N} \sum_{k = 0}^{N - 1} Z_k \twiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{N} \equiv \idft{n}{N}{Z_0}{Z_1}{Z_{N - 1}}.\]

Note

Due to

\[z_n = \frac{1}{N} \sum_{k = 0}^{N - 1} \left( Z_k \exp \left( \frac{\pi k}{N} I \right) \right) \twiddle{2 \pi}{n k}{N},\]

the inverse transform of a shifted signal can be computed using the normal inverse discrete Fourier transform where the input signal is multiplied by a pre-factor beforehand.

By utilizing the decimation in time, the inverse transform leads to

\[\begin{split}\begin{pmatrix} z_{n } \\ z_{n + N / 2} \\ \end{pmatrix} = \begin{pmatrix} \frac{1}{2} & \frac{1}{2} \twiddle{\pi}{2 n + 1}{N} \\ \frac{1}{2} & - \frac{1}{2} \twiddle{\pi}{2 n + 1}{N} \\ \end{pmatrix} \begin{pmatrix} \idft{n}{N / 2}{Z_0}{Z_2}{Z_{N - 2}} \\ \idft{n}{N / 2}{Z_1}{Z_3}{Z_{N - 1}} \\ \end{pmatrix}\end{split}\]

for \(\seq{n}{0}{1}{N / 2 - 1}\).

src/sdft.c

static int isdft (
    const size_t nitems,
    const size_t stride,
    const trig_t * const trigs,
    const double complex * const xs,
    double complex * const ys
) {
  if (1 == nitems) {
    ys[0] = xs[0];
  } else if (0 == nitems % 2) {
    const size_t nh = nitems / 2;
    // divide-and-conquer
    isdft(nh, stride * 2, trigs, xs         , ys     );
    isdft(nh, stride * 2, trigs, xs + stride, ys + nh);
    // unify two sub problem results
    for (size_t i = 0; i < nh; i++) {
      const size_t j = i + nh;
      const trig_t * const trig = trigs + (2 * i + 1) * stride;
      const double complex twiddle = trig->cos + I * trig->sin;
      double complex * const yi = ys + i;
      double complex * const yj = ys + j;
      const double complex ye = *yi;
      const double complex yo = *yj;
      *yi = 0.5 * ye + 0.5 * twiddle * yo;
      *yj = 0.5 * ye - 0.5 * twiddle * yo;
    }
  } else {
    // naive O(N^2) isdft
    for (size_t n = 0; n < nitems; n++) {
      double complex * const y = ys + n;
      *y = 0. + I * 0.;
      for (size_t k = 0; k < nitems; k++) {
        *y += xs[stride * k] * cexp(2. * pi * (n + 0.5) * k * I / nitems);
      }
      *y /= nitems;
    }
  }
  return 0;
}