FourierTrans.mw

> restart;

> #F:=x->piecewise(x=0,1/2,x>0,x);
#F:=x->piecewise(x<1/2,x,x>=1/2,1-x);
#F:=x->piecewise(x<1/2,-1,x>1/2,1);
#F:=x->piecewise(x<1/2,-1,x>=1/2,1);
F:=x->piecewise(x>0 and x<1/2,-1,x>1/2,1);
plot(F(x),x=0..1);

F := proc (x) options operator, arrow; piecewise(0 < x and x < 1/2, -1, 1/2 < x, 1) end proc

[Plot]

> KK:=2;
N:=2^KK;L:=1-0;

KK := 2

N := 4

L := 1

Direct Integrations for Fourier Series

> for n from 0 to N do
a[n]:=2/L*int(F(x)*cos(2*Pi*n/L*x),x=0..L);
end do;

a[0] := 0

a[1] := 0

a[2] := 0

a[3] := 0

a[4] := 0

> for n from 0 to N do
b[n]:=2/L*int(F(x)*sin(2*Pi*n/L*x),x=0..L);
end do;

b[0] := 0

b[1] := -4/Pi

b[2] := 0

b[3] := -4/3/Pi

b[4] := 0

> for n from 0 to N do
 c[n]:=1/L*int(F(x)*exp(-I*2*Pi*n/L*x),x=0..L);
end do;
for n from 1 to N do
 c[-n]:=1/L*int(F(x)*exp(I*2*Pi*n/L*x),x=0..L);
end do;

Note that:
a[n]=c]n]+c[-n], b[n]=I (c[n]-c[-n])
c[-n]= 1/2 (a[n] +  b[n]), c[n]=1/2 (a[n] - I b[n])

c[0] := 0

c[1] := 2*I/Pi

c[2] := 0

c[3] := 2/3*I/Pi

c[4] := 0

c[-1] := -2*I/Pi

c[-2] := 0

c[-3] := -2/3*I/Pi

c[-4] := 0

> F1:=unapply(sum(evalf(c[i]*exp(I*2*Pi*i/L*x)),i=-(N-1)..(N-1)),x):

> plot({Re(F1(x)),F(x)},x=0..1);

[Plot]

Check zero of imaginary part due to the cancellation

> NN:=20;
for k from 0 to NN do
 F1(k/NN);
end do:

NN := 20

Orthogonality for the integration

> KK:=3;
N:=2^KK;L:=1-0;

KK := 3

N := 8

L := 1

> for k from 0 to N-1 do
c[k]:=evalf(sum(F(i*L/N)*exp(-I*2*Pi*k*i/N),i=0..N-1));
end do;

c[0] := 0.

c[1] := 4.828427124*I

c[2] := 0.

c[3] := .828427124*I

c[4] := 0.

c[5] := -.828427124*I

c[6] := 0.

c[7] := -4.828427124*I

> F1:=unapply(sum(evalf(c[i]*exp(I*2*Pi*i/L*x)/N),i=0..(N/2-1))+
sum(evalf(c[N-i]*exp(-I*2*Pi*i/L*x)/N),i=1..(N/2-1)),x):

> plot({Re(F1(x)),F(x)},x=0..1);

[Plot]

> NN:=20;
for k from 0 to NN do
 F1(k/NN);
end do:

NN := 20

>

FFT

> x1:=array([evalf(seq(F(i/N),i=0..N-1))]);
y1:=array([evalf(seq(0,i=0..N-1))]);

x1 := vector([0., -1., -1., -1., 0., 1., 1., 1.])

y1 := vector([0., 0., 0., 0., 0., 0., 0., 0.])

> FFT(KK,x1,y1);

8

> print(x1);print(y1);

vector([0., -0.1e-8, 0., 0.1e-8, 0., 0.1e-8, 0., -0.1e-8])

vector([0., 4.828427122, 0., .828427124, 0., -.828427124, 0., -4.828427122])

> F2:=unapply(sum(evalf((x1[i]+I*y1[i])*exp(I*2*Pi*(i-1)/L*x)/N),i=1..N/2)+
sum(evalf((x1[N-i+2]+I*y1[N-i+2])*exp(-I*2*Pi*(i-1)/L*x)/N),i=2..N/2),x):

> plot({Re(F2(x)),F(x)},x=0..1);

[Plot]

>