Terzpegelspektrum
von Harald Kammerer maple@jademountain.de
> with(plots):with(plottools):
Warning, the name changecoords has been redefined
Die Farben werden zur Vereinfachung durchnummeriert
> cl[1]:=aquamarine:cl[2]:=black:cl[3]:=blue:cl[4]:=navy:cl[5]:=coral:cl[6]:=cyan:cl[7]:=brown:cl[8]:=gold:cl[9]:=green:cl[10]:=yellow:cl[11]:=grey:cl[12]:=khaki:cl[13]:=magenta:cl[14]:=maroon:cl[15]:=orange:cl[16]:=pink:cl[17]:=plum:cl[18]:=red:cl[19]:=sienna:cl[20]:=tan:cl[21]:=turquoise:cl[22]:=violet:cl[23]:=wheat:
> for i from 1 by 1 to 23 do
> cl[i+23]:=cl[i]:cl[i+46]:=cl[i]:cl[i+69]:=cl[i]:cl[i+92]:=cl[i]:cl[i+115]:=cl[i]:cl[i+138]:=cl[i]:cl[i+161]:=cl[i]:cl[i+184]:=cl[i]:cl[i+207]:=cl[i]:cl[i+230]:=cl[i]:cl[i+253]:=cl[i]:cl[i+276]:=cl[i]:cl[i+299]:=cl[i]:cl[i+322]:=cl[i]:cl[i+345]:=cl[i]:cl[i+368]:=cl[i]:cl[i+391]:=cl[i]:
> od:
Die zu analysierende Schwingung wird eingelesen. Verwendung findet hier ein gemessener Schwinggeschwindigkeitsverlauf. Im eingelesenen File sind die Daten einer gemessenen Schwingung an unterschiedlichen Messpunkten gespeichert. Es handelt sich um die berfahrt einer Straenbahn ber eine mit Stahlfederelementen schwingungsisolierte Weiche.
Im File sind folgende Werte enthalten:
Spalte 1: Zeit
Spalte 2: z-Komponente der gemessenen Schwingung neben der schwingungsisoliert gelagerten Gleistragplatte
Spalte 3: z-Komponente der gemessenen Schwingung auf der schwingungsisoliert gelagerten Gleistragplatte
> Datei:=`C:\\Messung.txt`;
Es werden spter zwei Ausgabefiles verwendet. Hier werden deren Namen definiert.
> Ausgabe1:=cat(Datei ,`teil`);
> Ausgabe2:=cat(Datei ,`gesamt`);
> vlistv:=readdata(Datei,3):
Im Beispiel wird die Schwingung neben der Platte betrachtet. Dieser Mepunkt liegt in einem Abstand von 2m neben der Weiche.
> wk:=2:
Die Anzahl der Zeitwerte wird aus der Zahl der Eintrge im Datenfeld bestimmt
> nt:=nops(vlistv);
Jetzt werden die Mewerte des betreffenden Mekanals aus dem Gesamtdatenfeld extrahiert.
> for i from 1 by 1 to nt do
> vlist[i]:=vlistv[i][wk]:
Die spter aus der Fourieranalyse erhaltenen Frequenzspektren werden in speziellen Bndern, den sogenannten Terzen dargestellt. Die Grenzen dieser Bnder sowie deren Mittelwerte sind der Literatur entnommen. Dieses Vorgehen hat seinen Ursprung in der Akustik. Im folgenden Feld sind die Grenzen und Mittelwerte der Terzen zusammengefat. Die niedrigste und die hchste bercksichtigte Terz sind hier willkrlich festgelegt. Whrend dem Pogrammablauf ergeben sich dijenigen Terzen, die tatschlich durch die Messung bercksichtigt werden.
> terz:=[[1.4, 1.8, 1.6], [1.8, 2.2, 2.0], [2.2, 2.8, 2.5], [2.8, 3.6, 3.2], [3.6, 4.5, 4.0], [4.5, 5.6, 5.0], [5.6, 7.1, 6.3], [7.1, 8.9, 8.0], [8.9, 11.2, 10.0], [11.2, 14.1, 12.5], [14.1, 17.8, 16.], [17.8, 22.4, 20.], [22.4, 28.2, 25.], [28.2, 35.5, 31.5], [35.5, 44.7, 40.], [44.7, 56.2, 50.], [56.2, 70.7, 63.], [70.7, 89.1, 80.], [89.1, 112., 100.], [112., 141., 125.], [141., 178., 160.], [178., 224., 200.], [224., 282., 250.], [282., 355., 315.], [355., 447., 400.], [447., 562., 500.], [562., 708., 630.], [708., 891., 800.], [891., 1122., 1000.], [1122., 1413., 1250.], [1413., 1778., 1600.], [1778., 2239., 2000.], [2239., 2818., 2500.], [2818., 3548., 3150.], [3548., 4467., 4000.], [4467., 5623., 5000.], [5623., 7079., 6300.], [7079., 8913., 8000.], [8913., 11220., 10000.], [11220., 14130., 12500.], [14130., 17780., 16000.], [17780., 22390., 20000.]]:
Anzahl der definierten Terzen
> nterz:=nops(terz);
Die Grenzen und Mittelwerte der Terzen werden in einzelne Felder geteilt
> for i from 1 by 1 to nterz do
> terzu[i]:=terz[i][1]:
> terzo[i]:=terz[i][2]:
> terzm[i]:=terz[i][3]:
Fr die spter folgende Analyse wird der Zeitverlauf in einzelne Segmente zerlegt. Hier wird die Lnge der Zeitintervalle festgelegt.
> tseg:=2;
Die Zeitschrittweite wird aus dem Datenfile ausgelesen
> dt:=vlistv[2][1]-vlistv[1][1];
Dauer der Aufzeichnung der Schwingung.
> T:=dt*nt;
Fr die Fourier-Analyse transienter Schwingungen ist die Verwendung unterschiedlicher Fensterfunktionen blich. Hier werden die beiden meist benutzten Funktionen wahlweise verwendet.
> Rechteck:=t->1;
> Hanning:=0.5*(1-cos(2*Pi*trel));
Es werden im folgenden sowohl der gesamte Zeitverlauf als auch die einzelnen, durch obige Segmentgre festgelegten Intervalle mittels Fourier-Analyse untersucht. Dementsprechend ist die Fensterfunktion sowohl fr den gesamten Zeitbereich als auch fr die Intervalle zu definieren. Das Rechteckfenster ist naturgem bereits entsprechend bestimmt, fr das Hanning-Fenster ist eine zustzliche Bestimmung notwendig.
Die Fensterfunktion fr die Betrachtung der Segmente wird mit "Fenster" bezeichnet, die fr die Betrachtung des Gesamtverlaufes mit "Fensters". Wahlweise ist die nichtverwendete Funktion im folgenden auszukommentieren.
Rechteckfenster
> #Fenster:=Rechteck:
> #Fensters:=Rechteck:
Hanningfenster
> Fenster:=t->subs(trel=t/tseg,Hanning):
> Fensters:=t->subs(trel=t/T,Hanning):
Grafische Darstellung der Fensterfunktion
> plot(Fensters(t),t=0..T, title="Gesamtzeitbereich");
> plot(Fenster(t),t=0..tseg, title="Zeitsegment");
Die Mewerte liegen in der Einheit mm/s vor. Es wird im Programm aber mit m/s gerechnet. Die Werte werden deshalb hier umgerechnet.
> v[i]:=vlist[i]/1000:
Der zu analysierende Geschwindigkeits-Zeit-Verlauf wird zunchst grafisch dargestellt.
> vT:=[seq([evalf((i-1)*dt),v[i]],i=1..nt)]:
> P0:=curve(vT, color=green):
> display([P0],title="Geschwindigkeit");
Die Zeitsegmente werden so festgelegt, da sie sich jeweils zur Hlfte berlappen. Damit ergibt sich fr die Anzahl der betrachteten Zeitsegmente
> ntb:=ceil(T/tseg)*2-1;
Fr die Zeitintervalle wird die Unter- und Obergrenze festgelegt
> for i from 1 by 1 to ntb do
> TU[i]:=(i-1)*(tseg/2);
> TO[i]:=TU[i]+tseg;
Die Anzahl der Werte innerhalb eines Intervalls ist
> nseg:=trunc(nt*tseg/T);
Die FFt-Analyse erfordert eine Anzahl von Werten mit einem ganzzahligen, positiven n. Wenn die vorliegende Anzahl an Werten dieser Regel nicht gengt, knnen entweder die restlichen Werte unbercksichtigt bleiben oder es werden die fehlenden Werte zu Null gesetzt. Die Entscheidung hngt von der konkreten Situation ab. Hier wird fr die Betrachtung des Gesamtzeitbereiches das erstgenannte Verfahren gewhlt. Zunchst wird deshalb die grte Zahl n bestimmt, fr die gilt .
> nFFTs:=floor(evalf(log[2](nt)));
Fr die Betrachtung der Zeitintervalle knnen auch Werte oberhalb der jeweiligen Intervallobergrenze benutzt werden, deshalb wird hier die kleinste Zahl n verwendet, fr die gilt .
> nFFT:=ceil(evalf(log[2](nseg)));
Anzahl der bercksichtigten Werte im gesamten Zeitbereich
> nFs:=2**nFFTs;
Anzahl der bercksichtigten Werte je Intervall
> nF:=2**nFFT;
Die Indizes der Intervallgrenzen werden festgelegt
> iU[i]:=trunc((i-1)*nseg/2+1);
> iO[i]:=trunc(iU[i]+nF-1);
Die Obergrenze des letzten Zeitintervalls ist hher als die Obergrenze des Gesamtbereichs. Deshalb wird der fehlende Bereich mit 0 gefllt.
> for i from nt+1 by 1 to iO[ntb] do
> v[i]:=0:
Fr die ntb Zeitintervalle werden Felder angelegt, welche die jeweils nF Mewerte enthalten. Die Teilintervalle werden im folgenden grafisch in verschiedenen Farben dargestellt.
> vTi[i]:=[seq([evalf((j-1)*dt),v[j]],j=iU[i]..iO[i])]:
> vband[i]:=[seq(v[j],j=iU[i]..iO[i])]:
> P1[i]:=curve(vTi[i], color=cl[i]):
> P1list:=[seq(P1[i],i=1..ntb)]:
> display(P1list);
Ein Intervall wird exemplarisch herausgegriffen.
> display(P1[2]);
Fr die konkreten Zeitschritte werden die Zahlenwerte der Fensterfunktion berechnet und als Feld abgelegt
Fr das Gesamtfeld:
> hs[i]:=evalf(Fensters(i*dt)):
Fr die Teilintervalle:
> for j from 1 by 1 to nF do
> h[i][j]:=evalf(Fenster(j*dt)):
Die Mewerte werden mit den diskreten Werten der Fensterfunktion multipliziert. Fr die Verwendung der FFT-Funktion von Maple sind die diskreten Mewerte in einem Feld abzulegen. Siehe hierzu die Maple-Hilfe zu FFT , wobei zu beachten ist, da die Mewerte hier real sind.
Berechnung der Datenfelder fr die Teilintervalle:
> xx[i][j]:=vband[i][j]*h[i][j]:
> yy[i][j]:=0:
> x[i]:=array([seq(xx[i][j],j=1..nF)]):
> y[i]:=array([seq(yy[i][j],j=1..nF)]):
Berechnung der Datenfelder fr das Gesamtzeitbereichs:
> for i from 1 by 1 to nFs do
> xs[i]:=v[i]*hs[i]:
> ys[i]:=0:
> xs:=array([seq(xs[i],i=1..nFs)]):
> ys:=array([seq(ys[i],i=1..nFs)]):
Auf Details der FFT-Analyse wird hier nicht nher eingegangen. Lediglich die Bestimmung einiger magebender Gren wird angegeben.
Die Ergebnisse der folgenden FFT-Analyse mssen mit Faktoren umgerechnet werden, da in der Analyse zunchst die Zeitschrittweite nicht bercksichtigt wird. Dazu wird die Frequenzschrittweite berechnet. Diese ergibt sich aus der Lnge des betrachteten Zeitintervalls. Fr das Gesamtdatenfeld gilt fr die Frequenzschrittweite
> dfs:=1/T;
Die maximale Frequenz, die bercksichtigt wird, ist in Abhngigkeit von der Zeitschrittweite
> fmax:=1/dt/2;
Fr die FFT -Analyse wird die vorgegeben Maple-Funktion verwendet.
> FFT(nFFTs,xs,ys):
Die Ergebnisse der FFT-Analyse mssen zunchst gewichtet werden.
> xsh[i]:=xs[i]/(nFs/2):
> ysh[i]:=ys[i]/(nFs/2):
> bsh[i]:=sqrt(xsh[i]**2+xsh[i]**2):
> psh[i]:=arctan(ysh[i]/xsh[i]):
Die Resultate der FFT-Analyse werden als Wertepaare mit der jeweils zugehrigen Frequenz als Feld abgespeichert.
Realteil:
> Fsr:=[seq([(i-1)*dfs,xsh[i]],i=1..nFs/2)]:
Imaginrteil:
> Fsi:=[seq([(i-1)*dfs,ysh[i]],i=1..nFs/2)]:
Betrag:
> Fsb:=[seq([(i-1)*dfs,bsh[i]],i=1..nFs/2)]:
Hier soll der Vollstndigkeit halber auch die Phase angegeben werden:
> Fsp:=[seq([(i-1)*dfs,psh[i]],i=1..nFs/2)]:
Hier werden betragsmigen Maxima der Spektren bestimmt
> WFsr:=seq(sqrt(xsh[i]**2),i=1..nFs/2):
> WFsi:=seq(sqrt(ysh[i]**2),i=1..nFs/2):
> WFsb:=seq(bsh[i],i=1..nFs/2):
> WFsp:=seq(sqrt(psh[i]**2),i=1..nFs/2):
> rmax:=max(WFsr);imax:=max(WFsi);bmax:=max(WFsb);pmax:=max(WFsp);
Bestimmung der Indezes der obigen Maxima
> imaxr:=0:
> imaxi:=0:
> imaxb:=0:
> imaxp:=0:
> for i from 1 by 1 to nFs/2 do
> if (WFsr[i] >= rmax) then imaxr:=i fi:
> if (WFsi[i] >= imax) then imaxi:=i fi:
> if (WFsb[i] >= bmax) then imaxb:=i fi:
> if (WFsp[i] >= pmax) then imaxp:=i fi:
> imaxr;imaxi;imaxb;imaxp;
Fr die sptere grafische Darstellung wird ein Bereich definiert, in dem das Resultat vergrert dargestellt wird.
> kU:=imaxb-trunc(nFs/100):
> kO:=imaxb+trunc(nFs/100):
> if (kU < 1) then kU:=1 fi;
> if (kO > nFs) then kO:=nFs fi;
> kU;kO;
> Fsba:=[seq([(i-1)*dfs,sqrt(ysh[i]**2+xsh[i]**2)],i=kU..kO)]:
Bestimmung der Frequenz mit dem extremalen Spektralwert
> Fpeak:=(imaxb-1)*dfs;
Grafische Darstellung
> PFsr:=curve(Fsr, color=green):
> PFsi:=curve(Fsi, color=red):
> PFsb:=curve(Fsb, color=blue):
> PFsp:=curve(Fsp, color=magenta, style=point):
> PFsba:=curve(Fsba, color=blue):PFsbap:=curve(Fsba, color=blue,style=point):
> display([PFsr],title="Realteil");
> display([PFsi],title="Imaginrteil");
> display([PFsb],title="Betrag");
> display([PFsp],title="Phase");
> display([PFsba,PFsbap],title="Bereich des Maximums");
In entsprechender Weise wird die FFT-Analyse fr die Teilintervalle wiederholt
> FFT(nFFT,x[i],y[i]):
Die Frequenzschrittweite ist hier
> df:=1/(nF*dt);
Fr die maximale Frequenz gilt wieder
Die Resultate der FFT-Annalyse mssen wieder gewichtet werden
> xh[i][j]:=x[i][j]/(nF/2):
> yh[i][j]:=y[i][j]/(nF/2):
Aus den oben definierten Terzen wird diejenige gesucht, die als letztes noch bercksichtigte Frequenzen beinhaltet
> terzmaxo:=0:
> if (terzo[i] > fmax) then
> if (terzu[i] < fmax) then ioterz:=i fi:
> fi:
> ioterz;
Fr alle Teilintervalle werden aus dem Resultat der FFT-Analyse die Betrge der Spektren berechnet. Gleichzeitig wird der Absolute Maximalbetrag ermittelt.
> bmax:=0:
> bh[i][j]:=sqrt(yh[i][j]**2+xh[i][j]**2):
> if (bh[i][j] > bmax) then bmax:=bh[i][j] fi:
> bmax;
Fr die grafische Darstellung werden die Grenz- (gelb) und Mittellinien (cyan) der Terzen festgelegt
> for i from 1 by 1 to ioterz do
> grnzoterz||i:=[[terzo[i],0],[terzo[i],bmax*1.1]]:
> mitteterz||i:=[[terzm[i],0],[terzm[i],bmax*1.1]]:
> grnzuterz||i:=[[terzu[i],0],[terzu[i],bmax*1.1]]:
> Pgreno[i]:=curve(grnzoterz||i, color=yellow):
> Pmitte[i]:=curve(mitteterz||i, color=cyan):
> Pgrenu[i]:=curve(grnzuterz||i, color=yellow):
> Pterz:=[seq(Pgreno[i],i=1..ioterz),seq(Pmitte[i],i=1..ioterz),seq(Pgrenu[i],i=1..ioterz)]:
Datenfeld fr die Grafik der Real - und Imaginrteile sowie die Betrge fr alle Zeitintervalle
> Fr||i:=[seq([(j-1)*df,xh[i][j]],j=1..nF/2)]:
> Fi||i:=[seq([(j-1)*df,yh[i][j]],j=1..nF/2)]:
> Fb||i:=[seq([(j-1)*df,bh[i][j]],j=1..nF/2)]:
> PFr[i]:=curve(Fr||i, color=cl[i]):
> PFi[i]:=curve(Fi||i, color=cl[i]):
> PFb[i]:=curve(Fb||i, color=cl[i]):
> PFb2[i]:=curve(Fb||i, color=blue):
Die Grafen der Spektren werden als Feld zusammengefat
> PFrlist:=[seq(PFr[i],i=1..ntb)]:
> PFilist:=[seq(PFi[i],i=1..ntb)]:
> PFblist:=[seq(PFb[i],i=1..ntb)]:
> PFb2list:=[seq(PFb2[i],i=1..ntb)]:
> display(PFrlist,title="Realteil");
> display(PFilist,title="Imaginrteil");
In der Darstellung der Betrge werden die Terzgrenzen und -mittellinien eingetragen
> display(PFblist,Pterz,title="Betrag");
Die Betrge werden als Animation dargestellt, um den Eindruck der vorbeifahrenden Straenbahn zu vermitteln
> display(PFb2list,insequence=true,view=0..1.1*bmax);
Jetzt wird eine Art "Peak-Hold" simuliert, indem aus der Folge der FFT-Analysen die Maximalwerte fr jede Frequenz bestimmt werden
> for i from 1 by 1 to nF/2 do
> vfbmax[i]:=0:
> for j from 1 by 1 to ntb do
> if (bh[j][i] > vfbmax[i]) then vfbmax[i]:=bh[j][i] fi:
Die Maximalwerte werden als Plot dargestellt (rot). Zum Vergleich ist das Spektrum des Gesamtsignals in blau mit eingezeichnet.
> Fbmax:=[seq([(i-1)*df,vfbmax[i]],i=1..nF/2)]:
> PFbmax:=curve(Fbmax, color=red):
> display({PFsb,PFbmax},view=0..1.1*bmax);
Die Darstellung des Frequenzspektrums in Terzen erfolgt in dB. Dazu mssen die Spektralwerte logaritmiert werden. Um dies zu ermglichen, mssen sie normiert werden. Als Bezug wird blicherweise (die Geschwindigkeit) gewhlt:
> v0:=5*10**(-8);
Aus den Spektralwerten wird der Pegel in dB bestimmt
Fr die Betrachtung der Teilintervalle
> LvFmax[i]:=20*log[10](vfbmax[i]/v0):
Fr das Gesamtfeld
> LvFmaxs[i]:=20*log[10](bsh[i]/v0):
Bestimmung der Minimal- und Maximalwerte
> Lmmax:=0:
> Lmmin:=0:
> if (Lmmax < LvFmax[i]) then Lmmax:=LvFmax[i] fi:
> if (Lmmin > LvFmax[i]) then Lmmin:=LvFmax[i] fi:
> Lmsmax:=0:
> Lmsmin:=0:
> if (Lmsmax < LvFmaxs[i]) then Lmsmax:=LvFmaxs[i] fi:
> if (Lmsmin > LvFmaxs[i]) then Lmsmin:=LvFmaxs[i] fi:
Im folgenden werden die grafischen Darstellungen vorbereitet. Es wird hier nicht jedes Detail beschrieben
Teilintervalle
> FLv:=[seq([(i-1)*df,LvFmax[i]],i=1..nF/2)]:
> PFLv:=curve(FLv, color=red):
> PFLvp:=curve(FLv, color=red,style=point):
> grnzoL||i:=[[terzo[i],Lmmin*1.1],[terzo[i],Lmmax*1.1]]:
> mitteL||i:=[[terzm[i],Lmmin*1.1],[terzm[i],Lmmax*1.1]]:
> grnzuL||i:=[[terzu[i],Lmmin*1.1],[terzu[i],Lmmax*1.1]]:
> PgrenoL[i]:=curve(grnzoL||i, color=yellow):
> PmitteL[i]:=curve(mitteL||i, color=cyan):
> PgrenuL[i]:=curve(grnzuL||i, color=yellow):
> PterzL:=[seq(PgrenoL[i],i=1..ioterz),seq(PmitteL[i],i=1..ioterz),seq(PgrenuL[i],i=1..ioterz)]:
Gesamtfeld
> FLsv:=[seq([(i-1)*dfs,LvFmaxs[i]],i=1..nFs/2)]:
> PFLsv:=curve(FLsv, color=blue):
> PFLsvp:=curve(FLsv, color=blue,style=point):
> grnzoLs||i:=[[terzo[i],Lmsmin*1.1],[terzo[i],Lmsmax*1.1]]:
> mitteLs||i:=[[terzm[i],Lmsmin*1.1],[terzm[i],Lmsmax*1.1]]:
> grnzuLs||i:=[[terzu[i],Lmsmin*1.1],[terzu[i],Lmsmax*1.1]]:
Darstellung der Spektralwerte in logarithmischem Mastab
> display(PFLv,PterzL,view=Lmmin*1.1..1.1*Lmmax);
> display(PFLsv,PterzL,view=Lmsmin*1.1..1.1*Lmsmax);
Gesamtfeld und Teilintervalle
> display(PFLsv,PFLv,PterzL,view=min(Lmsmin,Lmmin)*1.1..1.1*max(Lmsmax,Lmmax));
Die Darstellung der Spektren in Form von Terzen erfordert weitere Schritte: Die Frequenzachse wird zunchst gem den vordefinierten Terzen geteilt. Die Werte innerhalb der Terzbnder werden nach dem jeweiligen Maximalwert durchsucht und dieser wird schlielich in dB an der Stelle der Terzmittenfrequenz aufgetragen. Zuletzt werden die Terzmittenfrequenzen in gleichmigem Abstand entlang der x-Achse aufgetragen. Hier folgen die einzelnen Teilschritte diesen vorgehens ohne ausgiebige weitere Dokumentation.
Anlegen von Feldern fr die diskreten Frequenzwerte
> freq[i]:=(i-1)*df:
> freqs[i]:=(i-1)*dfs:
Unterste Grenzfrequenz
> tstart:=terzu[1];
Minimale Bandbreite
> dterzmin:=terzo[1]-terzu[1];
Anzahl der Terzen, welche Frequenzwerte des betrachteten Bereiches enthalten, angefangen mit der oben definierten niedrigsten Terz
> kterz:=0:
> if (terzu[i] < freq[nF/2]) then
> kterz:=kterz+1:
> kterz;
Gesamtfeld (sollte den gleichen Wert wie kterz ergeben)
> ksterz:=0:
> if (terzu[i] < freqs[nFs/2]) then
> ksterz:=ksterz+1:
> ksterz;
Zhlvariable fr die Betrachtung der Teilintervalle
iterz[i]: Index der hchsten diskreten Frequenz innerhalb des i-ten Terzbandes
> for i from 1 by 1 to kterz do iterz[i]:=0 od:
> j:=1:
> if (freq[i] <= terzo[j] and freq[i+1] > terzo[j] and i < nF/2) then iterz[j]:=i fi:
> if (freq[i+1] >= terzo[j] and i < nF/2) then j:=j+1 fi:
> iterz[kterz]:=nF/2:
Zhlvariable fr die Betrachtung des Gesamtzeitfeldes
isterz[i]: Index der hchsten diskreten Frequenz innerhalb des i-ten Terzbandes
> for i from 1 by 1 to kterz do isterz[i]:=0 od:
> if (freqs[i] <= terzo[j] and freqs[i+1] > terzo[j] and i < nFs/2) then isterz[j]:=i fi:
> if (freqs[i+1] >= terzo[j] and i < nFs/2) then j:=j+1 fi:
> isterz[kterz]:=nFs/2:
Index der niedrigsten und der hchsten diskreten Frequenz innerhalb der Terzbnder
> itu[1]:=0:
> if (freq[i] > terzu[1]) then
> itu[1]:=i:
> i:=nF/2:
> ito[1]:=iterz[1]:
> for i from 2 by 1 to kterz do
> itu[i]:=iterz[i-1]+1:
> ito[i]:=iterz[i]:
Fr den Gesamtzeitbereich:
> itsu[1]:=0:
> if (freqs[i] > terzu[1]) then
> itsu[1]:=i:
> i:=nFs/2:
> itso[1]:=isterz[1]:
> itsu[i]:=isterz[i-1]+1:
> itso[i]:=isterz[i]:
Die Maximalwerte innerhalb jedes Terzbandes werden aus den Spektren (in dB) gesucht
> for i from 1 by 1 to kterz do
> Lvterz[i]:=Lmmin:
> Lvsterz[i]:=Lmsmin:
> for j from itu[i] by 1 to ito[i] do
> if (LvFmax[j] > Lvterz[i]) then Lvterz[i]:=LvFmax[j] fi:
> for j from itsu[i] by 1 to itso[i] do
> if (LvFmaxs[j] > Lvsterz[i]) then Lvsterz[i]:=LvFmaxs[j] fi:
Die Maximalwerte innerhalb der Terzbnder werden zunchst zustzlich in der identischen Darstellung wie oben eingetragen
> terzsp:=[seq([terzm[i],Lvterz[i]],i=1..kterz)]:
> Pvterz:=curve(terzsp, color=magenta, thickness=3):
> Pvterzp:=curve(terzsp, color=magenta, style=point):
Gesamtzeitbereich
> terzsps:=[seq([terzm[i],Lvsterz[i]],i=1..kterz)]:
> Pvterzs:=curve(terzsps, color=brown, thickness=3):
> Pvterzps:=curve(terzsps, color=brown, style=point):
Darstellung Teilintervalle
> display(PFLv,PFLvp,PterzL,Pvterz,Pvterzp,Pterz);
Darstellung Gesamtzeitbereich
> display(PFLsv,PFLsvp,PterzL,Pvterzs,Pvterzps,Pterz);
Schlielich werden die ermittelten Maximalwerte in dB innerhalb der Terzbnder in gleichmigen Abstnden entlang der x-Achse aufgetragen.
> Lmmin:=0:Lmmax:=0:
> if (Lvsterz[i] > Lmmax) then Lmmax:=Lvsterz[i] fi:
> if (Lvsterz[i] < Lmmin) then Lmmin:=Lvsterz[i] fi:
> if (Lvterz[i] > Lmmax) then Lmmax:=Lvterz[i] fi:
> if (Lvterz[i] < Lmmin) then Lmmin:=Lvterz[i] fi:
> Pm||i:=[[i,Lmmin*1.1],[i,Lmmax*1.1]]:
> Pm[i]:=curve(Pm||i, color=cyan):
> PmL:=[seq(Pm[i],i=1..kterz)]:
> terzspc:=[seq([i,Lvterz[i]],i=1..kterz)]:
> terzspcs:=[seq([i,Lvsterz[i]],i=1..kterz)]:
> Pvterzc:=curve(terzspc, color=magenta, thickness=2):
> Pvterzcs:=curve(terzspcs, color=brown, thickness=2):
> Pvterzpc:=curve(terzspc, color=red, style=point, symbolsize=300):
> Pvterzpcs:=curve(terzspcs, color=blue, style=point, symbolsize=300):
> TXTu[i]:=textplot([i,0.03*(Lmmax-Lmmin),terzm[i]],align={ABOVE}):
> TXT[i]:=textplot([i,Lvterz[i]+0.03*(Lmmax-Lmmin),round(Lvterz[i])],align={ABOVE,RIGHT}):
> TXTs[i]:=textplot([i,Lvsterz[i]-0.03*(Lmmax-Lmmin),round(Lvsterz[i])],align={BELOW,RIGHT}):
> TLu:=seq(TXTu[i],i=1..kterz):
> TL:=seq(TXT[i],i=1..kterz):
> TLs:=seq(TXTs[i],i=1..kterz):
> Tmin:=0:
> if (terzm[i] = 4) then Tmin:=i fi:
> Tmax:=kterz:
Betrachtung der Teilintervalle
> display(Pvterzc,Pvterzpc,PmL,TL,TLu,view=[Tmin..Tmax,Lmmin*1.1..1.1*Lmmax]);
Betrachtung des Gesamtzeitbereiches
> display(Pvterzcs,Pvterzpcs,PmL,TLs,TLu,view=[Tmin..Tmax,Lmmin*1.1..1.1*Lmmax]);
Vergleich
> display(Pvterzc,Pvterzpc,Pvterzcs,Pvterzpcs,PmL,TLs,TL,TLu,view=[Tmin..Tmax,Lmmin*1.1..1.1*Lmmax]);
> fd := fopen(Ausgabe1,WRITE,TEXT):
> writedata(fd,terzsp);
> fclose(fd);
> fd := fopen(Ausgabe2,WRITE,TEXT):
> writedata(fd,terzsps);
>