Sliding DFT analysis in python (transcribed simulink model)

sliding_dft.py

import numpy
from matplotlib import pyplot
import math
import random
from collections import deque


class DelayLine:
    def __init__(self, length):
        self.buffer = deque()
        self.length = length

    def next(self, value):
        self.buffer.append(value)
        if len(self.buffer) == self.length:
            return self.buffer.popleft()

        return 0


class SubSystem:
    def __init__(self):
        self.input_delay_line = DelayLine(32)
        self.prev_output = 0

    def next(self, value):
        delayed_value = self.input_delay_line.next(value)
        output = value - delayed_value + self.prev_output
        self.prev_output = output
        return output


class SDFT:
    def __init__(self):
        self.sub_system = SubSystem()
        self.multiplier_delay = DelayLine(3)

    def __next(self, dac, feedback):
        value = self.multiplier_delay.next(dac * feedback)
        return self.sub_system.next(value)

    def compute(self, excitation_values, feedback_values, scale):
        values = []
        for i in range(0, len(excitation_values)):
            values.append(self.__next(excitation_values[i], feedback_values[i]) * scale)

        return values


def generate_wave(clock_frequency, signal_frequency, num_cycles, fn, latency=0, scale=1.0):
    print("latency {}".format(latency))
    total_time = num_cycles / signal_frequency
    clock_cycles = int(clock_frequency * total_time)
    clock_timestep = 1.0 / clock_frequency
    times = [c * clock_timestep for c in range(0, clock_cycles)]

    data = [0] * clock_cycles
    for c, t in enumerate(times):
        phase = (t + latency) * signal_frequency * 2.0 * math.pi
        data[c] = fn(phase)

    return times, data


def do_sdft(excitation_sin, excitation_cos, excitation_feedback):
    real_sdft = SDFT()
    real_values = real_sdft.compute(excitation_cos, excitation_feedback, -2.0)

    imag_sdft = SDFT()
    imag_values = imag_sdft.compute(excitation_sin, excitation_feedback, 1.0)

    return real_values, imag_values


def main_sdft():
    clock_frequency = 100000.0
    signal_frequency = 3125.0
    num_cycles = 5

    times, excitation_sin = generate_wave(clock_frequency, signal_frequency, num_cycles, math.sin)
    _, excitation_cos = generate_wave(clock_frequency, signal_frequency, num_cycles, math.cos)

    us_times = numpy.array(times) * 1.0e6

    def plot_sdft(feedback_latency, show_magnitudes):
        scale = 0.5
        _, excitation_feedback = generate_wave(clock_frequency,
                                               signal_frequency,
                                               num_cycles,
                                               math.sin,
                                               feedback_latency,
                                               scale)

        real_values, imag_values = do_sdft(excitation_sin, excitation_cos, excitation_feedback)
        if show_magnitudes:
            values = []
            for i in range(0, len(real_values)):
                r = real_values[i]
                i = imag_values[i]
                values.append(math.sqrt(r * r + i * i))
            pyplot.ylabel("Magnitude")
        else:
            values = real_values
            pyplot.ylabel("Real part")

        pyplot.plot(us_times, values)
        pyplot.xlabel("Time [us]")

        stable_start_index = 40  # index where signal stabilizes
        mean_val = numpy.mean(values[stable_start_index])
        pyplot.title("Latency {} ns (mean of stable part {})".format(
            feedback_latency * 1.0e9, mean_val))

    pyplot.subplot(2, 1, 1)
    plot_sdft(0.0, False)

    pyplot.subplot(2, 1, 2)
    plot_sdft(1.0 / 50.0e6, False)

    pyplot.tight_layout()
    pyplot.show()


if __name__ == "__main__":
    main_sdft()