Application Center - Maplesoft

App Preview:

Terzpegelspektrum

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

terzpegelspektrum.mws

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`;

Datei := `C:\\Messung.txt`

Es werden spter zwei Ausgabefiles verwendet. Hier werden deren Namen definiert.

> Ausgabe1:=cat(Datei ,`teil`);

Ausgabe1 := `C:\\Messung.txtteil`

> Ausgabe2:=cat(Datei ,`gesamt`);

Ausgabe2 := `C:\\Messung.txtgesamt`

> 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);

nt := 8192

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]:

> od:

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);

nterz := 42

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]:

> od:

Fr die spter folgende Analyse wird der Zeitverlauf in einzelne Segmente zerlegt. Hier wird die Lnge der Zeitintervalle festgelegt.

> tseg:=2;

tseg := 2

Die Zeitschrittweite wird aus dem Datenfile ausgelesen

> dt:=vlistv[2][1]-vlistv[1][1];

dt := .15625000e-2

Dauer der Aufzeichnung der Schwingung.

> T:=dt*nt;

T := 12.80000000

Fr die Fourier-Analyse transienter Schwingungen ist die Verwendung unterschiedlicher Fensterfunktionen blich. Hier werden die beiden meist benutzten Funktionen wahlweise verwendet.

> Rechteck:=t->1;

Rechteck := 1

> Hanning:=0.5*(1-cos(2*Pi*trel));

Hanning := .5-.5*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");

[Maple Plot]

> plot(Fenster(t),t=0..tseg, title="Zeitsegment");

[Maple Plot]

Die Mewerte liegen in der Einheit mm/s vor. Es wird im Programm aber mit m/s gerechnet. Die Werte werden deshalb hier umgerechnet.

> for i from 1 by 1 to nt do

> v[i]:=vlist[i]/1000:

> od:

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");

[Maple Plot]

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;

ntb := 13

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;

> od:

Die Anzahl der Werte innerhalb eines Intervalls ist

> nseg:=trunc(nt*tseg/T);

nseg := 1280

Die FFt-Analyse erfordert eine Anzahl von 2^n 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 2^n <= nT .

> nFFTs:=floor(evalf(log[2](nt)));

nFFTs := 13

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 nseg <= 2^n .

> nFFT:=ceil(evalf(log[2](nseg)));

nFFT := 11

Anzahl der bercksichtigten Werte im gesamten Zeitbereich

> nFs:=2**nFFTs;

nFs := 8192

Anzahl der bercksichtigten Werte je Intervall

> nF:=2**nFFT;

nF := 2048

Die Indizes der Intervallgrenzen werden festgelegt

> for i from 1 by 1 to ntb do

> iU[i]:=trunc((i-1)*nseg/2+1);

> iO[i]:=trunc(iU[i]+nF-1);

> od:

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:

> od:

Fr die ntb Zeitintervalle werden Felder angelegt, welche die jeweils nF Mewerte enthalten. Die Teilintervalle werden im folgenden grafisch in verschiedenen Farben dargestellt.

> for i from 1 by 1 to ntb do

> 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]):

> od:

> P1list:=[seq(P1[i],i=1..ntb)]:

> display(P1list);

[Maple Plot]

Ein Intervall wird exemplarisch herausgegriffen.

> display(P1[2]);

[Maple Plot]

Fr die konkreten Zeitschritte werden die Zahlenwerte der Fensterfunktion berechnet und als Feld abgelegt

Fr das Gesamtfeld:

> for i from 1 by 1 to nt do

> hs[i]:=evalf(Fensters(i*dt)):

> od:

Fr die Teilintervalle:

> for i from 1 by 1 to ntb do

> for j from 1 by 1 to nF do

> h[i][j]:=evalf(Fenster(j*dt)):

> od:

> od:

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:

> for i from 1 by 1 to ntb do

> for j from 1 by 1 to nF do

> xx[i][j]:=vband[i][j]*h[i][j]:

> yy[i][j]:=0:

> od:

> x[i]:=array([seq(xx[i][j],j=1..nF)]):

> y[i]:=array([seq(yy[i][j],j=1..nF)]):

> od:

Berechnung der Datenfelder fr das Gesamtzeitbereichs:

> for i from 1 by 1 to nFs do

> xs[i]:=v[i]*hs[i]:

> ys[i]:=0:

> od:

> 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;

dfs := .7812500000e-1

Die maximale Frequenz, die bercksichtigt wird, ist in Abhngigkeit von der Zeitschrittweite

> fmax:=1/dt/2;

fmax := 320.0000000

Fr die FFT -Analyse wird die vorgegeben Maple-Funktion verwendet.

> FFT(nFFTs,xs,ys):

Die Ergebnisse der FFT-Analyse mssen zunchst gewichtet werden.

> for i from 1 by 1 to nFs do

> 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]):

> od:

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);

rmax := .3862520942e-5

imax := .5623293852e-5

bmax := .5462429502e-5

pmax := 1.570682316

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:

> od:

> imaxr;imaxi;imaxb;imaxp;

175

125

175

3603

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;

94

256

> 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;

Fpeak := 13.59375000

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");

[Maple Plot]

> display([PFsi],title="Imaginrteil");

[Maple Plot]

> display([PFsb],title="Betrag");

[Maple Plot]

> display([PFsp],title="Phase");

[Maple Plot]

> display([PFsba,PFsbap],title="Bereich des Maximums");

[Maple Plot]

In entsprechender Weise wird die FFT-Analyse fr die Teilintervalle wiederholt

> for i from 1 by 1 to ntb do

> FFT(nFFT,x[i],y[i]):

> od:

Die Frequenzschrittweite ist hier

> df:=1/(nF*dt);

df := .3125000000

Fr die maximale Frequenz gilt wieder

> fmax:=1/dt/2;

fmax := 320.0000000

Die Resultate der FFT-Annalyse mssen wieder gewichtet werden

> for i from 1 by 1 to ntb do

> for j from 1 by 1 to nF do

> xh[i][j]:=x[i][j]/(nF/2):

> yh[i][j]:=y[i][j]/(nF/2):

> od:

> od:

Aus den oben definierten Terzen wird diejenige gesucht, die als letztes noch bercksichtigte Frequenzen beinhaltet

> terzmaxo:=0:

> for i from 1 by 1 to nterz do

> if (terzo[i] > fmax) then

> if (terzu[i] < fmax) then ioterz:=i fi:

> fi:

> od:

> ioterz;

24

Fr alle Teilintervalle werden aus dem Resultat der FFT-Analyse die Betrge der Spektren berechnet. Gleichzeitig wird der Absolute Maximalbetrag ermittelt.

> bmax:=0:

> for i from 1 by 1 to ntb do

> for j from 1 by 1 to nF do

> bh[i][j]:=sqrt(yh[i][j]**2+xh[i][j]**2):

> if (bh[i][j] > bmax) then bmax:=bh[i][j] fi:

> od:

> od:

> bmax;

.1374767087e-4

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):

> od:

> 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

> for i from 1 by 1 to ntb do

> 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)]:

> od:

> for i from 1 by 1 to ntb do

> 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):

> od:

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");

[Maple Plot]

> display(PFilist,title="Imaginrteil");

[Maple Plot]

In der Darstellung der Betrge werden die Terzgrenzen und -mittellinien eingetragen

> display(PFblist,Pterz,title="Betrag");

[Maple Plot]

Die Betrge werden als Animation dargestellt, um den Eindruck der vorbeifahrenden Straenbahn zu vermitteln

> display(PFb2list,insequence=true,view=0..1.1*bmax);

[Maple Plot]

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:

> od:

> od:

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);

[Maple Plot]

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);

v0 := 1/20000000

Aus den Spektralwerten wird der Pegel in dB bestimmt

Fr die Betrachtung der Teilintervalle

> for i from 1 by 1 to nF/2 do

> LvFmax[i]:=20*log[10](vfbmax[i]/v0):

> od:

Fr das Gesamtfeld

> for i from 1 by 1 to nFs/2 do

> LvFmaxs[i]:=20*log[10](bsh[i]/v0):

> od:

Bestimmung der Minimal- und Maximalwerte

Fr die Betrachtung der Teilintervalle

> Lmmax:=0:

> Lmmin:=0:

> for i from 1 by 1 to nF/2 do

> if (Lmmax < LvFmax[i]) then Lmmax:=LvFmax[i] fi:

> if (Lmmin > LvFmax[i]) then Lmmin:=LvFmax[i] fi:

> od:

Fr das Gesamtfeld

> Lmsmax:=0:

> Lmsmin:=0:

> for i from 1 by 1 to nFs/2 do

> if (Lmsmax < LvFmaxs[i]) then Lmsmax:=LvFmaxs[i] fi:

> if (Lmsmin > LvFmaxs[i]) then Lmsmin:=LvFmaxs[i] fi:

> od:

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):

> for i from 1 by 1 to ioterz do

> 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):

> od:

> 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):

> for i from 1 by 1 to ioterz do

> 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]]:

> od:

Darstellung der Spektralwerte in logarithmischem Mastab

Teilintervalle

> display(PFLv,PterzL,view=Lmmin*1.1..1.1*Lmmax);

[Maple Plot]

Gesamtfeld

> display(PFLsv,PterzL,view=Lmsmin*1.1..1.1*Lmsmax);

[Maple Plot]

Gesamtfeld und Teilintervalle

> display(PFLsv,PFLv,PterzL,view=min(Lmsmin,Lmmin)*1.1..1.1*max(Lmsmax,Lmmax));

[Maple Plot]

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

Teilintervalle

> for i from 1 by 1 to nF/2 do

> freq[i]:=(i-1)*df:

> od:

Gesamtfeld

> for i from 1 by 1 to nFs/2 do

> freqs[i]:=(i-1)*dfs:

> od:

Unterste Grenzfrequenz

> tstart:=terzu[1];

tstart := 1.4

Minimale Bandbreite

> dterzmin:=terzo[1]-terzu[1];

dterzmin := .4

Anzahl der Terzen, welche Frequenzwerte des betrachteten Bereiches enthalten, angefangen mit der oben definierten niedrigsten Terz

Teilintervalle

> kterz:=0:

> for i from 1 by 1 to nterz do

> if (terzu[i] < freq[nF/2]) then

> kterz:=kterz+1:

> fi:

> od:

> kterz;

24

Gesamtfeld (sollte den gleichen Wert wie kterz ergeben)

> ksterz:=0:

> for i from 1 by 1 to nterz do

> if (terzu[i] < freqs[nFs/2]) then

> ksterz:=ksterz+1:

> fi:

> od:

> ksterz;

24

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:

> for i from 1 by 1 to nF/2 do

> 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:

> od:

> 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:

> j:=1:

> for i from 1 by 1 to nFs/2 do

> 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:

> od:

> isterz[kterz]:=nFs/2:

Index der niedrigsten und der hchsten diskreten Frequenz innerhalb der Terzbnder

Fr die Teilintervalle:

> itu[1]:=0:

> for i from 1 by 1 to nF/2 do

> if (freq[i] > terzu[1]) then

> itu[1]:=i:

> i:=nF/2:

> fi:

> od:

> ito[1]:=iterz[1]:

> for i from 2 by 1 to kterz do

> itu[i]:=iterz[i-1]+1:

> ito[i]:=iterz[i]:

> od:

Fr den Gesamtzeitbereich:

> itsu[1]:=0:

> for i from 1 by 1 to nFs/2 do

> if (freqs[i] > terzu[1]) then

> itsu[1]:=i:

> i:=nFs/2:

> fi:

> od:

> itso[1]:=isterz[1]:

> for i from 2 by 1 to kterz do

> itsu[i]:=isterz[i-1]+1:

> itso[i]:=isterz[i]:

> od:

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:

> od:

> for j from itsu[i] by 1 to itso[i] do

> if (LvFmaxs[j] > Lvsterz[i]) then Lvsterz[i]:=LvFmaxs[j] fi:

> od:

> od:

Die Maximalwerte innerhalb der Terzbnder werden zunchst zustzlich in der identischen Darstellung wie oben eingetragen

Teilintervalle

> 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);

[Maple Plot]

Darstellung Gesamtzeitbereich

> display(PFLsv,PFLsvp,PterzL,Pvterzs,Pvterzps,Pterz);

[Maple Plot]

Schlielich werden die ermittelten Maximalwerte in dB innerhalb der Terzbnder in gleichmigen Abstnden entlang der x-Achse aufgetragen.

> Lmmin:=0:Lmmax:=0:

> for i from 1 by 1 to kterz do

> 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:

> od:

> for i from 1 by 1 to kterz do

> Pm||i:=[[i,Lmmin*1.1],[i,Lmmax*1.1]]:

> Pm[i]:=curve(Pm||i, color=cyan):

> od:

> 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):

> for i from 1 by 1 to kterz do

> TXTu[i]:=textplot([i,0.03*(Lmmax-Lmmin),terzm[i]],align={ABOVE}):

> od:

> for i from 1 by 1 to kterz do

> TXT[i]:=textplot([i,Lvterz[i]+0.03*(Lmmax-Lmmin),round(Lvterz[i])],align={ABOVE,RIGHT}):

> od:

> for i from 1 by 1 to kterz do

> TXTs[i]:=textplot([i,Lvsterz[i]-0.03*(Lmmax-Lmmin),round(Lvsterz[i])],align={BELOW,RIGHT}):

> od:

> TLu:=seq(TXTu[i],i=1..kterz):

> TL:=seq(TXT[i],i=1..kterz):

> TLs:=seq(TXTs[i],i=1..kterz):

> Tmin:=0:

> for i from 1 by 1 to kterz do

> if (terzm[i] = 4) then Tmin:=i fi:

> od:

> Tmax:=kterz:

Betrachtung der Teilintervalle

> display(Pvterzc,Pvterzpc,PmL,TL,TLu,view=[Tmin..Tmax,Lmmin*1.1..1.1*Lmmax]);

[Maple Plot]

Betrachtung des Gesamtzeitbereiches

> display(Pvterzcs,Pvterzpcs,PmL,TLs,TLu,view=[Tmin..Tmax,Lmmin*1.1..1.1*Lmmax]);

[Maple Plot]

Vergleich

> display(Pvterzc,Pvterzpc,Pvterzcs,Pvterzpcs,PmL,TLs,TL,TLu,view=[Tmin..Tmax,Lmmin*1.1..1.1*Lmmax]);

[Maple Plot]

> fd := fopen(Ausgabe1,WRITE,TEXT):

> writedata(fd,terzsp);

> fclose(fd);

> fd := fopen(Ausgabe2,WRITE,TEXT):

> writedata(fd,terzsps);

> fclose(fd);

>

>