lasp/test/scipy_sos_thirdoctavebank.py

"""
This script checks whether the 3rd octave equalizer is norm compliant.

Created on Tue Dec 31 14:43:19 2019
@author: anne
"""

from lasp.filter import (ThirdOctaveBankDesigner)
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import freqz, iirdesign, sosfreqz, butter, cheby2


plt.close('all')


# %% Run check
def check_one(fs):
    """Check for fs = fs"""
    # Create filter bank
    nsections = 5  # filter property
    designer = ThirdOctaveBankDesigner(fs)

    # Test
    xmin = designer.nominal_txt_tox('25')  # lowest highest band, nominal name
    xmax = designer.nominal_txt_tox('20k')  # highest band

    compliant = True  # remains True if all filters comply
    fnyq = fs/2
    fll = designer.fl(xmin) / 1.1  # plot limits
    fuu = designer.fu(xmax) * 1.1
    fig = plt.figure()
    plt.xlabel('f (Hz)')
    plt.ylabel('Magnitude (dB)')
    for x in range(xmin, xmax+1):  # loop over bands

        fl = designer.fl(x)  # lower cut off
        fu = designer.fu(x)  # upper

        fln = fl/fnyq  # lower cutoff, normalized [0, 1]
        fun = fu/fnyq  # upper

        sos = designer.createSOSFilter(x)
        fnorm, h = sosfreqz(sos, worN=2**18, whole=False)
        freq = fnorm*fnyq/np.pi

        h_dB = 20*np.log10(np.abs(h) + np.finfo(float).eps)
        if not designer.testStandardCompliance(x, freq, h_dB):
            print('==============================')
            print(designer.nominal_txt(x), ' does not comply\n')
            compliant = False

        plt.semilogx(freq, h_dB)

        # Upper / lower limits
        freqlim, ulim, llim = designer.band_limits(x)
        plt.semilogx(freqlim, ulim)
        plt.semilogx(freqlim, llim)

        plt.ylim(-5, 1)
        plt.xlim(fll, fuu)

    print(f"Passed test : {compliant} for fs = {fs} Hz")
    return compliant


if __name__ == "__main__":
    # Make list of fs to be tested
    # Taken from https://en.wikipedia.org/wiki/Sampling_(signal_processing)
    fs_all = [8000, 11025, 16000, 22050, 32000, 37800, 44056, 44100, 47250,
              48000, 50000, 50400, 64000, 88200, 96000, 176400, 192000, 352800]
    allcompliant = True
    for fs in fs_all:
        compliant = check_one(fs)
        if not compliant:
            allcompliant = False
print(f"\n\nAll passed test : {allcompliant}")