DSP Auswertung

Date: october 2015 Author: ESR

Notwendige Modulen



In [2]:

    
%reset -f
%matplotlib notebook
%load_ext autoreload
%autoreload 1
%aimport functions
import numpy as np
import copy
import acoustics
from functions import *
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
mpl.rcParams['lines.linewidth']=0.5
# uncomment next line to connect a qtconsole to the same session
#%qtconsole









    



The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Daten und Signale Importieren

Die passby Tabelle ist direct vom EMPA übernommen. Diese wird durch den script Ereigniss_to_json.py generiert. Folgende änderungen sind gemacht:

passby 14 wird als 'Regio' bezeichnet
passby 13 von 'Regio' to 'IC'



In [3]:

    
%%capture 
import json
passby = json.load(open('Tabellen\passby.json','r+'))
fill_passby_with_signals(passby)
#test
#passby = {k:passby[k] for k in['1','12','3']}

Verschiedene Grössen werden für jeden passby bzw. Abschnit berechnet

Für Details auswertungLS anschauen.

Für alle passby:

1. Anzahl Achsen 
2. Achsen Abstand

Für alle Abschnitte:

1. berechne bandpass signal

Für die Abschnitte Q1 und Q2

2. berechne tPeaks
3. berechne Speeds



In [4]:

    
%%capture
for pbID, pb in passby.items():
    #fill axledistance
    #https://de.wikipedia.org/wiki/Drehgestelltypen_%28Schweiz%29
    if pb['Zugstyp']=='IC':
        pb['axleDistance'] = 2.5
    else:
        pb['axleDistance'] = 2.5
        
    for A in ['Q1','Q2','Q3','Q4']:
        v = pb[A]
        #v['signals']['bandpass'] = v['signals']['MIC'].bandpass(20,20000)
        v['signals']['A'] = v['signals']['MIC'].weigh()
        if A in ['Q1','Q4']:
            v['tPeaks'] = detect_weel_times(v['signals']['LS'])
            iterator = zip(['vAv', 'dv', 'ti_vi'], train_speed(v['tPeaks'], axleDistance = pb['axleDistance']))
            v.update( {k:v for k,v in iterator} )
    assert(len(pb['Q1']['tPeaks'])==len(pb['Q4']['tPeaks']))
    pb['axleN'] = len(pb['Q1']['tPeaks'])

Spektren Berechnen

Für Details SpektrenBerechnen anschauen.

Auswertung Intervalle: Aufteilung des Zug definieren

für Abschnitten Q1 und Q4

full: gesamte signal
n: Wagenmitte der n-te Wagen bis zur Wagenmitte der n+1 Wagen

Die nächste funktion implementiert die Intervalle aus die Vorbeifahrtszeiten jedes Drehgestell



In [4]:

    
def defIntervals(tp):
    Intervals = {
                'full': (-np.inf,np.inf)
                #'vorbei': (tp.min(),tp.max()),
                }
    t = tp.reshape(len(tp)//4,4).mean(1)

    for n, (t1,t2) in enumerate(zip(t[:-1],t[1:])):
        Intervals[int(n+1)] = (t1,t2)
    return Intervals

für Abschnitten Q2 und Q3

full: gesamte signal



In [5]:

    
for pbID, pb in passby.items():
    # Intervals füllen
    for A in ['Q1','Q4']:
        v = pb[A]
        v['intervals'] = defIntervals(v['tPeaks'])
        
    for A in ['Q2','Q3']:
        v = pb[A]
        v['intervals'] = {'full': (-np.inf,np.inf)}

Passby mit die Terzpektren ausfüllen

1. passby mit intervalle füllen
2. passby mit `leq` Spektren für die definierte intervalle



In [6]:

    
Bands = acoustics.signal.OctaveBand(fstart = 100, fstop = 8000, fraction=3)#AB 100Hz wegen nahfeldeffekte
for pbID, pb in passby.items():
    for A in ['Q1','Q2','Q3','Q4']:
        # calc Octave
        v = pb[A]
        # wähle intervalle
        intervals=v['intervals']
        # berechne gefiterte signale
        sn = pb[A]['signals']['A']
        f , octFilterBank =  sn.third_octaves(frequencies = Bands.nominal)
        # leq
        spektrum, leq = cut_third_oct_spectrum( octFilterBank, intervals, lType= 'leq')
        v['spektrum'].update(spektrum)
        v.setdefault('leq',{}).update(leq)
        v['spektrum']['f'] = f.nominal

Werte in .json exportieren

die Classe MyEncoder ermöglicht das speichern von das dict passby wie folgt:

ohne Zeitsignalen
np.array werden als listen gespeichert



In [7]:

    
%%capture
passby['Beschreibung'] = "postprocessing der Vorbeifarten.\n Die postprocessing ist im `DSPAuswertung.ipynb` zu sehen. "
with open('results.json', 'w+') as data:
    json.dump( passby, data, cls = MyEncoder)