Tinnitus_Musik_Therapie_Projektarbeit_Repo/TinnitusAnalyse/SoundGenerator.py

import wave  # bearbeiten von .wav-Dateien
import struct
import sounddevice as sd
import numpy as np
import sys  # für Fehlermeldungen
from scipy import signal
import csv
import time
import matplotlib.pyplot as plt

"""---------------------------------------------------------------------------------------------------------------------
In .wav-Dateien wird der Ton in absoluten Werte eingetragen. Die Standart-framerate ist 44100
das heißt für jede Sekunde an Ton gibt es 44100 Werte, die die Tonwelle über die Zeit beschreiben
---------------------------------------------------------------------------------------------------------------------"""


class Tinnitus:  # beinhaltet alle Werte, die vom Nutzer eingestellt werden
    def __init__(self, l_freq=0, r_freq=0, l_amp=0, r_amp=0, l_rausch=0, r_rausch=0, ohr=0, l_rausch_ug=10, r_rausch_ug=10, l_rausch_og=20000, r_rausch_og=20000):
        self.vorname = ""
        self.nachname = ""
        self.kommentar = ""
        # [Alle Frequenzangaben in Hz, alle Amplitudenangaben von 0-1 (float)]
        self.linksFrequenz = l_freq
        self.rechtsFrequenz = r_freq
        self.linksLautstaerke = l_amp
        self.rechtsLautstaerke = r_amp
        self.linksRauschenLautstaerke = l_rausch
        self.rechtsRauschenLautstaerke = r_rausch
        self.linksRauschenUntereGrenzfrequenz = l_rausch_ug
        self.rechtsRauschenUntereGrenzfrequenz = r_rausch_ug
        self.linksRauschenObereGrenzfrequenz = l_rausch_og
        self.rechtsRauschenObereGrenzfrequenz = r_rausch_og

    def speichern(self):  # speichert die Nutzerdaten in eine .csv-Datei
        datei = open("TinnitusDaten.csv", "w")

        daten = "Vorname;" + self.vorname + "\n"
        daten += "Nachname;" + self.nachname + "\n"
        daten += "linke Frequenz;" + str(self.linksFrequenz) + "\n"
        daten += "linke Lautstärke;" + str(self.linksLautstaerke) + "\n"
        daten += "linkes Rauschen Lautstärke;" + str(self.linksRauschenLautstaerke) + "\n"
        daten += "linkes Rauschen untere Grenzfrequenz;" + str(self.linksRauschenUntereGrenzfrequenz) + "\n"
        daten += "linkes Rauschen obere Grenzfrequenz;" + str(self.linksRauschenObereGrenzfrequenz) + "\n"

        daten += "rechte Frequenz;" + str(self.rechtsFrequenz) + "\n"
        daten += "rechte Lautstärke;" + str(self.rechtsLautstaerke) + "\n"
        daten += "rechtes Rauschen Lautstärke;" + str(self.rechtsRauschenLautstaerke) + "\n"
        daten += "rechtes Rauschen untere Grenzfrequenz;" + str(self.rechtsRauschenUntereGrenzfrequenz) + "\n"
        daten += "rechtes Rauschen obere Grenzfrequenz;" + str(self.rechtsRauschenObereGrenzfrequenz) + "\n"

        daten += "Kommentar;" + str(self.kommentar) + "\n"

        datei.write(daten)
        datei.close()


"""---------------------------------KLASSE: SOUND-----------------------------------------------------------------------
Sound beinhaltet alle Variablen, die zum erstellen einer .wav-Datei benötigt werden (siehe soun.wav_speichern())
   Das 'sound_obj' ist für das dynamische abspielen zuständig (siehe sound.play())
   Beim Initialisieren muss ein Tinnitus-Objekt übergeben werden
---------------------------------------------------------------------------------------------------------------------"""


class Sound:
    def __init__(self, tinnitus, wav_name="MeinTinnitus.wav", audio=None, nchannels=2, sampwidth=2, framerate=44100,
                 comptype="NONE", compname="not compressed", mute=True):
        if audio is None:
            audio = []
        self.tinnitus = tinnitus
        self.wav_name = wav_name #Der Dateiname
        self.audio = audio  # ein Array, in das die Sound-Werte geschrieben werden (von -1, bis +1)
        self.nchannels = nchannels  # Zahl der audio channels (1:mono 2:stereo)
        self.sampwidth = sampwidth  # Größe eines einzelnen Sound-Werts (in bytes)
        self.framerate = framerate  # Abtastrate
        self.nframes = len(audio)  # Anzahl der Sound-Werte -> Muss bei jeder audio-Änderung aktuallisiert werden
        self.comptype = comptype
        self.compname = compname
        self.mute = mute  # wenn der mute boolean auf true gesetzt ist, sollte kein Ton ausgegeben werden
        self.sound_obj = sd.OutputStream(channels=2, callback=self.callback,
                                         samplerate=self.framerate)  # Objekt fürs Abspielen (siehe sound.play())
        self.start_idx = 0  # wird für sound_obj benötigt
        self.music_samplerate = 0     # die samplerate der ausgewählten Musikdatei
        self.music_data = 0           # das Numpy Array der ausgewählten Musikdatei
        self.filterfortschritt = 0, 0.  # für feedback-fkt, 1.Position: Abschnitt-Nmr, 2.Positon: Speicherfortschritt


    def wav_speichern(self):  # ezeugt/aktuallisiert die .wav-Datei
        print("Sound wird als .wav-Datei gespeichert. Bitte warten...\nDer Vorgang kann ca. 10 Sekunden dauern")

        wav_obj = wave.open(self.wav_name, "w")

        # Die Callback-Funktion aufrufen, um die Audiodaten zu bekommen
        frames = self.framerate * 10  # entspricht 10 Sekunden
        status = ""  # für den Funktionsaufruf benötigt, sonst keine Funktion
        audio = np.ones((frames, 2))
        audio = self.callback(audio, frames, self.sound_obj.time, status)

        # Rahmenparameter für die .wav-Datei setzen
        self.nframes = len(audio)
        wav_obj.setparams((self.nchannels, self.sampwidth, self.framerate, self.nframes, self.comptype, self.compname))

        packedMusic = []  # Liste an die wir die einzelnen frames bereits in binär umgewandelt aneinanderreihen

        #Die Audiosamples schreiben
        for x in range(self.nframes):  # geht jeden Sample-Wert der Musikdatei einzeln durch
            # Die Audiodaten müssen von float in einen passenden int-Wert umgerechnet werden
            packedMusic.append(struct.pack('h', int(audio[x][0] * 32767.0)))
            packedMusic.append(struct.pack('h', int(audio[x][1] * 32767.0)))

        value_str = b"".join(packedMusic)
        wav_obj.writeframes(value_str)

        wav_obj.close()
        print("Speichern beendet.")

    """Die Objekt-Funktion 'start()' startet die asynchrone Soundwiedergabe. Sie ruft dabei immer wieder die Funktion
       sound.callback() auf. Daher können die dort genutzten Variablen dynamisch geändert werden. """

    def play(self):
        if not self.mute:  # Nie abspielen, wenn die GUI auf stumm geschaltet ist
            self.sound_obj.start()  # öffnet thread der immer wieder callback funktion aufruft und diese daten abspielt


    def stop(self):
        self.sound_obj.stop()  # beendet die asynchrone Soundwiedergabe

    """Die Funktion callback() erzeugt bei jedem Aufruf die Audiodaten, abhängig von den aktuellen Tinnitus-Variablen.
       Die errechneten  Werte werden in das NumPy-Array 'outdata' geschrieben """

    def callback(self, outdata, frames, time, status):
        if status:  # Warnungen, wenn das Soundobj. auf Fehler stößt (hauptsächlich over/underflow wg. Timingproblemen)
            print(status, file=sys.stderr)


        # Whitenoise erzeugen
        for x in range(len(outdata)):
            rand = (np.random.rand() - 0.5)  # Zufallszahl zwischen -0.5 und 0.5
            # links:
            outdata[x][0] = rand * self.tinnitus.linksRauschenLautstaerke
            # rechts:
            outdata[x][1] = rand * self.tinnitus.rechtsRauschenLautstaerke

        # Whitenoise durch Bandpass laufen lassen
        if self.tinnitus.linksRauschenLautstaerke:
            # (-3dB Grenzen) bzw was der Bandpass durchlässt
            fGrenz = [self.tinnitus.linksRauschenUntereGrenzfrequenz,
                         self.tinnitus.linksRauschenObereGrenzfrequenz]
            # sos (=second order sections = Filter 2. Ordnung) ist ein Array der Länge (filterOrder) und beinhaltet
            # die Koeffizienten der IIR Filter 2. Ordnung (b0, b1, b2 & a0, a1, a2)
            sos = signal.butter(5, fGrenz, 'bandpass', fs=self.framerate, output='sos')
            # sosfilt filtert das Signal mittels mehrerer 'second order sections' (= Filter 2. Ordnung) die über sos definiert sind
            outdata[:, 0] = signal.sosfilt(sos, outdata[:, 0])

            # Plotten des Filters für Filterentwicklung und Dokumentation nützlich---------
            # w, h = signal.sosfreqz(sos, worN=1500)
            # plt.subplot(2, 1, 1)
            # db = 20 * np.log10(np.maximum(np.abs(h), 1e-5))
            # plt.plot(w / np.pi, db)
            # plt.ylim(-75, 5)
            # plt.grid(True)
            # plt.yticks([0, -20, -40, -60])
            # plt.ylabel('Gain [dB]')
            # plt.title('Frequency Response')
            # plt.subplot(2, 1, 2)
            # plt.plot(w / np.pi, np.angle(h))
            # plt.grid(True)
            # plt.yticks([-np.pi, -0.5 * np.pi, 0, 0.5 * np.pi, np.pi],
            #            [r'$-\pi$', r'$-\pi/2$', '0', r'$\pi/2$', r'$\pi$'])
            # plt.ylabel('Phase [rad]')
            # plt.xlabel('Normalized frequency (1.0 = Nyquist)')
            # plt.show()
            # -------------------------------------------------------------------------------
        if self.tinnitus.rechtsRauschenLautstaerke:
            # (-3dB Grenzen) bzw was der Bandpass durchlässt
            fGrenz = [self.tinnitus.rechtsRauschenUntereGrenzfrequenz,
                      self.tinnitus.rechtsRauschenObereGrenzfrequenz]
            # sos (=second order sections = Filter 2. Ordnung) ist ein Array der Länge (filterOrder) und beinhaltet
            # die Koeffizienten der IIR Filter 2. Ordnung (b0, b1, b2 & a0, a1, a2)
            sos = signal.butter(5, fGrenz, 'bandpass', fs=self.framerate, output='sos')
            # sosfilt filtert das Signal mittels mehrerer 'second order sections' (= Filter 2. Ordnung) die über sos definiert sind
            outdata[:, 1] = signal.sosfilt(sos, outdata[:, 1])


        # Sinus addieren: f(t) = A * sin(2 * pi * f * t)
        for x in range(len(outdata)):
            # links: rauschen und sinus wieder zusammen addieren
            outdata[x][0] += self.tinnitus.linksLautstaerke * np.sin(2 * np.pi * self.tinnitus.linksFrequenz *
                                                                    ((x + self.start_idx) / self.framerate))
            # rechts: rauschen und sinus wieder zusammen addieren
            outdata[x][1] += self.tinnitus.rechtsLautstaerke * np.sin(2 * np.pi * self.tinnitus.rechtsFrequenz *
                                                                     ((x + self.start_idx) / self.framerate))

        self.start_idx += frames

        return outdata

    def musik_filtern(self):
        """
        Diese Funktion filtert die Tinnitus Frequenz aus einer gewählten Musikdatei. Dabei geht sie in 3 großen
        Schritten vor:
        1. Die nötigen Informationen über den Tinnitus aus der .csv Datei herausholen
        2. Die digitalen Filter erstellen und die Tinnitus Frequenz aus der Audiodatei "herausschneiden"
        3. Die fertigen Audiodatei als .wav Datei speichern
        """
        nframes = len(self.music_data)  # Gesamtanzahl der Frames in der Musikdatei

        # ------------1. Die nötigen Informationen über den Tinnitus aus der .csv Datei herausholen---------------------
        self.filterfortschritt = 1, 0  # der erste schritt
        csvstring = open("TinnitusDaten.csv").read()
        tinnitus_data = csvstring.split("\n")

        # linke Frequenz aus csv Datei holen
        lf = tinnitus_data[2]
        lf = lf.split(";")
        lf = float(lf[1])


        # rechte Frequenz aus csv Datei holen
        rf = tinnitus_data[7]
        rf = rf.split(";")
        rf = float(rf[1])


        # linke Lautstärke aus cvs Datei holen
        ll = tinnitus_data[3]
        ll = ll.split(";")
        ll = float(ll[1])

        # rechte Lautstärke aus cvs Datei holen
        rl = tinnitus_data[8]
        rl = rl.split(";")
        rl = float(rl[1])

        # -------- 2. Die digitalen Filter erstellen und die Tinnitus Frequenz aus der Audiodatei "herausschneiden------
        self.filterfortschritt = 2, 0  # der zweite schritt
        start_time = time.time()    # einen Timer laufen lassen um zu sehen wie lange Filterung dauert

        self.music_data = self.music_data/32767 # convert array from int16 to float
        """ OLD IIR Notch Filter 2nd Order----------------------------------------------------------------------
        w0 = float(lf / (self.music_samplerate / 2))  # Frequency to remove from a signal. If fs is specified, this is 
        in the same units as fs. By default, it is a normalized scalar that must satisfy 0 < w0 < 1, with w0 = 1 
        corresponding to half of the sampling frequency.
        Q = 30.0  # Quality factor. Dimensionless parameter that characterizes notch filter -3 dB bandwidth bw relative 
        to its center frequency, Q = w0/bw.
        b, a = signal.iirnotch(lf, Q, fs=self.music_samplerate)
        ----------------------------------------------------------------------------------------------------------------
        """

        """ New IIR Notch Filter 5th order--------------------------------------------------------------------------"""
        # ------------------------------------------LEFT EAR FILTERING-------------------------------------------------
        # Filterparameter festlegen------------
        order = 5   # Filterordnung
        bandwidth = 175     # Bandbreite des Sperrbereichs in Hz
        #stop_attenuation = 100  # minimum Attenuation (Damping, Reduction) in stop Band [only for elliptic filter necessary]
        cutoff_frequencies = [(lf - (bandwidth / 2)),(lf + (bandwidth / 2))]  # the cutoff frequencies (lower and upper)
        max_ripple_passband = 50    # Maximal erlaubte Welligkeit im Passbereich
        # -------------------------------------

        if ll != 0.0:   # nur wenn die Lautstärke des linken Tinnitus ungleich 0 ist, wird auf diesem Ohr auch gefiltert
            b, a = signal.iirfilter(order, cutoff_frequencies, rp=max_ripple_passband, btype='bandstop', ftype='butter',
                                    fs=self.music_samplerate)  # Diese Funktion erstellt den IIR-Bandpassfilter (links)

            music_links = signal.lfilter(b, a, self.music_data[:, 0])  # diese Funktion wendet den Filter an

            # FIR Filterversuch
            # b = signal.firwin(order, cutoff_frequencies, pass_zero="bandstop", fs=self.music_samplerate, width=bandwidth,
            #                   window="hamming")
            #
            # music_links = signal.lfilter(b, [1.0], self.music_data[:, 0])


        else:
            music_links = self.music_data[:, 0]  # ungefiltert, wenn kein Tinnitus angegeben wurde

        # ------------------------------------- RIGHT EAR FILTERING ------------------------------------------------
        if rl != 0.0:  # nur wenn die Lautstärke des rechten Tinnitus ungleich 0 ist, wird auf diesem Ohr auch gefiltert
            cutoff_frequencies = [(rf - (bandwidth / 2)), (
                        rf + (bandwidth / 2))]  # change the cutoff frequencies to the tinnitus of the RIGHT EAR

            b, a = signal.iirfilter(order, cutoff_frequencies, rp=max_ripple_passband, btype='bandstop', ftype='butter',
                                    fs=self.music_samplerate)  # Diese Funktion erstellt den IIR-Bandpassfilter (rechts)

            music_rechts = signal.lfilter(b, a, self.music_data[:, 1])    # rechts

            # FIR Filterversuch
            # b = signal.firwin(order, cutoff_frequencies, pass_zero="bandstop", fs=self.music_samplerate, width=bandwidth,
            #                   window="hamming")
            #
            # music_rechts = signal.lfilter(b, [1.0], self.music_data[:, 1])
        else:
            music_rechts = self.music_data[:, 1]  # diese Funktion filtert die Audiodaten(die Tinnitusfreq wird entfernt)

        endTimeFiltering = time.time()
        print("benötigte Zeit zum Filtern rechts Ohr =", endTimeFiltering - start_time, "s")

        #----------------------- 3. Maxima finden und Samples auf 1 normieren ----------------------------
        self.filterfortschritt = 3, 0  # der dritte schritt

        #Maximum finden (Funktion max(...) ist minimal schneller, macht aber Probleme beim Feedback)
        start_time = time.time()
        max_ges = 0
        fortschritt = 0
        for i in range(nframes):
            if max_ges < abs(music_links[i]):
                max_ges = abs(music_links[i])
            if max_ges < abs(music_rechts[i]):
                max_ges = abs(music_rechts[i])
            if i % int(nframes/10) == 0:    # glaub hier stand 10000 davor oder 50000
                # fortschritt = i / nframes * 100
                # self.filterfortschritt = 3, round(fortschritt, 1)
                # print("  max: ", self.filterfortschritt[1], "%")

                fortschritt += 10
                self.filterfortschritt = 3, round(fortschritt, 1)
                print("  max: ", self.filterfortschritt[1], "%")
        end_time = time.time()
        print("Zeitaufwand Maxima-Suche: ", end_time - start_time)

        #auf 4 Nachkommastellen aufrunden
        start_time = time.time()
        max_ges = int(max_ges * 10000)
        max_ges += 1
        max_ges /= 10000
        end_time = time.time()
        print("Zeitaufwand Maximum runden: ", end_time - start_time)

        #Alle samples normieren
        start_time = time.time()
        music_links /= max_ges
        music_rechts /= max_ges
        end_time = time.time()
        print("Zeitaufwand samples normieren: ", end_time - start_time)


        # ------------------------- 4. Die fertigen Audiodatei als .wav Datei speichern --------------------------------
        self.filterfortschritt = [4, 0]  # der vierte Schritt
        start_time = time.time()
        wav_obj = wave.open("MyTinnitusFreeSong.wav", "w")

        # Rahmenparameter für die .wav-Datei setzen
        wav_obj.setparams((self.nchannels, self.sampwidth, self.music_samplerate, nframes, self.comptype, self.compname))

        """The values are stored in a temporary list, and when the process is finished, they are joined together into 
        an string which is then sent to the output with the traditional writeframes method."""

        packedMusic = []    # Liste an die wir die einzelnen frames bereits in binär umgewandelt aneinanderreihen

        # Die Audiosamples schreiben
        print("Musikdatei wird erstellt...")
        fortschritt = 0
        for tinnitus_data in range(nframes):  #geht jeden Sample-Wert der Musikdatei einzeln durch
            # Die Audiodaten müssen von float in einen passenden int-Wert umgerechnet werden
            packedMusic.append(struct.pack('h', int(music_links[tinnitus_data] * 32767.0)))
            packedMusic.append(struct.pack('h', int(music_rechts[tinnitus_data] * 32767.0)))

            # wav_obj.writeframes(struct.pack('h', int(music_links[x] * 32767.0)))  # Werte für links und rechts werden bei
            # wav_obj.writeframes(struct.pack('h', int(musicRechts[x] * 32767.0)))  # wav abwechselnd eingetragen
            if tinnitus_data % int(nframes/10) == 0:
                fortschritt += 10
                self.filterfortschritt = 4, round(fortschritt, 1)
                print("  samples: ", self.filterfortschritt[1], "%")

        end_time = time.time()
        print("Zeitaufwand für das packen der einzelnen Samples =", end_time - start_time, "s")


        value_str = b"".join(packedMusic)

        start = time.time()
        wav_obj.writeframes(value_str)
        end = time.time()
        print("Zeitaufwand für das schreiben aller Frames in die wav Datei =", end - start, "s")
        wav_obj.close()

        print("Speichern beendet.")
        self.filterfortschritt = 5, 0  #Nach erfolgreichem Filtern Fortschritt zur Bestätigung auf 5 setzen
        # Plot (hilfreich für Filterentwurf)
        freq, h = signal.freqz(b, a, fs=self.music_samplerate)
        fig, ax = plt.subplots(2, 1, figsize=(8, 6))
        ax[0].plot(freq, 20 * np.log10(abs(h)), color='blue')
        ax[0].set_title("Frequency Response")
        ax[0].set_ylabel("Amplitude (dB)", color='blue')
        ax[0].set_xlim([0, 10000])
        ax[0].set_ylim([-120, 10])
        ax[0].grid()
        ax[1].plot(freq, np.unwrap(np.angle(h)) * 180 / np.pi, color='green')
        ax[1].set_ylabel("Angle (degrees)", color='green')
        ax[1].set_xlabel("Frequency (Hz)")
        ax[1].set_xlim([0, 10000])
        ax[1].set_yticks([-90, -60, -30, 0, 30, 60, 90])
        ax[1].set_ylim([-90, 90])
        ax[1].grid()
        plt.show()