#!/usr/bin/env python """ MPAI CAE-ARP Tape Audio Restoration. Implements MPAI CAE-ARP Tape Audio Restoration Technical Specification. It identifies and restore portions of the Preservation Audio File, providing: - Restored Audio Files; - Editing List """ import matplotlib.pyplot as plt import numpy as np import os import shutil import yaml from control import c2d, TransferFunction from scipy.io import wavfile from scipy.signal import freqs, freqz, tf2zpk, zpk2tf, lfilter __author__ = "Nadir Dalla Pozza" __copyright__ = "Copyright 2022, Audio Innova S.r.l." __credits__ = ["Niccolò Pretto", "Nadir Dalla Pozza", "Sergio Canazza"] __license__ = "GPL v3.0" __version__ = "1.0.1" __maintainer__ = "Nadir Dalla Pozza" __email__ = "nadir.dallapozza@unipd.it" __status__ = "Production" # Class for customizing Console Colors class CC: PURPLE = '\033[95m' CYAN = '\033[96m' DARK_CYAN = '\033[36m' BLUE = '\033[94m' GREEN = '\033[92m' YELLOW = '\033[93m' RED = '\033[91m' BOLD = '\033[1m' UNDERLINE = '\033[4m' END = '\033[0m' def save_file(file): if not os.path.exists(temp_path): # Create directory os.mkdir(temp_path) print("temp directory '% s' created" % temp_path) raf_path = os.path.join(temp_path, 'RestoredAudioFiles') make_raf = False if not os.path.exists(raf_path): # Create directory os.mkdir(raf_path) make_raf = True print("Restored Audio Files directory '% s' created" % raf_path) else: print((CC.PURPLE + "Restored Audio Files directory '% s' already exists!" + CC.END) % raf_path) overwrite = input('Do you want to overwrite it? [y/n]: ') if overwrite.casefold() == 'y': # Overwrite directory shutil.rmtree(raf_path) os.mkdir(raf_path) make_raf = True print('Restored Audio Files directory overwritten') elif overwrite.casefold() != 'n': print(CC.RED + 'Unknown command, exiting' + CC.END) quit(os.EX_USAGE) if make_raf: print("Saving Restored Audio File to: '%s' ..." % raf_path) wavfile.write(os.path.join(raf_path, '1.wav'), FS, file) if __name__ == '__main__': # Read configuration file config = object try: config = yaml.safe_load(open('config.yaml', 'r')) if 'WORKING_PATH' not in config: print(CC.RED + 'WORKING_PATH key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) if 'PRESERVATION_FILE_NAME' not in config: print(CC.RED + 'PRESERVATION_FILE_NAME key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) if 'STANDARD_W' not in config: print(CC.RED + 'STANDARD_W key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) if 'SPEED_W' not in config: print(CC.RED + 'SPEED_W key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) if 'STANDARD_R' not in config: print(CC.RED + 'STANDARD_R key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) if 'SPEED_R' not in config: print(CC.RED + 'SPEED_R key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) if 'FS' not in config: print(CC.RED + 'FS key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) if 'BPS' not in config: print(CC.RED + 'BPS key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) if 'PLOTS' not in config: print(CC.RED + 'PLOTS key not found in config.yaml!' + CC.END) quit(os.EX_CONFIG) except FileNotFoundError: print(CC.RED + 'config.yaml file not found!' + CC.END) quit(os.EX_NOINPUT) # config.yaml variables WORKING_PATH = config['WORKING_PATH'] PRESERVATION_FILE_NAME = config['PRESERVATION_FILE_NAME'] STANDARD_W = config['STANDARD_W'] SPEED_W = config['SPEED_W'] STANDARD_R = config['STANDARD_R'] SPEED_R = config['SPEED_R'] FS = config['FS'] BPS = config['BPS'] PLOTS = config['PLOTS'] # Output path temp_path = os.path.join(WORKING_PATH, 'temp', PRESERVATION_FILE_NAME) # Configuration parameters check # Recording tape speed check if SPEED_W != 3.75 and SPEED_W != 7.5 and SPEED_W != 15 and SPEED_W != 30: print( CC.RED + 'Incorrect SPEED_W: \'' + str(SPEED_W) + '\'. Accepted value are: 3.75, 7.5, 15, 30.' + CC.END ) quit(os.EX_CONFIG) # Reading tape speed check. if SPEED_R != 3.75 and SPEED_R != 7.5 and SPEED_R != 15 and SPEED_R != 30: print( CC.RED + 'Incorrect SPEED_R: \'' + str(SPEED_R) + '\'. Accepted value are: 3.75, 7.5, 15, 30.' + CC.END ) quit(os.EX_CONFIG) # Equalization standard check. if not (STANDARD_R == 'CCIR' or STANDARD_R == 'NAB'): print( CC.RED + 'Incorrect STANDARD_R: \'' + STANDARD_R + '\'. Accepted values are: CCIR, NAB.' + CC.END ) quit(os.EX_CONFIG) if not (STANDARD_W == 'CCIR' or STANDARD_W == 'NAB'): print( CC.RED + 'Incorrect STANDARD_W: \'' + STANDARD_W + '\'. Accepted values are: CCIR, NAB.' + CC.END ) quit(os.EX_CONFIG) # CCIR speed check. if STANDARD_W == 'CCIR' and SPEED_W == 3.75: print( CC.YELLOW + 'CCIR is undefined at 3.75 ips. Recording equalization standard is set to NAB.' + CC.END ) STANDARD_W = 'NAB' if STANDARD_R == 'CCIR' and SPEED_R == 3.75: print( CC.YELLOW + 'CCIR is undefined at 3.75 ips. Reading equalization standard is set to NAB.' + CC.END ) STANDARD_R = 'NAB' # NAB speed check. if STANDARD_W == 'NAB' and SPEED_W == 30: print( CC.YELLOW + 'NAB is undefined at 30 ips. Recording equalization standard is set to CCIR.' + CC.END ) STANDARD_W = 'CCIR' if STANDARD_R == 'NAB' and SPEED_R == 30: print( CC.YELLOW + 'NAB is undefined at 30 ips. Reading equalization standard is set to CCIR.' + CC.END ) STANDARD_R = 'CCIR' # Sampling frequency check. if len(PRESERVATION_FILE_NAME) == 0: if FS != 44100 and FS != 48000 and FS != 96000: print( CC.RED + 'Incorrect FS: \'' + str(FS) + '\'. Accepted values are: 44100, 48000, 96000.' + CC.END ) quit(os.EX_CONFIG) # Bits per sample check. if BPS != 8 and BPS != 16 and BPS != 24 and BPS != 32: print( CC.RED + 'Incorrect BPS: \'' + str(BPS) + '\'. Accepted values are: 8, 16, 24, 32.' + CC.END ) quit(os.EX_CONFIG) # Display input parameters print('Input parameters:') print(' WORKING_PATH: ' + WORKING_PATH) print(' PRESERVATION_FILE_NAME: ' + PRESERVATION_FILE_NAME) print(' STANDARD_W: ' + STANDARD_W) print(' SPEED_W: ' + str(SPEED_W) + ' ips') print(' STANDARD_R: ' + STANDARD_R) print(' SPEED_R: ' + str(SPEED_R) + ' ips') if len(PRESERVATION_FILE_NAME) == 0: print(' FS: ' + str(FS) + ' Hz') print(' BPS: ' + str(BPS) + '\n') # Preservation Audio File check paf = [] if len(PRESERVATION_FILE_NAME) > 0: audio_file = PRESERVATION_FILE_NAME + '.wav' paf_path = os.path.join(WORKING_PATH, 'PreservationAudioFile', audio_file) try: print("Opening '%s'..." % paf_path) FS, paf = wavfile.read(paf_path) print('Preservation Audio File opened!') print('Overwritten parameters:') print(' FS: ' + str(FS) + ' Hz\n') except OSError: print(CC.RED + "Preservation Audio File not found!" + CC.END) quit(os.EX_NOINPUT) # Equalization standard time constants # CCIR time constants. t2_30 = 17.5 * 10**(-6) # time constant CCIR_30 t2_15 = 35 * 10**(-6) # time constant CCIR_15 t2_7 = 70 * 10**(-6) # time constant CCIR_7.5 # NAB time constants. t3 = 3180 * 10**(-6) t4_15 = 50 * 10**(-6) # time constant NAB_15 t4_7 = 50 * 10**(-6) # time constant NAB_7.5 t4_3 = 90 * 10**(-6) # time constant NAB_3.75 # Decision stage a = [] # Filter numerator b = [] # Filter denominator case = 0 # Reference case number # This section will establish which time constants must be modified to obtain the desired equalisation standard. if STANDARD_W == 'CCIR': if SPEED_W == 30: if STANDARD_R == 'NAB': # Case 1 if SPEED_R == 15: FS = 2 * FS # Doubling the sampling frequency # Correction filter: NABw15_mod + CCIRr30 # - NAB constants divided by 2 t3 = t3 / 2 t4 = t4_15 / 2 # - CCIR_30 constant not altered t2 = t2_30 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 1 # Case 2 elif SPEED_R == 7.5: FS = 4 * FS # Quadrupling the sampling frequency # Correction filter: NABw7.5_mod + CCIRr30 # - NAB constants divided by 4 t3 = t3 / 4 t4 = t4_7 / 4 # - CCIR_30 constant not altered t2 = t2_30 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 2 # Case 3 else: # SPEED_R == 3.75 FS = 8 * FS # Multiplying by 8 the sampling frequency # Correction filter: NABw3.75_mod + CCIRr30 # - NAB constants divided by 8 t3 = t3 / 8 t4 = t4_3 / 8 # - CCIR_30 constant not altered t2 = t2_30 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 3 else: # STANDARD_R == 'CCIR' # Case 31 if SPEED_R == 30: print('Reference case: 31') print(CC.GREEN + 'Nothing to do!' + CC.END) quit(os.EX_OK) # Case 15 elif SPEED_R == 15: FS = 2 * FS # Doubling sampling frequency # Plot information case = 15 # Case 16 else: # SPEED_R == 7.5 FS = 4 * FS # Quadrupling the sampling frequency # Plot information case = 16 elif SPEED_W == 15: if STANDARD_R == 'NAB': # Case 28 if SPEED_R == 15: # No speed change # Correction filter: NABw15 + CCIRr15 # - NAB_15 constants not altered t4 = t4_15 # - CCIR_30 constant not altered t2 = t2_15 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 28 # Case 6 elif SPEED_R == 7.5: FS = 2 * FS # Doubling the sampling frequency # Correction filter: NABw7.5_mod + CCIRr15 # - NAB constants divided by 2 t3 = t3 / 2 t4 = t4_7 / 2 # - CCIR_15 constant not altered t2 = t2_15 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 6 # Case 7 else: # SPEED_R == 3.75 FS = 4 * FS # Quadrupling the sampling frequency # Correction filter: NABw3.75_mod + CCIRr15 # - NAB constants divided by 4 t3 = t3 / 4 t4 = t4_3 / 4 # - CCIR_15 constant not altered t2 = t2_15 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 7 else: # STANDARD_R == 'CCIR' # Case 19 if SPEED_R == 30: FS = FS / 2 # Halving the sampling frequency # Plot information case = 19 # Case 33 elif SPEED_R == 15: print('Reference case: 33') print(CC.GREEN + 'Nothing to do!' + CC.END) quit(os.EX_OK) # Case 20 else: # SPEED_R == 7.5 FS = 2 * FS # Doubling the sampling frequency # Plot information case = 20 else: # SPEED_W == 7.5 if STANDARD_R == 'NAB': # Case 10 if SPEED_R == 15: FS = FS / 2 # Halving the sampling frequency # Correction filter: NABw15_mod + CCIRr7.5 # - NAB constants multiplied by 2 t3 = t3 * 2 t4 = t4_15 * 2 # - CCIR_7.5 constant not altered t2 = t2_7 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 10 # Case 30 elif SPEED_R == 7.5: # No speed change # Correction filter: NABw7.5 + CCIRr7.5 # - NAB_7.5 constant not altered t4 = t4_7 # - CCIR_7.5 constant not altered t2 = t2_7 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 30 # Case 11 else: # SPEED_R == 3.75 FS = 2 * FS # Doubling the sampling frequency # Correction filter: NABw3.75_mod + CCIRr7.5 # - NAB constants divided by 2 t3 = t3 / 2 t4 = t4_3 / 2 # - CCIR_7.5 constant not altered t2 = t2_7 # Filter coefficients a = [t2 * t3, t2 + t3, 1] b = [t3 * t4, t3, 0] # Plot information case = 11 else: # STANDARD_R == 'CCIR' # Case 23 if SPEED_R == 30: FS = FS / 4 # Quartering the sampling frequency # Plot information case = 23 # Case 24 elif SPEED_R == 15: FS = FS / 2 # Halving the sampling frequency # Plot information case = 24 # Case 35 else: # SPEED_R == 7.5 print('Reference case: 35') print(CC.GREEN + 'Nothing to do!' + CC.END) quit(os.EX_OK) else: # STANDARD_W == 'NAB' if SPEED_W == 15: if STANDARD_R == 'NAB': # Case 32 if SPEED_R == 15: print('Reference case: 32') print(CC.GREEN + 'Nothing to do!' + CC.END) quit(os.EX_OK) # Case 17 elif SPEED_R == 7.5: FS = 2 * FS # Doubling the sampling frequency # Correction filter: NABw7.5_mod + NABr15 # - NABw constants divided by 2 t3_mod = t3 / 2 t4_mod = t4_7 / 2 # - NABr constant not altered t4 = t4_15 # Filter coefficients a = [t3 * t3_mod * t4, t3 * (t3_mod + t4), t3] b = [t3 * t3_mod * t4_mod, t3_mod * (t3 + t4_mod), t3_mod] # Plot information case = 17 # Case 18 else: # SPEED_R == 3.75 FS = 4 * FS # Quadrupling the sampling frequency # Correction filter: NABw3.75_mod + NABr15 # - NAB constants divided by 4 t3_mod = t3 / 4 t4_mod = t4_3 / 4 # - NABr constant not altered t4 = t4_15 # Filter coefficients a = [t3 * t3_mod * t4, t3 * (t3_mod + t4), t3] b = [t3 * t3_mod * t4_mod, t3_mod * (t3 + t4_mod), t3_mod] # Plot information case = 18 else: # STANDARD_R == 'CCIR' # Case 4 if SPEED_R == 30: FS = FS / 2 # Halving the sampling frequency # Correction filter: CCIRw30_mod + NABr15 # - CCIR_30 constant multiplied by 2 t2 = t2_30 * 2 # - NAB_15 constant not altered t4 = t4_15 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 4 # Case 27 elif SPEED_R == 15: # No speed change # Correction filter: CCIRw15 + NABr15 # - CCIR_15 constant not altered t2 = t2_15 # - NAB_15 constant not altered t4 = t4_15 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 27 # Case 5 else: # SPEED_R == 7.5 FS = FS * 2 # Doubling the sampling frequency # Correction filter: CCIRw7.5_mod + NABr15 # - CCIR_7.5 constant divided by 2 t2 = t2_7 / 2 # - NAB_15 constant not altered t4 = t4_15 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 5 elif SPEED_W == 7.5: if STANDARD_R == 'NAB': # Case 21 if SPEED_R == 15: FS = FS / 2 # Halving the sampling frequency # Correction filter: NABw15_mod + NABr7.5 # - NABw constants multiplied by 2 t3_mod = t3 * 2 t4_mod = t4_15 * 2 # - NABr constant not altered t4 = t4_7 # Filter coefficients a = [t3 * t3_mod * t4, t3 * (t3_mod + t4), t3] b = [t3 * t3_mod * t4_mod, t3_mod * (t3 + t4_mod), t3_mod] # Plot information case = 21 # Case 34 elif SPEED_R == 7.5: print('Reference case: 34') print(CC.GREEN + 'Nothing to do!' + CC.END) quit(os.EX_OK) # Case 22 else: # SPEED_R == 3.75 FS = 2 * FS # Doubling the sampling frequency # Correction filter: NABw3.75_mod + NABr7.5 # - NABw constants divided by 2 t3_mod = t3 / 2 t4_mod = t4_3 / 2 # - NABr constant not altered t4 = t4_7 # Filter coefficients a = [t3 * t3_mod * t4, t3 * (t3_mod + t4), t3] b = [t3 * t3_mod * t4_mod, t3_mod * (t3 + t4_mod), t3_mod] # Plot information case = 22 else: # STANDARD_R == 'CCIR' # Case 8 if SPEED_R == 30: FS = FS / 4 # Quartering the sampling frequency # Correction filter: CCIRw30_mod + NABr7.5 # - CCIR_30 constant multiplied by 4 t2 = t2_30 * 4 # - NAB_7.5 constant not altered t4 = t4_7 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 8 # Case 9 elif SPEED_R == 15: FS = FS / 2 # Halving the sampling frequency # Correction filter: CCIRw15_mod + NABr7.5 # - CCIR_15 constant multiplied by 2 t2 = t2_15 * 2 # - NAB_7.5 constant not altered t4 = t4_7 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 9 # Case 29 else: # SPEED_R == 7.5 # No speed change # Correction filter: CCIRw7.5 + NABr7.5 # - CCIR_7.5 constant not altered t2 = t2_7 # - NAB_7.5 constant not altered t4 = t4_7 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 29 else: # SPEED_W == 3.75 if STANDARD_R == 'NAB': # Case 25 if SPEED_R == 15: FS = FS / 4 # Quartering the sampling frequency # Correction filter: NABw15_mod + NABr3.75 # - NAB constants multiplied by 4 t3_mod = t3 * 4 t4_mod = t4_15 * 4 # - NABr constant not altered t4 = t4_3 # Filter coefficients a = [t3 * t3_mod * t4, t3 * (t3_mod + t4), t3] b = [t3 * t3_mod * t4_mod, t3_mod * (t3 + t4_mod), t3_mod] # Plot information case = 25 # Case 26 elif SPEED_R == 7.5: FS = FS / 2 # Halving the sampling frequency # Correction filter: NABw7.5_mod + NABr3.75 # - NAB constants multiplied by 2 t3_mod = t3 * 2 t4_mod = t4_7 * 2 # - NABr constant not altered t4 = t4_3 # Filter coefficients a = [t3 * t3_mod * t4, t3 * (t3_mod + t4), t3] b = [t3 * t3_mod * t4_mod, t3_mod * (t3 + t4_mod), t3_mod] # Plot information case = 26 # Case 36 else: # SPEED_R == 3.75 print('Reference case: 36') print(CC.GREEN + 'Nothing to do!' + CC.END) quit(os.EX_OK) else: # STANDARD_R == 'CCIR' # Case 12 if SPEED_R == 30: FS = FS / 8 # Dividing by 8 the sampling frequency # Correction filter: CCIRw30_mod + NABr3.75 # - CCIR_30 constant multiplied by 8 t2 = t2_30 * 8 # - NAB_3.75 constant not altered t4 = t4_3 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 12 # Case 13 elif SPEED_R == 15: FS = FS / 4 # Quartering the sampling frequency # Correction filter: CCIRw15_mod + NABr3.75 # - CCIR_15 constant multiplied by 4 t2 = t2_15 * 4 # - NAB_3.75 constant not altered t4 = t4_3 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 13 # Case 14 else: # SPEED_R == 7.5 FS = FS / 2 # Halving the sampling frequency # Correction filter: CCIRw7.5_mod + NABr3.75 # - CCIR_7.5 constant multiplied by 2 t2 = t2_7 * 2 # - NAB_3.75 constant not altered t4 = t4_3 # Filter coefficients a = [t3 * t4, t3, 0] b = [t2 * t3, t2 + t3, 1] # Plot information case = 14 # Casting FS to int because wavfile.write() is stupid FS = round(FS) print('Reference case: ' + str(case)) print('Operational sampling frequency: ' + str(FS) + ' Hz.') # Correction filter # Not all cases present a correction filter! if len(a) != 0: # Analog transfer function H_a = TransferFunction(a, b) # Analog frequency vector w_a = np.logspace(np.log10(1), np.log10(FS * np.pi), 5000) if PLOTS: # Analog filter frequency response w_t, h_t = freqs(a, b, worN=w_a) # Plot analog graph # - Magnitude plt.subplot(2, 1, 1) plt.semilogx(w_t / (2 * np.pi), 20 * np.log10(abs(h_t))) plt.xlim([1, 24000]) plt.xlabel('Frequency') plt.ylim([-40, 40]) plt.ylabel('Amplitude response [dB]') plt.grid(True) # - Phase plt.subplot(2, 1, 2) plt.semilogx(w_t / (2 * np.pi), np.angle(h_t) * 180 / np.pi) plt.xlim([1, 24000]) plt.xlabel('Frequency') plt.ylabel('Phase [deg]') plt.grid(True) # Digital transfer function through bilinear digitisation H_d = c2d(H_a, 1/FS, 'bilinear') num_d = H_d.num[0][0] # Inspect Hd.num to see why [0][0] is needed... den_d = H_d.den[0][0] # Same story here # Digital frequency vector w_d = np.logspace(np.log10(1), np.log10(FS / 2), 5000) if PLOTS: # Digital filter frequency response w_n, h_n = freqz(num_d, den_d, worN=w_d, fs=FS) # Plot digital graph # - Magnitude plt.subplot(2, 1, 1) plt.semilogx(w_n, 20 * np.log10(abs(h_n)), '--') plt.legend(['Analog', 'Bilinear']) # - Phase plt.subplot(2, 1, 2) plt.semilogx(w_n, np.angle(h_n) * 180 / np.pi, '--') plt.legend(['Analog', 'Bilinear']) # Pole check # New pole frequency pole_frequency = 2 # Move to zero-pole representation z, p, k = tf2zpk(a, b) # Check if the function presents a pole at 0 Hz for i in range(len(p)): if p[i] == 0: # Replace pole p[i] = -pole_frequency * 2 * np.pi print('\n' + CC.PURPLE + 'Pole at 0 Hz replaced!' + CC.END) # Back to transfer function representation ap, bp = zpk2tf(z, p, k) # Analog transfer function Hp_a = TransferFunction(ap, bp) if PLOTS: # Analog filter frequency response wp_t, hp_t = freqs(ap, bp, worN=w_a) # Plot analog graph # - Magnitude plt.subplot(2, 1, 1) plt.semilogx(wp_t / (2 * np.pi), 20 * np.log10(abs(hp_t))) # - Phase plt.subplot(2, 1, 2) plt.semilogx(wp_t / (2 * np.pi), np.angle(hp_t) * 180 / np.pi) # Digital transfer function through bilinear digitisation Hp_d = c2d(Hp_a, 1 / FS, 'bilinear') num_d = Hp_d.num[0][0] den_d = Hp_d.den[0][0] if PLOTS: # Digital filter frequency response wp_n, hp_n = freqz(num_d, den_d, worN=w_d, fs=FS) # Plot digital graph # - Magnitude plt.subplot(2, 1, 1) plt.semilogx(wp_n, 20 * np.log10(abs(hp_n)), '--') plt.legend(['Analog', 'Bilinear', 'Pole - Analog', 'Pole - Digital']) # - Phase plt.subplot(2, 1, 2) plt.semilogx(wp_n, np.angle(hp_n) * 180 / np.pi, '--') plt.legend(['Analog', 'Bilinear', 'Pole - Analog', 'Pole - Digital']) if PLOTS: plt.show() if len(PRESERVATION_FILE_NAME) > 0: print('\nFiltering Preservation Audio File...') # Filter Preservation Audio File raf = lfilter(num_d, den_d, paf, axis=0) # Again, wavfile.write() is stupid, and you must cast everything to not destroy your ears... raf = np.rint(raf).astype(paf.dtype) # Save Restored Audio File save_file(raf) else: # Save Restored Audio File save_file(paf) # End print(CC.GREEN + CC.BOLD + "Success!" + CC.END + '\n')