A Fast Guaranteed-Stable Sliding DFT Algorithm

fast_sdft_response.py

# SDFT network transfer function example according to
# "A Fast Guaranteed-Stable Sliding DFT Algorithm"
# by Richard Lyons
# https://www.dsprelated.com/showarticle/1533.php


import matplotlib.pyplot as plot
import numpy as np

from scipy.fft import fftshift
from scipy.signal import freqz


def fast_sdft_response(N, k, show=False):

    I = 1 if k == int(k) else -1
    P = np.exp(2j * np.pi * k / N)

    B = [P, -1] + [0] * (N - 2) + [-P * I, I]
    A = [1, -2 * np.real(P), 1]

    with np.errstate(divide='ignore', invalid='ignore'):
        _, H = freqz(B, A, whole=True)

    H[int(k * len(H) / N)] = np.nan

    K = np.linspace(-N / 2, N / 2, len(H))
    H = fftshift(H)

    if show:

        figure, axis1 = plot.subplots()
        axis2 = axis1.twinx()

        axis1.plot(K, np.abs(H), color='b')
        axis1.axvline(k, color='m', linestyle='--')
        axis2.plot(K, np.angle(H), color='g', alpha=0.5)

        axis1.set_ylim(0 - 0.5, N + 0.5)
        axis1.set_ylabel('$|H(z)|$', color='b')
        axis1.set_xlabel('k', color='m')

        axis2.set_ylim(-np.pi - 0.5, +np.pi + 0.5)
        axis2.set_ylabel('$\\angle H(z)$', color='g')
        axis2.set_xlabel('k', color='m')

        axis1.set_title(f'N={N}, k={k}')
        axis1.grid(axis='x')
        axis2.grid(axis='y')

    return K, H


if __name__ == '__main__':

    fast_sdft_response(8, 0,   show=True)
    fast_sdft_response(8, 0.5, show=True)
    plot.show()