Discrete Fourier Transform

Clear[x];
seed=3;
(* pick an even integer *)
Num=88;
GCD[seed,Num]

(* igrid is grid of Num values seperated by Δ=1 *)
(* note we start index at i=0, and proceed to N-1 *)
(* so matrix igrid[[1]] contains zer component etc. *)
igrid=Table[i,{i,0,Num-1}];
ordertable=Table[Mod[seed^igrid[[i]],Num],{i,1,Num}];
(* ordertable  contains the list of functions for the grid points in igrid *)
Dimensions[ordertable]
ListPlot[ordertable,Joined->True]

(* construct fourier space grid  (Subscript[ω, n]) = ((n Δ)/N)   where n=-N/2,-N/2+1, ..0, ... N/2 :   Note this grid contains N+1 values *)

fgrid=Table[n/Num,{n,-Num/2,Num/2}];

Dimensions[fgrid]

f[k_]:=1/Sqrt[Num] Sum[ordertable[[i]] Exp[2 Pi I (i-1) k /Num],{i,1,Num}]

fft=Chop[Fourier[ordertable]];
MatrixForm[fft];

mytable=Table[Chop[N[f[k]]],{k,-Num/2,Num/2}];
MatrixForm[mytable];

powerspectrum=Thread[{fgrid,Abs[mytable]^2}];

ListPlot[powerspectrum,Filling->Axis,PlotRange->{0,8900}]

ListPlot[ordertable,Joined->True]

ordertable

Mod[seed^10,Num]

Periodicity in the signal and the DFT

We noted above that if the function f(r) is periodic so that f(r)=f(r+a), then the DFT seems to show
enhanced peaks at the locations, in frequency space,  that are integer multiples of 1/a. In order to investigate
this further we again consider our expression for the DFT of  (Subscript[f, k])


 (Subscript[h, n]) =   (1/(Sqrt[N])) (∑_(k=0)^(N-1))  (Subscript[f, k]) (e^(2 π i n k/N))     
 
 Lets partition this sum in the following order 
 
 
 (1/(Sqrt[N]))(∑_(k=0)^(r-1))   (Subscript[f, k]) (e^(2 π i n k/N))  +  (1/(Sqrt[N]))(∑_(k=r)^(2r-1))   (Subscript[f, k]) (e^(2 π i n k/N))  +  (1/(Sqrt[N]))(∑_(k=2r)^(3r-1))   (Subscript[f, k]) (e^(2 π i n k/N)) + ...
 
 where r is the perid of f.  In the second sum we make a change of the sumation variable so that  k=m+r, also in the third sum
 we let k= m+ 2r and so on. Thus we get
 
 (1/(Sqrt[N]))(∑_(k=0)^(r-1))   (Subscript[f, k]) (e^(2 π i n k/N))  +  (1/(Sqrt[N]))(∑_(m=0)^(r-1))   (Subscript[f, m+r])    (e^(2 π i n r /N)) (e^(2 π i n m/N))  +  (1/(Sqrt[N]))(∑_(m=0)^(r-1))   (Subscript[f, m+2 r ])  (e^(2 π i n  (2r) /N)) (e^(2 π i n m /N)) +  ...

but since   (Subscript[f, k]) = (Subscript[f, k+r])   we can factor the following terms and simplify this expression as


 (1/(Sqrt[N]))(∑_(k=0)^(r-1))   (Subscript[f, k]) (e^(2 π i n k/N))  (1 +  (e^(2 π i n r /N)) + (e^(2 π i  (2n) r /N)) + ...  )
 
 The number of terms in the bracket is given by the smallest number of times r divides into N-1. Let’s assume that number is M  ⩽ N-1
 
  The terms in brackets  have the form of a geometric series
 
 1+ z + (z^2) + ... (z^(M-1))    where  z = (e^(2 π i n r /N))  we can rewrite that expression as  ( we neglect the remaining number (< r) of terms in the sum)
 
  (1/(Sqrt[N]))(∑_(k=0)^(r-1))   (Subscript[f, k]) (e^(2 π i n k/N)) ( (((e^(2 π i n r  M /N))-1)/( (e^(2 π i n r /N))-1)))  =   (1/(Sqrt[N]))(∑_(k=0)^(r-1))   (Subscript[f, k]) (e^(2 π i n k/N)) ( (((e^(2 π i n r M/N))  -1)/( (e^(2 π i n r /N))-1)))  
  
  where we have used the result  
  
  (∑_(k=0)^(M-1))(z^k) =  (((Z^M)-1)/(z-1))
  
  The term in brackets can be tacken outside of the summation sign and when calculating the power spectrum we
  obtain
  
|  ( (((e^(2 π i n r M/N))  -1)/( (e^(2 π i n r /N))-1))) ((|)^2)  | (1/(Sqrt[N]))(∑_(k=0)^(r-1))   (Subscript[f, k]) (e^(2 π i n k/N)) ((|)^2)


We can simplify the term in brackets by expanding the exponential factors in terms of trigonometric functions. We then get
for it 

((((Sin )^2)(πn r M/N))/(((Sin )^2)(π n r /N)))        which is always  less than or equal to unity, unless the factor in the denominator vanishes. This
occurs when   ((n r)/N)  is an integer.  In that case we can go to the original expression 

1+ z + (z^2) + ... (z^(M-1))  and insert z=1, in which case the sum adds to the value M, and the power spectrum

looks like

|  M |2   | (1/(Sqrt[N]))(∑_(k=0)^(r-1))   (Subscript[f, k]) (e^(2 π i n k/N)) ((|)^2)    for   ((n r)/N)  an integer

and

| (1/(Sqrt[N]))(∑_(k=0)^(r-1))   (Subscript[f, k]) (e^(2 π i n k/N)) ((|)^2)              for   ((n r)/N) not an integer


Therefore the power spectrum is amplified by a factor (M^2) if   ((n r)/N)=p   where p is an integer, or  the location of the 
peak  (for n ≠ 0)  occurs at   n= ((N p)/r) = (q/r)   where q is some integer.

(* take power spectrum of our fft transform and drop the first point coorepsonding to n=0 *)
fftpower=Drop[Thread[{igrid,Abs[fft]^2}],1];

ListPlot[fftpower,Filling->Axis,PlotRange->All]

There are local peaks at  n=35,26,18,9    

according to our formula r  is a multiple of 35/88=0.397,  26/88=0.295,  18/88 = 0.204,  and 9/88 = 0.102  (we know that r =10 so 1/r =0.1).
If we choose a Fourier array with larger length we should get more accurate results.

Bignum=280;
igrid=Table[i,{i,0,Bignum-1}];
ordertable=Table[Mod[seed^igrid[[i]],Num],{i,1,Bignum}];
fft=Chop[Fourier[ordertable]];
fftpower=Drop[Thread[{igrid,Abs[fft]^2}],1];

ListPlot[fftpower,Filling->Axis,PlotRange->All]

Table[fftpower[[i]],{i,15,30}]

Table[fftpower[[i]],{i,50,60}]

Table[fftpower[[i]],{i,80,90}]

So at n=28 we have a local maximum,  and therefore  28/280 = 0.1, also at n=56  56/280=0.2  and at n=84, 84/280=0.30