PDF-Buch Moderne Physik mit Maple

Update auf Maple 10

Kapitel 3.1.2

Worksheet fft1_10.mws

>

c International Thomson Publishing Bonn        1995         filename: fft1.ms

Autor: Komma                                                                Datum:  10.8.94                       

Thema: Schnelle Fouriertransformation

Zur harmonischen Analyse diskreter Meßwerte kann man die schnelle Fourieranalyse verwenden. Der entsprechende Maple-Befehl lautet <FFT  und kann zur Beschleunigung der Berechnung zusammen mit <evalhf eingesetzt werden. Das trigonometrische Interpolationspolynom lautet (vgl. Bronstein):

>	restart;with(plots):#readlib(FFT):

>

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

>

Wobei die komplexen Koeffizienten von <FFT geliefert werden. Zwei weitere Formulierungsmöglichkeiten 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);

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

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:

>	#tk;ftk;data;

Stückweise definierte Funktionen mit Heaviside und nicht mit if (boolean-error).

>	T:=3: m:=6:

>

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

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

>

>

>	#g:=t->t*Heaviside(t)-2*(t-1)*Heaviside(t-1);

>	plot({f});

>

>

#ftk;

>

>

Real- und Imaginärteil der Funktionswerte

>	rx:=array([ftk]);

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

Ggf. "Meßfehler" einbauen (die dann aber nicht in data stehen)

>	#rx[7]:=9/20;rx[15]:=1/2;

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:

>	#c();a(); b();

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

read `fig.m`:winpl():opt2d;

vtitle:=FFT: vxlab:=Zeit: vylab:=Ampl:fontgr:=30:

Vorgegebene Funktion mit "Meßpunkten" und Fourierreihe  

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

#pspl(`ffft1.ps`):

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

>

#display({plot({f(t),TA},t=-0.2*T..1.1*T,opt2d)});

Nur zur Probe

>	display({plot({f(t),Tab},t=0..1.1*T),p2});

>

Darstellung der komplexen Fourierreihe (nur für m<4 sonst zu lange Rechenzeit bei starken Oszillationen)

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

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

>

Vergleich von Tab und TA

>	#plot(Tab-TA,t=0..T,y=0..10^(-9));

>

Darstellung des Spektrums

>	with(stats[statplots]):

>	#histogram([seq(Weight(i,A[i]),i=1..p)]);

>	# und hier muß noch dieses 0.I weg...

>	#seq(Weight(i,A[i]),i=1..p);

>

>	#seq(Weight(i,Re(A[i])),i=1..p);

>

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

>

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

>

vtitle:=Spektrum:vxlab:=Frequenz: vylab:=G:

pspl(`ffft2.ps`):

replot(histogram([seq(Weight(i,A[i]),i=1..p)]),opt2d);

Es gibt auch die Rücktransformation

>	evalhf(iFFT(m,var(rx),var(ix)));  #op(rx);op(ix);

>

komma@oe.uni-tuebingen.de