Kapitel 25
Fourierreihen und Fouriertransformation
> 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);
> simplify(%);
> plot(form, t=0..3);
>
Ein zweites Beispiel
> form2:=formseries(t, t, 1, 15);
> simplify(%);
>
> plot(form2, t=0..3);
>
>
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);
>
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);
>
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);
> 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);
>
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);
>
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);
>
mit
und
>
> TA:=a[0]/2+sum(A[k]*sin(2*Pi*k*t/T+phi[k]),k=1..p);
mit und
>
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;
> f:=t->sprung(t)-sprung(t-1/2);
> plot(f,0..2);
>
> T:=2: m:=5:
> #ftk;
>
Real- und Imaginärteil der Funktionswerte
> rx:=array([ftk]);
> ix:=array([0 $ j=1..n]);
>
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);
> seq(a[i],i=0..5);
> seq(b[i],i=1..5);
Achtung: Wenn mit Floats gerechnet wird, bleit bei komplexen Zahlen 0.*I stehen:
> 2+0*I, 2.+0*I, 2+0.*I;
>
Vorgegebene Funktion mit "Meßpunkten" und Fourierreihe mit
> p2:=plot([data],style=point,color=red,symbol=box):
> plots[display]({plot({f(t),TA},t=-0.2*T..1.1*T),p2});
Und zur Probe mit
> plots[display]({plot({f(t),Tab},t=-0.2*T..1.1*T),p2});
>
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});
>
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);
>
> 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;
Berechnung der Koeffizienten
> bi:=seq( evalf(value(subs(k_=i, f_=u, t_=t, T0_=T, bk))), i=1..15);
> ai:=seq( evalf(value(subs(k_=i, f_=u, t_=t, T0_=T, ak))), i=1..15);
(Stimmt - es handelt sich um eine ungerade Funktion.)
>
Reihenentwicklung
> u1:=evalfseries(u, t, T, 15);
>
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);
> de:= diff(eq, t);
>
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;
> #sol;
2. Laplacemethode
> st:=time():
> lsol:=dsolve( {de, i(0)=0, D(i)(0)=0}, i(t), method=laplace ):
> time()-st;
> #lsol;
3. Fouriermethode
> st:=time():
> fsol:=dsolve( {de}, i(t), method=fourier ):
> time()-st;
> #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(%);
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);
>
> R:=10: C:=1e-4: L:=2e-2:
> T,evalf(Tr),evalf(Lr);
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);
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);
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);
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);
>
> T:=1/20:R:=100:
> plot([u1,ilsu*R,evalc(Re(ifsu))*R], t=0..4*T,colors);
>
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);
Aufruf von FFT zur Berechnung von 2^10 Koeffizienten
> FFT(10, datar, datai);
Und hier ist das Ergebnis der Transformation:
Realteil
> dataxy := [seq( [i, datar[i]],i=1..1024)]:
> plot(dataxy,0..1024);
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;
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];
>
> plot(dataxy, -50..50);
>
Rücktransformation
> iFFT(10, datar, datai);
Genauigkeit?
> k:='k':datar[k] $ k=460..470;
>
Geschwindigkeitstests
> st:=time():
> FFT(10, datar, datai);
> iFFT(10, datar, datai);
> time()-st;
Mit Hardwarefloats
> st:=time():
> evalhf(FFT(10, var(datar), var(datai) ) );
> evalhf(iFFT(10, var(datar), var(datai) ) );
> time()-st;
>
Genauikgeit?
> k:='k':datar[k] $ k=460..470;
(Läßt sich auch mehrmals ausführen, z.B. in einer Schleife...)
>
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);
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);
In Zahlen (nicht normiert):
> seq(abs(ndatar[i]+I*ndatai[i]),i=40..60);
>
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);
> 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);
(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);
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);
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);
>
> 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);
>
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);
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);
> invfourier(%,w,t);
Impuls zur Zeit 0 erfordert alle Frequenzen mit konstantem Gewicht
> fourier(a*Dirac(t),t,w);
> invfourier(%,w,t);
>
Genau eine Frequenz
> invfourier(a*Dirac(w-w0),w,t);
> fourier(%,t,w);
>
Fünf diskrete Frequenzen
> g:=add(Dirac(w-i)/i,i=1..5);
> f:=invfourier(g,w,t);
> plot(evalc(Re(f)),t);
>
Das eigentliche Anwendungsgebiet der Fouriertransformation sind nichtperiodische Ereignisse
> alias(sigma=Heaviside):
> f:=sigma(t+1)-sigma(t-1);
> plot( f, t=-2..2, numpoints=500);
> g:=fourier( f, t, w);
> g:=simplify(g);
zu denen ein kontinuierliches Spektrum gehört
> plot({abs(g),abs(2/w)}, w=-4*Pi..4*Pi,0..2);
>
> invfourier(g, w, t);
> simplify(%);
>
Frequenzband
> g:=sigma(w-10)-sigma(w-20);
> f:=invfourier(g,w,t);
> f:=evalc(Re(f));
> combine(simplify(%),trig);
Es treten nur die Randfrequenzen auf!
> plot(2*Pi*evalc(Re(f)),t=-2..2);
>
Gaußverteilung
> g:=exp(-(w-15)^2/2);
> f:=invfourier(g,w,t);
> f:=simplify(evalc(Re(f)));
> huelle:=remove(has,f,cos);
>
> plot({f,huelle,-huelle},t=-3..3);
>
>
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);
>
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);
in eine lineare Gleichung für F(omega)
> bild:=F(dgl,t,omega);
mit der Lösung
> bildsol:=solve(bild,F(omega));
und der Rücktransformierten
> sol:=invfourier(bildsol,omega,t);
oder kürzer:
> dsolve(dgl,f(t),method=fourier);
>
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);
Das Filter, von dem jeder Experimentalphysiker träumt, ist die Diracfunktion
> filter:=t->Dirac(t);faltung=value(faltung);filter:='filter':
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);
oder rückwärts
> expand(invfourier(bild,omega,t));
>
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;
> bildsol:=solve(bildgl,Fd(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);
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));
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:=t->5*cos(t)+sin(5*t)^2; filter:=x->Dirac(x);
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):
plot(filter(x),x=-2..2,0..2);
>
Werte für die Ansprechwahrscheinlichkeit und die Breite. Zeichnung der Messung, der Originaldaten, und der Filterfunktion:
> plot([evalc(Re(sol)),messung(t),filter(t)],t=-1..3,color=[red,blue,black]);
>
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]);
Bei den anderen Filtern (außer Dirac) muß bei der Integration mit assume etwas nachgeholfen werden.
>
>