> #Aufgabe 30
> restart;
> with(linalg): with(plots):
Warning, new definition for norm
Warning, new definition for trace
> # Prozedur zur Berechnung des Fourierpolynoms S_n(f)
> # M Anzahl der verwendeten Funktionswerte zur Berechnung der Fourierkoeffizienten (Zweierpotenz)
> # Feld der Funktionswerte f(2 Pi k/M), k=0,...,M-1
> # n Grad des Fourierpolynoms
> # r Anzahl der Stellen an denen das Fourierpolynom ausgewertet wird.
> fourier:=proc(n,r,f)
> # option trace;
> global p;
> local dft,SandeTukey1,c,c1,j,M,q;
> SandeTukey1:=proc(a)
> # a Vektor von Zweierpotenzlaenge 
> # option trace;
> local N,m,m1,w,l,r,j,s,k,t,b,u,i,aa;
> N:=vectdim(a); b:=vector(N); aa:=vector(N); 
> for j from 1 to N do aa[j]:=a[j] od;
> t:=round(log[2](N)); m:=2*N;
> for i from t  by (-1) to 1 do
>   m:=m/2;  r:=round(N/m); w:=evalf(cos(2*Pi/m)-I*sin(2*Pi/m));
>   for j from 0 to r-1 do
>       for l from 0 to round(m/2)-1 do
>          s:=2*l*r; u:=round(s/2);
>          b[j+s+1]:=aa[j+u+1]+aa[j+u+round(N/2)+1];
>          b[j+s+r+1]:=(aa[j+u+1]-aa[j+u+round(N/2)+1])*w^l;
>       od
>   od;
>   for j from 1 to N do aa[j]:=b[j] od
> od;
> evalm(aa)
> end;
> p:=vector(r); q:=vector(r);
> M:=vectdim(f); c:=vector(M);
> c:=SandeTukey1(f);
> for j from 1 to M do c[j]:=c[j]/M od;
> c1:=vector(r,0);
> for j from 1 to n+1 do c1[j]:=c[j] od;
> for j from 0 to n-1 do c1[r-j]:=c[M-j] od;
> q:=SandeTukey1(c1);
> p[1]:=q[1];
> for j from 1 to r-1 do p[j+1]:=q[r+1-j] od;
> listplot([seq([2*Pi*j/r,Re(p[j+1])],j=0..r-1)]);
> end: 
> f:=proc(x) if (x>=0) and (x <=evalf(Pi/2)) then evalf(1-2*x/Pi) else 
> if (x>evalf(3*Pi/2)) and (x <= evalf(2*Pi)) then evalf(-3+2*x/Pi) else
> if (x>=evalf(Pi/2)) and (x <= evalf(3*Pi/2)) then 0 fi fi fi
> end:
> M:=16:
> f1:=vector(M):
> for j from 1 to M do f1[j]:=f(evalf(2*Pi*(j-1)/M)) od:
> r:=64:
> fourier(8,r,f1);
> F1:=listplot([seq([2*Pi*j/r,Re(p[j+1])],j=0..r-1)]):
> 

> M:=32:
> f1:=vector(M):
> for j from 1 to M do f1[j]:=f(evalf(2*Pi*(j-1)/M)) od:
> r:=64:
> fourier(8,r,f1);
> F2:=listplot([seq([2*Pi*j/r,Re(p[j+1])],j=0..r-1)]):

> with(plots): display({F1,F2});

>