kill(all)$

maxima : discrete fourier transform 'gauss function',differentation,version 1.0
peter.vlasschaert@gmail.com,14/01/2021

ref : maple code : Denis C,Prof UGent. 'see website',see below code.

ratprint:false$

numerical : evaluation
-----------------------------

/*if numer#false then numer:false else numer:true;*/

given : function: (gauss function)
============

declare(N,integer)$

N:18 $ /* number of points : even number N<18 */
u:makelist(exp(-10 + 40*i*(1 - i/N)/N), i , 1, N)$
xy:makelist([i,u[i]],i,1,N)$
q1:['discrete,xy]$

example : N = 16, error floating-point-overflow.
use :
if numer#false then numer:false else numer:true;

wxplot2d([q1,q1],[style,points,lines],[color,red,blue]);

definition : "Discrete Fourier Transform". (2), inside exp
===============================

ss[uu,ii,jj,NN]:=uu[i+1]/NN*exp(2*%i*ii*jj*%pi/NN);

s1:ev(sum(ss[u,i,j,N],i,0,N-1),nouns)$

→ list : "Discrete Fourier Transform".

Fu:makelist(ev(s1,float,ratexpand),j,0,N-1)$

plot : discrete fourier transform (gauss function )
----------------------------------------------------------------

real values :

rFu:realpart(Fu)$

plot :

xyft:makelist([i,rFu[i]],i,1,N)$
qdft:['discrete,xyft]$
wxplot2d([qdft,qdft],[style,points,lines],[color,red,blue]);

rem :lot of deformation in the center ( = N/2 ),see below why.

approximation : differentation 'function'.
==============================

→ kFu 
----------------------------------------------------------
mathematics :   FT(df(x)/dx) = %i*k *FT(f(x))
----------------------------------------------------------

kFu:makelist(0,i1,1,N);

for k:1 - N/2 thru N/2 - 1 do(
   if 0 <= k then j : k else j :N + k, 
        kFu[j + 1] : ev(-2*%i*Fu[j + 1]*k*%pi/N,float,
        ratexpand)
                              )$

info:
ev(-2*%i*Fu[j + 1]*k*%pi/N,float,ratexpand) , correct=(-1) : k ,sequence.
       

→ correction : nyquist frequency

rr:N/2$

kFu[rr+1]:0;

→ final : kFu

ikFu:imagpart(kFu)$

rkFu:realpart(kFu)$

plot : discrete fourier transform (differentation of the function )
----------------------------------------------------------------------------------

xyf:makelist([i,rkFu[i]],i,1,N)$
qdf:['discrete,xyf]$
wxplot2d([qdf,qdf],[style,points,lines],[color,red,blue]);

plot :  not very good,between 4 & 14, jump .' see below'

definition : "(I)nverse Discrete Fourier Transform".  'differentation', (-2) exp
===================================================

si[uuu,i1,j1,NNN]:=uuu[i1+1]/NNN*exp(-2*%i*i1*j1*%pi/NNN);

si1:ev(sum(si[kFu,ii1,jj1,N],ii1,0,N-1),nouns)$
FikFu:realpart(makelist(ev(si1,float,ratexpand),jj1,0,N-1))$

plot : differentation of function
----------------------------------------

xyd:makelist([i,FikFu[i]],i,1,N)$

qd1:['discrete,xyd]$

wxplot2d([qd1,qd1],[style,points,lines],[color,red,blue]);

info:combination plot :exp(-x^2),diff(exp(-x^2),x)
-------------------------------------------------------------

wxplot2d([exp(-x^2),diff(exp(-x^2),x)],[x,-10,10]);

rem : use fft , ifft ,check definition can be different in other computer source.
------------------------------------------------------------------------------------------------------------------------------
see: website
appendix:
new: 04 /01/2021:maxima: relation between fft & inverse_fft ,version 1.0
download  (Maxima,code)       :maximabasicfftandinverseifftversion1
download  (Maxima,png)         :maximabasicfftandinverseifftversion1outp
new: 06 /01/2021:maxima: relation between fft & inverse_fft (1D,2D),version 1.0a
download  (Maxima,code)       :maximabasicfftandinverseifftversion1a
download  (Maxima,png)         :maximabasicfftandinverseifftversion1aoutp
new: 13 /01/2021:maxima : approximation 1e order differentation f(x) = f'(x) with fft,version 1.0
download  (Maxima,code)       :maxidifferentationapproxwithfftversion1final
download  (Maxima,png)         :maxidifferentationapproxwithfftversion1finaloutp
---------------------------------------------------------------------------------------------------------------------------------

maple code :
=========:
-----------------------------------------------------------------------------------
restart:
    N := 64:
    u := Array([seq(exp(-10 + 40*i*(1 - i/N)/N), i = 1 .. N)]):
    u0 := Array([seq(4*i/N, i = 1 .. N)]):
    print(plots[pointplot]([seq([i, u[i]], i = 1 .. N)])):
 Fu := Array([seq(
        evalf(add(u[i + 1]*exp(2*I*i*j*Pi/N)/N, i = 0 .. N - 1)),
        j = 0 .. N - 1)]);
    kFu := Array(1 .. N):
    for k from 1 - N/2 to N/2 - 1 do
        if 0 <= k then j := k else j := N + k end if:
        kFu[j + 1] := -2*I*Fu[j + 1]*k*Pi/N
    end do:
    kFu[N/2 + 1] := 0:
    FikFu := Array([seq(
        evalc(Re(add(kFu[i + 1]*exp(-2*I*i*j*Pi/N), i = 0 .. N - 1))),
        j = 0 .. N - 1)]):
    print(plots[pointplot]([seq([i, FikFu[i]], i = 1 .. N)])):

----------------------------------------------------------------------------------