Zurück zum Inhaltsverzeichnis

Kapitel 25

Fourierreihen und Fouriertransformation

Fourierreihen

> restart:

Die Koeffizienten

> ak:= 2/T0_ * Int( f_* cos(2*Pi*t_*k_/T0_), t_=0..T0_):

> bk:= 2/T0_ * Int( f_* sin(2*Pi*t_*k_/T0_), t_=0..T0_):

Int wird träge formuliert, weil sonst das Integral mit f_=const berechnet wird (eine Alternative wäre, f_(t_) zu verwenden, oder eine Prozedur zu schreiben...).

>

Bildung der Reihe

> formseries:=(f, t, T0, nmax) ->

> value(subs(k_=0, f_=f, t_=t, T0_=T0, ak/2)) +

> add( value(subs(k_=n, f_=f, t_=t, T0_=T0, ak)) * cos(2*Pi*n*t/T0), n=1..nmax)+

> add( value(subs(k_=n, f_=f, t_=t, T0_=T0, bk)) * sin(2*Pi*n*t/T0), n=1..nmax):

>

Ein Beispiel

> alias(sigma=Heaviside):

> form:=formseries(sigma(t)-sigma(t-1/2),t,2,7);

form := 1/4+cos(Pi*t)/Pi-1/3*cos(3*Pi*t)/Pi+1/5*cos...
form := 1/4+cos(Pi*t)/Pi-1/3*cos(3*Pi*t)/Pi+1/5*cos...

> simplify(%);

1/420*(105*Pi+420*cos(Pi*t)-140*cos(3*Pi*t)+84*cos(...
1/420*(105*Pi+420*cos(Pi*t)-140*cos(3*Pi*t)+84*cos(...

> plot(form, t=0..3);

[Maple Plot]

>

Ein zweites Beispiel

> form2:=formseries(t, t, 1, 15);

form2 := 1/2-sin(2*Pi*t)/Pi-1/2*sin(4*Pi*t)/Pi-1/3*...
form2 := 1/2-sin(2*Pi*t)/Pi-1/2*sin(4*Pi*t)/Pi-1/3*...
form2 := 1/2-sin(2*Pi*t)/Pi-1/2*sin(4*Pi*t)/Pi-1/3*...

> simplify(%);

-1/360360*(-180180*Pi+360360*sin(2*Pi*t)+180180*sin...
-1/360360*(-180180*Pi+360360*sin(2*Pi*t)+180180*sin...
-1/360360*(-180180*Pi+360360*sin(2*Pi*t)+180180*sin...
-1/360360*(-180180*Pi+360360*sin(2*Pi*t)+180180*sin...

>

> plot(form2, t=0..3);

[Maple Plot]

>

>

Verschiedene Varianten zur numerischen Berechnung der Koeffizienten.

1. Das Integral wird zuerst mit value symbolisch ausgewertet und dann numerisch dargestellt (kann weitere Rechnungen beschleunigen)

> evalfseries:=(f, t, T0, nmax) ->

> evalf(value(subs(k_=0, f_=f, t_=t, T0_=T0, ak/2))) +

> add( evalf(value(subs(k_=n, f_=f, t_=t, T0_=T0, ak))) * cos(2*Pi*n*t/T0), n=1..nmax)+

> add( evalf(value(subs(k_=n, f_=f, t_=t, T0_=T0, bk))) * sin(2*Pi*n*t/T0), n=1..nmax):

>

> evalfseries(sigma(t)-sigma(t-1/2),t,2,7);

.2500000000+.3183098861*cos(Pi*t)-.1061032954*cos(3...
.2500000000+.3183098861*cos(Pi*t)-.1061032954*cos(3...
.2500000000+.3183098861*cos(Pi*t)-.1061032954*cos(3...

>

2. Häufig ist eine symbolische Integration nicht möglich. Dann muß mit evalf(Int(...)) numerisch integriert werden:

> numseries:=(f, t, T0, nmax) ->

> evalf(subs(k_=0, f_=f, t_=t, T0_=T0, ak/2)) +

> add( evalf(subs(k_=n, f_=f, t_=t, T0_=T0, ak)) * cos(2*Pi*n*t/T0), n=1..nmax)+

> add( evalf(subs(k_=n, f_=f, t_=t, T0_=T0, bk)) * sin(2*Pi*n*t/T0), n=1..nmax):

>

> numseries(sigma(t)-sigma(t-1/2),t,2,7);

.2500000000+.3183098862*cos(Pi*t)+.3537785065e-16*c...
.2500000000+.3183098862*cos(Pi*t)+.3537785065e-16*c...
.2500000000+.3183098862*cos(Pi*t)+.3537785065e-16*c...
.2500000000+.3183098862*cos(Pi*t)+.3537785065e-16*c...
.2500000000+.3183098862*cos(Pi*t)+.3537785065e-16*c...

>

3. Numerische Integration mit Unterdrückung 'verschwindender Koeffizienten', die bei der numerischen Integration durch Rundungsfehler auftreten (s.o.)

> chop:=x->evalf(round( evalf(x)* 10^5 )/10^5, 5):

> chop(1/3), chop(1e-5), chop(1e-6);

.33333, .10000e-4, 0.

> chopseries:=(f, t, T0, nmax) ->

> chop(subs(k_=0, f_=f, t_=t, T0_=T0, ak/2)) +

> add( chop(subs(k_=n, f_=f, t_=t, T0_=T0, ak)) * cos(2*Pi*n*t/T0), n=1..nmax)+

> add( chop(subs(k_=n, f_=f, t_=t, T0_=T0, bk)) * sin(2*Pi*n*t/T0), n=1..nmax):

>

> chopseries(sigma(t)-sigma(t-1/2),t,2,7);

.25000+.31831*cos(Pi*t)-.10610*cos(3*Pi*t)+.63660e-...
.25000+.31831*cos(Pi*t)-.10610*cos(3*Pi*t)+.63660e-...

>

Komplexe Koeffizienten und FFT

> restart:

Das trigonometrische Interpolationspolynom mit komplexen Koeffizienten lautet

> Tn:=sum(c[k]*exp(2*Pi*I*k*t/T),k=0..n-1);

Tn := sum(c[k]*exp(2*I*Pi*k*t/T),k = 0 .. n-1)

>

Zwei Formulierungsmöglichkeiten für Fourierreihen mit reellen Koeffizienten sind:

> Tab:=a[0]/2+sum(a[k]*cos(2*Pi*k*t/T)+b[k]*sin(2*Pi*k*t/T),k=1..p-1)+a[p]/2*cos(2*Pi*p*t/T);

Tab := 1/2*a[0]+sum(a[k]*cos(2*Pi*k*t/T)+b[k]*sin(2...

>

mit p = floor((n+1)/2), a[0] = 2*c[0], a[p] = 2*c[p]

und a[i] = 2*Re(c[i]), b[i] = -2*Im(c[i])

>

> TA:=a[0]/2+sum(A[k]*sin(2*Pi*k*t/T+phi[k]),k=1..p);

TA := 1/2*a[0]+sum(A[k]*sin(2*Pi*k*t/T+phi[k]),k = ...

mit A[i] = sqrt(a[i]^2+b[i]^2), phi[i] = arctan(a[i],b[... und A[p] = a[p]/2, phi[p] = Pi/2

>

Zur Demonstration der Vorgehensweise werden die komplexen Koeffizienten mit der schnellen Fouriertransformation (FFT s.u.) berechnet. Dazu wird die zu analysierende Funktion 'digitalisiert':

Bereitstellen der benötigten Größen (tn ist der Zeitschritt, (tk|ftk) die Punkte)

> p:=floor((n+1)/2):

> tn:=k*T/n:

> tk:=tn $ k=0..n-1:

> ftk:='seq(f(tn),k=0..n-1)':

> data:=[tn,f(tn)] $ k=0..n-1:

> n:=2^m:

>

Für das Beispiel von oben wird zunächst eine stückweise definierte Funktion mit einer etwas modifizierten Heavisidefunktion aufgestellt (sonst 'undefined' für x=0).

> sprung:=proc(x) if x=0 then 1 else Heaviside(x) fi end;

sprung := proc (x) if x = 0 then 1 else Heaviside(x...

> f:=t->sprung(t)-sprung(t-1/2);

f := proc (t) options operator, arrow; sprung(t)-sp...

> plot(f,0..2);

[Maple Plot]

>

> T:=2: m:=5:

> #ftk;

>

Real- und Imaginärteil der Funktionswerte

> rx:=array([ftk]);

rx := vector([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0...

> ix:=array([0 $ j=1..n]);

ix := vector([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

>

Schnelle Transformation. Das Ergebnis steht wieder in den arrays rx, ix.

> evalhf(FFT(m,var(rx),var(ix))):

> #op(rx);op(ix);

Bilden der komplexen Koeffizienten c[i]. Mit etwas Nachhilfe bei Rundungsfehlern.

> for i to n do c[i-1]:= (rx[i]+I*ix[i])/n:

> if abs(Im(c[i-1]))<10^(-10) then c[i-1]:=Re(c[i-1]); fi ;

> if abs(Re(c[i-1]))<10^(-10) then c[i-1]:=I*Im(c[i-1]); fi

> od:

Die reellen Koeffizienten a[i], b[i] und A[i]

> a[0]:=2*c[0]: a[p]:=2*c[p]: phi[p]:=Pi/2: A[p]:=a[p]/2:

> for i to p-1 do
#a[i]:=c[i]+c[n-i];

> #b[i]:=I*(c[i]-c[n-i]);
a[i]:=2*Re(c[i]); b[i]:=-2*Im(c[i]); # vermeidet x+0.*I, neu in Maple 6 (?complex)

> A[i]:=sqrt(a[i]^2+b[i]^2);

> phi[i]:=arctan(a[i],b[i]);

> od:

> #'Tn'=Tn;'Tab'=Tab;'TA'=TA;

>

> seq(c[i],i=0..5);

.2500000000, .1742682873-.1430182873*I, .3125000000...
.2500000000, .1742682873-.1430182873*I, .3125000000...

> seq(a[i],i=0..5);

.5000000000, .3485365746, .6250000000e-1, -.7176744...

> seq(b[i],i=1..5);

.2860365746, .3142087182, .1342674440, -0., .272146...

Achtung: Wenn mit Floats gerechnet wird, bleit bei komplexen Zahlen 0.*I stehen:

> 2+0*I, 2.+0*I, 2+0.*I;

2, 2., 2.+0.*I

>

Vorgegebene Funktion mit "Meßpunkten" und Fourierreihe mit A[i]

> p2:=plot([data],style=point,color=red,symbol=box):

> plots[display]({plot({f(t),TA},t=-0.2*T..1.1*T),p2});

[Maple Plot]

Und zur Probe mit a[i], b[i]

> plots[display]({plot({f(t),Tab},t=-0.2*T..1.1*T),p2});

[Maple Plot]

>

Darstellung der komplexen Fourierreihe

> p1:=plot({evalc(Re(evalc(Tn))),evalc(Im(Tn))},t=0..1.1*T): #p1;

> plots[display]({plot({f(t),TA},t=0..1.1*T),p1,p2});

[Maple Plot]

>

Das Spektrum als Histogramm (Realteil, weil in A[p] 0.*I stehen kann)

> with(stats[statplots]):

> histogram([seq(Weight(i-1,Re(A[i])),i=1..p)],numbars=p,area=1);

[Maple Plot]

>

> histogram([seq(Weight(i-1,a[i]),i=1..p-1)],numbars=p-1,area=1);

> histogram([seq(Weight(i-1,b[i]),i=1..p-1)],numbars=p-1,area=1);

>

In den Phasen muss für den Plot erst Float(undefined) entfernt werden

> for i to p do if type(phi[i],Float(undefined)) then phi[i]:=0 fi od;

> histogram([seq(Weight(i-1,phi[i]),i=1..p-1)],numbars=p-1,area=1);

>

Lösung von Differentialgleichungen mit Fourierreihen

Beispiel: Ein Reihenschwingkreis wird mit einer Sägezahnspannung betrieben. Gesucht ist der Strom (oder die Spannung am Widerstand).

> restart:

Die Koeffizienten und die Reihenentwicklung werden wie oben definiert

> ak:= 2/T0_ * Int( f_* cos(2*Pi*t_*k_/T0_), t_=0..T0_):

> bk:= 2/T0_ * Int( f_* sin(2*Pi*t_*k_/T0_), t_=0..T0_):

> evalfseries:=(f, t, T0, nmax) ->

> evalf(value(subs(k_=0, f_=f, t_=t, T0_=T0, ak/2))) +

> add( evalf(value(subs(k_=n, f_=f, t_=t, T0_=T0, ak))) * cos(2*Pi*n*t/T0), n=1..nmax)+

> add( evalf(value(subs(k_=n, f_=f, t_=t, T0_=T0, bk))) * sin(2*Pi*n*t/T0), n=1..nmax):

>

Generatorspannung

> u:=(t-1/2*T)/T;

u := (t-1/2*T)/T

Berechnung der Koeffizienten

> bi:=seq( evalf(value(subs(k_=i, f_=u, t_=t, T0_=T, bk))), i=1..15);

bi := -.3183098861, -.1591549430, -.1061032954, -.7...
bi := -.3183098861, -.1591549430, -.1061032954, -.7...
bi := -.3183098861, -.1591549430, -.1061032954, -.7...

> ai:=seq( evalf(value(subs(k_=i, f_=u, t_=t, T0_=T, ak))), i=1..15);

ai := 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0...

(Stimmt - es handelt sich um eine ungerade Funktion.)

>

Reihenentwicklung

> u1:=evalfseries(u, t, T, 15);

u1 := -.3183098861*sin(2*Pi*t/T)-.1591549430*sin(4*...
u1 := -.3183098861*sin(2*Pi*t/T)-.1591549430*sin(4*...
u1 := -.3183098861*sin(2*Pi*t/T)-.1591549430*sin(4*...
u1 := -.3183098861*sin(2*Pi*t/T)-.1591549430*sin(4*...
u1 := -.3183098861*sin(2*Pi*t/T)-.1591549430*sin(4*...

>

Aufstellen der DGL für eine sinusförmige Generatorspannung (die Teillösungen werden später summiert). Der Ansatz mit u1 als Generatorspannung ist nicht empfehlenswert!

> eq := b*sin(omega*t)=R*i(t)+L*diff(i(t),t)+1/C*int(i(t),t);

eq := b*sin(omega*t) = R*i(t)+L*diff(i(t),t)+int(i(...

> de:= diff(eq, t);

de := b*cos(omega*t)*omega = R*diff(i(t),t)+L*diff(...

>

Verschiedene Lösungsmethoden (wenn Sie das #-Zeichen vor sol, lsol und fsol entfernen, können Sie die Lösungen ansehen):

1. Maple entscheidet selbst

> st:=time():

> sol:=dsolve( {de, i(0)=0, D(i)(0)=0}, i(t)):

> time()-st;

4.005

> #sol;

2. Laplacemethode

> st:=time():

> lsol:=dsolve( {de, i(0)=0, D(i)(0)=0}, i(t), method=laplace ):

> time()-st;

.902

> #lsol;

3. Fouriermethode

> st:=time():

> fsol:=dsolve( {de}, i(t), method=fourier ):

> time()-st;

.280

> #fsol;

Haben Sie den Unterschied in den Lösungen gefunden (die Rechenzeit ist nicht gemeint :-))?

Probe

> #subs(sol, de): simplify(%); # für Tüftler

> subs(lsol, de): simplify(%);

b*cos(omega*t)*omega = b*cos(omega*t)*omega

Umwandlung der Gleichungen (sol, lsol, fsol) in Terme

> isol:=subs( sol,i(t)):

> ilsol:=subs( lsol,i(t)):

> ifsol:=subs( fsol,i(t)):

> #ifsol;

Zur Orientierung (Kontrolle) für Zahlenwerte die Eigenfrequenz

> T:='T':R:='R':C:='C':L:='L':

> omr:=1/sqrt(L*C): Tr:=2*Pi/omr: om:=2*Pi/T:Lr:=1/(om^2*C):

Aufsummieren der Teillösungen

> isu:=add( subs( b=bi[i], omega=2*Pi/T*i, isol), i=1..15):

> ilsu:=add( subs( b=bi[i], omega=2*Pi/T*i, ilsol), i=1..15):

> ifsu:=add( subs( b=bi[i], omega=2*Pi/T*i, ifsol), i=1..15):

Konkrete Werte

> T:=1/100:

> plot(u1, t=0..4*T);

[Maple Plot]

>

> R:=10: C:=1e-4: L:=2e-2:

> T,evalf(Tr),evalf(Lr);

1/100, .8885765875e-2, .2533029590e-1

Darstellung der Generatorspannung und der Spannung am Widerstand, die mit der Laplacemethode und der Fouriermethode berechnet wurde

> colors:=color=[red,green,blue,black]:

> plot([u1,ilsu*R,evalc(Re(ifsu))*R], t=0..4*T,colors);

[Maple Plot]

Die 'Fourierlösung' schwingt im Gegensatz zur 'Laplacelösung' nicht ein, sondern ist periodisch. Das muß so sein, weil die Fouriertransformation (s.u.) einen rein imaginären Kern hat (deshalb auch keine Anfangsbedingungen).

>

Experimentieren Sie, unter welchen Bedingungen die Differenz der beiden Lösungen kann vernachlässigt werden kann (vielleicht hängt es mit der Dämpfung (R) zusammen?):

> plot(ilsu-evalc(Re(ifsu)), t=0..4*T);

[Maple Plot]

Falls sich jemand für die Lösung ohne spezielle Methode interessiert:

> #plot([u1,evalc(Re(isu))*R,ilsu*R,evalc(Re(ifsu))*R], t=0..4*T,colors);

Hier ist die Differenz zur 'Laplacelösung'

> plot(evalc(Re(isu))-ilsu, t=0..4*T);

[Maple Plot]

Also im Wesentlichen nur Rundungsfehler.

Noch zwei Experimentiervorschläge:

>

> T:=1/200:

> plot([u1,ilsu*R,evalc(Re(ifsu))*R], t=0..4*T,colors);

[Maple Plot]

>

> T:=1/20:R:=100:

> plot([u1,ilsu*R,evalc(Re(ifsu))*R], t=0..4*T,colors);

[Maple Plot]

>

Schnelle Fouriertransformation (FFT)

Es wird wieder ein Rechtecksimpuls simuliert

> restart:

> datar := array( [0$462, 1$100, 0$462]):

> datai := array( [0$1024] ):

> dataxy := [seq( [i, datar[i]], i=1..1000)]:

> plot(dataxy);

[Maple Plot]

Aufruf von FFT zur Berechnung von 2^10 Koeffizienten

> FFT(10, datar, datai);

1024

Und hier ist das Ergebnis der Transformation:

Realteil

> dataxy := [seq( [i, datar[i]],i=1..1024)]:

> plot(dataxy,0..1024);

[Maple Plot]

Imaginärteil (Plot bitte selbst erzeugen)

> dataxy := [seq( [i, datai[i]],i=1..1024)]:

> #plot(dataxy,0..1024);

>

Gibt es da eine Systematik?

> evalf(datar[1]+I*datai[1],5);

> for i from 2 to 5 do:

> evalf(datar[i]+I*datai[i],5), evalf(datar[1026-i]+I*datai[1026-i], 5);

> od;

100.

-98.438-.30201*I, -98.438+.30201*I

93.841+.57581*I, 93.841-.57581*I

-86.465-.79584*I, -86.465+.79584*I

76.720+.94154*I, 76.720-.94154*I

Ja! Der erste Koeffizient ist reell, der 'Rest' konjugiert-komplex bezüglich der Mitte des Arrays. Summiert man bei der Fourierreihe nicht von 0 bis n-1 (im Array von 1 bis n), sondern von -n/2 bis n/2, so läßt sich das auch so darstellen (Absolutbetrag der Koeffizienten):

>

> dataxy := [seq( [i-1, abs(datar[1024+i]+I*datai[1024+i]) ],i=-50..0), seq( [i-1, abs(datar[i]+I*datai[i]) ],i=1..50) ]:

> dataxy[45..55];

[[-7, 39.03425427], [-6, 52.35998277], [-5, 65.1482...
[[-7, 39.03425427], [-6, 52.35998277], [-5, 65.1482...

>

> plot(dataxy, -50..50);

[Maple Plot]

>

Rücktransformation

> iFFT(10, datar, datai);

1024

Genauigkeit?

> k:='k':datar[k] $ k=460..470;

.3450266415e-9, -.4562319860e-9, -.1444915806e-8, ....
.3450266415e-9, -.4562319860e-9, -.1444915806e-8, ....

>

Geschwindigkeitstests

> st:=time():

> FFT(10, datar, datai);

> iFFT(10, datar, datai);

> time()-st;

1024

1024

2.594

Mit Hardwarefloats

> st:=time():

> evalhf(FFT(10, var(datar), var(datai) ) );

> evalhf(iFFT(10, var(datar), var(datai) ) );

> time()-st;

>

1024.

1024.

.310

Genauikgeit?

> k:='k':datar[k] $ k=460..470;

.503248470624767865e-9, -.692478897998084088e-9, -....
.503248470624767865e-9, -.692478897998084088e-9, -....
.503248470624767865e-9, -.692478897998084088e-9, -....

(Läßt sich auch mehrmals ausführen, z.B. in einer Schleife...)

>

Glättung von Daten

Beispiel: Simulation eines Filters zur Rauschunterdrückung

> restart:

Zunächst eine reine Sinusschwingung

> datar:= [ seq( evalf(sin(50*2*Pi*t/1024)), t=1..1024 ) ]:

> dataxy:= [ seq( [i,datar[i]], i=1..1024) ]:

> plot(dataxy, 1..100);

[Maple Plot]

Und ihr Spektrum

> datar:=array(datar):

> datai:=array( [0$1024] ):

> ndatar:=copy(datar): ndatai:=copy(datai):

> FFT(10, ndatar, ndatai):

> dataxy:=[seq( [ i-1, abs(ndatar[i]+I*ndatai[i]) ],i=1..1024)]:

> plot(dataxy, 1..512);

[Maple Plot]

In Zahlen (nicht normiert):

> seq(abs(ndatar[i]+I*ndatai[i]),i=40..60);

0., 0., 0., .5692407186e-6, 0., 0., 0., .8735131003...
0., 0., 0., .5692407186e-6, 0., 0., 0., .8735131003...

>

Rauschen wird überlagert

> datar:= [ seq( evalf(datar[i]+rand(1000)()/100), i=1..1024) ]:

> dataxy:=[seq( [i,datar[i]], i=1..1024)]:

> plot(dataxy, 1..100);

[Maple Plot]

> datar:=array(datar):

> datai:=array( [0$1024] ):

> FFT(10, datar, datai):

> fdatar:=copy(datar):

> fdatai:=copy(datai):

> dataxy:=[seq( [ i, abs(datar[i]+I*datai[i]) ],i=2..1024)]:

(ohne c[0])

>

Nun sieht das Spektrum so aus

> plot(dataxy, 1..512);

[Maple Plot]

(die zweite Hälfte (der Koeffizienten) ist symmetrisch dazu)

>

Mit einem Filter werden Frequenzen oberhalb 100 abgeschnitten

> filter:=array([1$100, 0$924]):

> for i from 1 to 1024 do:

> datar[i] := datar[i]*filter[i]:

> datai[i] := datai[i]*filter[i]:

> od:

Und das Spektrum zurücktransformiert

> iFFT(10, datar, datai):

> dataxy:=[seq( [i,datar[i]], i=1..1024)]:

> plot(dataxy, 1..200);

[Maple Plot]

Filter von 40 bis 60

> filter:=array([0$40, 1$20, 0$964]):

> datar:=copy(fdatar):

> datai:=copy(fdatai):

> for i from 1 to 1024 do:

> datar[i] := datar[i]*filter[i]:

> datai[i] := datai[i]*filter[i]:

> od:

> iFFT(10, datar, datai):

> dataxy:=[seq( [i,datar[i]], i=1..1024)]:

> plot(dataxy, 1..200);

[Maple Plot]

Gaußfilter

> filter:=array( [0$40, seq( evalf(exp(-(i-10)^2/19)), i=1..19), 0$965] ):

> dataxy:=[seq( [i,filter[i]], i=1..1024)]:

> plot(dataxy, 30..70);

[Maple Plot]

>

> datar:=copy(fdatar):

> datai:=copy(fdatai):

> for i from 1 to 1024 do:

> datar[i] := datar[i]*filter[i]:

> datai[i] := datai[i]*filter[i]:

> od:

> iFFT(10, datar, datai):

>

> dataxy:=[seq( [i,datar[i]], i=1..1024)]:

> plot(dataxy, 1..200);

[Maple Plot]

>

Fouriertransformation

Geht man von der Fourierreihe zum Fourierintegral über, so kann man mit einem kontinuierlichen Spektrum auch aperiodische Funktionen der Zeit beschreiben.

> restart:

> with(inttrans):

Die Definition der Fouriertransformation lautet

> int(f(t)*exp(-I*w*t),t=-infinity..infinity);

int(f(t)*exp(-I*w*t),t = -infinity .. infinity)

Anmerkung: Die Ausführung des obigen Befehls zu gegebenem f(t) liefert 'undefined' wird also _nicht_ als Fouriertransformation erkannt.

>

Ein paar Beispiele, die auch Aufschluß über die Normierung geben.

Konstante Funktion (Frequenz 0)

> fourier(a,t,w);

2*a*Pi*Dirac(w)

> invfourier(%,w,t);

a

Impuls zur Zeit 0 erfordert alle Frequenzen mit konstantem Gewicht

> fourier(a*Dirac(t),t,w);

a

> invfourier(%,w,t);

a*Dirac(t)

>

Genau eine Frequenz

> invfourier(a*Dirac(w-w0),w,t);

1/2*a*exp(I*w0*t)/Pi

> fourier(%,t,w);

a*Dirac(-w+w0)

>

Fünf diskrete Frequenzen

> g:=add(Dirac(w-i)/i,i=1..5);

g := Dirac(w-1)+1/2*Dirac(w-2)+1/3*Dirac(w-3)+1/4*D...

> f:=invfourier(g,w,t);

f := 1/120*(60*exp(I*t)+30*exp(2*I*t)+20*exp(3*I*t)...

> plot(evalc(Re(f)),t);

[Maple Plot]

>

Das eigentliche Anwendungsgebiet der Fouriertransformation sind nichtperiodische Ereignisse

> alias(sigma=Heaviside):

> f:=sigma(t+1)-sigma(t-1);

f := sigma(t+1)-sigma(t-1)

> plot( f, t=-2..2, numpoints=500);

[Maple Plot]

> g:=fourier( f, t, w);

g := exp(I*w)*(Pi*Dirac(w)-I/w)-exp(-I*w)*(Pi*Dirac...

> g:=simplify(g);

g := 2*sin(w)/w

zu denen ein kontinuierliches Spektrum gehört

> plot({abs(g),abs(2/w)}, w=-4*Pi..4*Pi,0..2);

[Maple Plot]

>

> invfourier(g, w, t);

1/2*sigma(t+1)-1/2*sigma(-t-1)-1/2*sigma(t-1)+1/2*s...

> simplify(%);

sigma(t+1)-sigma(t-1)

>

Frequenzband

> g:=sigma(w-10)-sigma(w-20);

g := sigma(w-10)-sigma(w-20)

> f:=invfourier(g,w,t);

f := -1/2*(Pi*Dirac(t)*t+I)*(-exp(10*I*t)+exp(20*I*...

> f:=evalc(Re(f));

f := -1/2*(Pi*Dirac(t)*t*(-cos(10*t)+cos(20*t))+sin...

> combine(simplify(%),trig);

1/2*(-sin(10*t)+sin(20*t))/(Pi*t)

Es treten nur die Randfrequenzen auf!

> plot(2*Pi*evalc(Re(f)),t=-2..2);

[Maple Plot]

>

Gaußverteilung

> g:=exp(-(w-15)^2/2);

g := exp(-1/2*(w-15)^2)

> f:=invfourier(g,w,t);

f := 1/2*exp(-225/2)*sqrt(2)*exp(15*I*t)*exp(225/2-...

> f:=simplify(evalc(Re(f)));

f := 1/2*sqrt(2)*exp(-1/2*t^2)*cos(15*t)/(sqrt(Pi))...

> huelle:=remove(has,f,cos);

huelle := 1/2*sqrt(2)*exp(-1/2*t^2)/(sqrt(Pi))

>

> plot({f,huelle,-huelle},t=-3..3);

[Maple Plot]

>

>

Wie bei allen Integraltransformationen sind zwei Eigenschaften der F-Transformation besonders nützlich:

1. Transformation von Ableitungen

> restart:with(inttrans):

> alias(F=fourier):alias(F(omega)=fourier(f(t),t,omega)):

> seq(F(diff(f(t),t$n),t,omega),n=1..4);

I*omega*F(omega), -omega^2*F(omega), -I*omega^3*F(o...

>

Das überführt z.B. die Schwingungsgleichung

> dgl:=diff(f(t),t$2)+p*diff(f(t),t)+q*f(t)=sin(b*t);

> #dgl:=diff(f(t),t$2)+p*diff(f(t),t)+q*f(t)=Heaviside(t);

dgl := diff(f(t),`$`(t,2))+p*diff(f(t),t)+q*f(t) = ...

in eine lineare Gleichung für F(omega)

> bild:=F(dgl,t,omega);

bild := -omega^2*F(omega)+I*p*omega*F(omega)+q*F(om...

mit der Lösung

> bildsol:=solve(bild,F(omega));

bildsol := I*Pi*(-Dirac(-omega+b)+Dirac(omega+b))/(...

und der Rücktransformierten

> sol:=invfourier(bildsol,omega,t);

sol := -1/2*I*(exp(I*b*t)*b^2+I*exp(I*b*t)*p*b-exp(...

oder kürzer:

> dsolve(dgl,f(t),method=fourier);

f(t) = -1/2*I*(exp(I*b*t)*b^2+I*exp(I*b*t)*p*b-exp(...

>

2. Faltung

Bei jeder Messung werden die Daten gefiltert (Auflösungsvermögen und Ansprechwahrscheinlichkeit der Apparatur), so daß das Ergebnis immer eine Faltung ist:

> restart: with(inttrans):alias(F=fourier):

> faltung:=Int(filter(t-tau)*daten(tau),tau=-infinity..infinity);

faltung := Int(filter(t-tau)*daten(tau),tau = -infi...

Das Filter, von dem jeder Experimentalphysiker träumt, ist die Diracfunktion

> filter:=t->Dirac(t);faltung=value(faltung);filter:='filter':

filter := Dirac

Int(Dirac(-t+tau)*daten(tau),tau = -infinity .. inf...

die auch implizit so definiert ist: Sie reproduziert bei der Faltung die Originalfunktion exakt (ist also das 'Einselement' der Faltung).

>

Wir können nun etwas realistischere Filter untersuchen. Dazu brauchen wir zunächst den Faltungssatz

> alias(Ff(omega)=fourier(filter(t),t,omega)):alias(Fd(omega)=fourier(daten(t),t,omega)):

> faltung:=int(filter(t-tau)*daten(tau),tau=-infinity..infinity):

> bild:=F(faltung,t,omega);

bild := Fd(omega)*Ff(omega)

oder rückwärts

> expand(invfourier(bild,omega,t));

int(daten(-_U1)*filter(t+_U1),_U1 = -infinity .. in...

>

Angenommen es wurde eine Messung messung(t) gemacht und die Filterfunktion filter(t) ist bekannt. Dann kann man die Funktion daten(t) so erhalten

> bildgl:=F(messung(t),t,omega)=bild;

bildgl := F(messung(t),t,omega) = Fd(omega)*Ff(omeg...

> bildsol:=solve(bildgl,Fd(omega));

bildsol := F(messung(t),t,omega)/Ff(omega)

Das sind die transformierten Daten Fd(omega) und wir können die Lösung in allgemeiner Form schon paratstellen:

> sol:=invfourier(bildsol,omega,t);

sol := invfourier(F(messung(t),t,omega)/Ff(omega),o...

Nun können verschiedene 'Messungen' und Filtertypen eingesetzt werden. Der rote Input ist aktiv, der schwarze Input inaktiv. Wenn Sie den Cursor in eine Inputzeile stellen, können Sie mit dem Maplebutton (oder der rechten Maustaste) zwischen aktiv und inaktiv (executable und non executable) toggeln.

a) Rechteckfilter mit der halben Breite b und der Ansprechwahrscheinlichkeit a<1:

> a:='a':b:='b':t:='t':
messung:=t->5*cos(t)+sin(5*t)^2;
filter:=x->a/(2*b)*(Heaviside(x+b)-Heaviside(x-b));

messung := proc (t) options operator, arrow; 5*cos(...

filter := proc (x) options operator, arrow; 1/2*a*(...

b) Gaußfilter

b:='b':t:='t':daten:='daten':
messung:=t->5*cos(t)+sin(5*t)^2; filter:=x->exp(-x^2/b^2)/(b*sqrt(Pi));

messung := proc (t) options operator, arrow; 5*cos(...

filter := proc (x) options operator, arrow; exp(-x^...

c) Diracfilter

messung:=t->5*cos(t)+sin(5*t)^2; filter:=x->Dirac(x);

messung := proc (t) options operator, arrow; 5*cos(...

filter := Dirac

d) Glockenkurve (mit anderer Messung)

a:='a':b:='b':messung:=t->5*cos(5*t); filter:=x->a/(x^2+b)/normf;
normf:=int(1/(t^2+b),t=-infinity..infinity):

messung := proc (t) options operator, arrow; 5*cos(...

filter := proc (x) options operator, arrow; a/((x^2...

plot(filter(x),x=-2..2,0..2);

>

Werte für die Ansprechwahrscheinlichkeit und die Breite. Zeichnung der Messung, der Originaldaten, und der Filterfunktion:

> a:=1:b:=0.2:

> plot([evalc(Re(sol)),messung(t),filter(t)],t=-1..3,color=[red,blue,black]);

[Maple Plot]

>

Neue Filterbreite, Rechteckfilter, Gaußfilter , Diracfilter, Glockenkurve

Wenn Sie mit den Zahlen etwas experimentieren, werden Sie feststellen, dass man bei dem einfachen Modell des Rechteckfilters (und der hier angenommenen 'Messkurve') gewisse Filterbreiten vermeiden sollte.

Zur Übung können Sie versuchen, die Messung durch Faltung der Originaldaten mit der Filterfunktion zu reproduzieren. Im Falle des Rechteckfilters gelingt das z.B. so:

> a:='a':b:='b':daten:=unapply(evalc(Re(sol)),t):

> a:=0.8:b:=0.2:

> plot([evalc(Re(sol)),messung(t),faltung+0.1,filter(t)],t=-1..3,color=[red,blue,black]);

[Maple Plot]

Bei den anderen Filtern (außer Dirac) muß bei der Integration mit assume etwas nachgeholfen werden.

>

>

Zurück zum Inhaltsverzeichnis