Tinnitus_Musik_Therapie_Projektarbeit_Repo/TinnitusAnalyse/DigitalFilterTest.py

import matplotlib.pyplot as plt  # For plotting
from math import sin, pi, cos  # For generating input signals
import numpy as np
import sys  # For reading command line arguments


fs = 44100  # sampling frequency (Abtastfrequenz)

koeff = {
    "b0": 0,
    "b1": 0,
    "b2": 0,
    "a1": 0,
    "a2": 0
}


def koeffizienten_berechnen(omega, r):
    # Koeffizientenberechnung nach Tobola VL S.107 - IIR Filter, 2. Ordnung
    koeff["b0"] = 1
    koeff["b1"] = 0
    koeff["b2"] = -1
    koeff["a1"] = -2*r*cos(omega)
    koeff["a2"] = r**2


def filter(x):
    y = [0]*len(x)
    for k in range(4, len(x)):
        y[k] = koeff["b0"]*x[k] + koeff["b1"]*x[k-1] + koeff["b2"]*x[k-2] - koeff["a1"]*y[k-1] - koeff["a2"]*y[k-2]

    return y




dauer_ms = 10000  # 10 Sekunden
num_samples = dauer_ms * (fs / 1000)  # framerate -pro Sekunde- umgerechnet in -pro Millisekunde-

t = np.linspace(0, 10, int(num_samples)) # array zum darstellen der x-Achse
amp = 1

f_input1 = 1
f_input2 = 50
f_input3 = 10

input1 = []
input2 = []
input3 = []
for x in range(int(num_samples)):  # einen einfachen Sinus ins array schreiben
    input1.append(amp * sin(2 * pi * f_input1 * (x / fs)))
    input2.append(amp * sin(2 * pi * f_input2 * (x / fs)))
    #input3.append(amp * sin(2 * pi * f_input3 * (x / fs)))


input_ges = np.add(input1, input2) # Sinus aufaddieren um ein halbwegs realistisches Audiosignal zu bekommen
#input_ges = np.add(input_ges, input3)

#Filterparameter hier einstellen
fr = 5
omega = 2*pi*fr/fs
r = 0.95
koeffizienten_berechnen(omega, r)   # Koeffizienten berechnen mit Mittelfrequenz 10Hz

output = filter(input_ges)

### Plot the signals for comparison
plt.figure(1)
plt.subplot(231)
plt.ylabel('Amplitude')
plt.xlabel('t [s]')
plt.title('f_1=' + str(f_input1) + "Hz")
plt.plot(t, input1)

plt.subplot(232)
plt.ylabel('Amplitude')
plt.xlabel('t [s]')
plt.title('f_2=' + str(f_input2) + "Hz")
plt.plot(t, input2)

# plt.subplot(233)
# plt.ylabel('Amplitude')
# plt.xlabel('t [s]')
# plt.title('f_3=' + str(f_input3) + "Hz")
# plt.plot(t, input3)

plt.subplot(234)
plt.ylabel('Amplitude')
plt.xlabel('t [s]')
plt.title('(input in Filter) f_ges = f_1 + f_2')
plt.plot(t, input_ges)

plt.subplot(235)
plt.ylabel('Amplitude')
plt.xlabel('Samples')
plt.title('gefiltertes Signal, mit f_r = ' + str(fr) + "Hz und r=" + str(r))
plt.plot(t, output)
plt.show()