Initial implementation of polyphase quadrature filterbank

1 files changed, 158 insertions, 0 deletions

diff --git a/tools/plot/main.py b/tools/plot/main.py
new file mode 100644
index 0000000..ddf0bce
--- /dev/null
+++ b/tools/plot/main.py
@@ -0,0 +1,158 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import argparse
+import sys
+import os
+import re
+
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.fftpack
+
+import pqf
+
+RUNS = 2
+SAMPLERATE = 44100
+
+def show_proto(proto):
+
+    fig = plt.figure()
+
+    yf = scipy.fftpack.fft(proto)
+
+    max_freq = SAMPLERATE // 2;
+    f_len = len(proto) // 2
+    freq_step = max_freq / f_len;
+
+    freq = np.empty(f_len)
+
+    for i in range(f_len):
+        freq[i] = i * freq_step
+
+    fig.add_subplot(111).plot(freq, 20 * np.log10(np.abs(yf[:f_len])))
+    fig.add_subplot(222).plot(proto)
+
+    plt.show()
+
+def calc_energy_db(data):
+    res = 0.0
+    for x in data:
+        res += x * x
+
+    return 10 * np.log10(res / len(data))
+
+def do_an(ctx, time, data, shift):
+    data_slice = data[shift : ctx.frame_sz() + shift]
+    time_slice = time[shift : ctx.frame_sz() + shift]
+
+    energy_in = calc_energy_db(data_slice)
+
+    res = ctx.do(data_slice)
+
+    subbands_num = ctx.subbands_num();
+    subband_sz = ctx.subband_sz();
+
+    energy_out = calc_energy_db(res)
+
+    rel = energy_out - energy_in
+
+    attenuations = []
+    for i in range(subbands_num):
+        start = i * subband_sz
+        end = start + subband_sz
+        tmp = calc_energy_db(res[start:end]) - energy_in
+        for j in range(subband_sz):
+            attenuations.append(tmp)
+
+    textstr = 'out - in : %.3f' % rel
+
+    fig = plt.figure()
+
+    color = 'tab:blue'
+    a = fig.add_subplot(111)
+    a.set_ylabel('value', color=color)
+    a.plot(res, color=color)
+    a.set_title("output")
+    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+    a.text(0.01, 0.95, textstr, transform=a.transAxes, fontsize=10,
+        verticalalignment='top', bbox=props)
+
+    ax2 = a.twinx()
+    color = 'tab:red'
+    ax2.set_ylabel('attenuation, db', color=color)
+    ax2.plot(attenuations, color=color)
+
+    b = fig.add_subplot(222)
+
+    textstr = 'in: %.3f' % energy_in
+    b.plot(time_slice, data_slice)
+    b.text(0.5, 0.95, textstr, transform=a.transAxes, fontsize=10,
+        verticalalignment='top', bbox=props)
+
+    b.set_title("input");
+
+    plt.show()
+
+def process_freq(freq, subband_sz, subbands_num, prototype):
+    ctx = pqf.AnalyzeCtx(subband_sz, subbands_num, prototype)
+    frame_sz = ctx.frame_sz()
+    time = np.arange(0, frame_sz * RUNS, 1)
+    amplitude = 1 * np.sin(2 * np.pi * freq * time / SAMPLERATE).astype(np.float32)
+
+    do_an(ctx, time, amplitude, 0)
+    do_an(ctx, time, amplitude, frame_sz)
+
+def show_response(pos, subband_sz, subbands_num, prototype):
+    ctx = pqf.AnalyzeCtx(subband_sz, subbands_num, prototype)
+
+    time = np.arange(0, ctx.frame_sz() * RUNS, 1)
+    amplitude = np.zeros(ctx.frame_sz() * RUNS)
+    amplitude[ctx.frame_sz() + pos] = 1.0
+
+    do_an(ctx, time, amplitude, 0)
+    do_an(ctx, time, amplitude, ctx.frame_sz())
+
+def load_prototype(path):
+    prototype = []
+    with open(path) as fp:
+        for line in fp:
+            if line[0] == "#":
+                continue
+            for val in re.findall(r"[-+]?\d*\.\d+|\d+",line):
+                prototype.append(np.float32(val))
+
+    return prototype
+
+def main():
+    parser = argparse.ArgumentParser()
+
+    modes = ['freq', 'delta', 'proto']
+    parser.add_argument('--mode', choices=modes, required=True, help='run mode')
+    parser.add_argument('--file', type=str, help='path to prototype file', required=True, action='store')
+    parser.add_argument('--subband-size', type=int, help='size of subband', action='store')
+    parser.add_argument('--subbands-num', type=int, help='number of subbands to split', action='store')
+
+    args, free_args = parser.parse_known_args()
+
+    if not os.path.isfile(args.file):
+        print("Prototype file {} does not exist.".format(args.file), file=sys.stderr)
+        sys.exit()
+
+    prototype = load_prototype(args.file)
+
+    if (args.mode == 'proto'):
+        show_proto(prototype)
+    if (args.mode == 'delta'):
+        if len(free_args) != 1:
+            print('one free argument (position) must be provided in delta mode')
+            sys.exit()
+        show_response(int(free_args[0]), args.subband_size, args.subbands_num, prototype)
+    else:
+        if len(free_args) != 1:
+            print('one free argument (frequency) must be provided in frequency mode')
+            sys.exit()
+        process_freq(int(free_args[0]), args.subband_size, args.subbands_num, prototype)
+
+if __name__ == "__main__":
+    main()