Fourier Transform

FFT - Implementation 7.1

Here is the recursive Fast Fourier Transform from Implementation 7.1

FFT[invec_] := Module[{k,elist,olist,netlist,eft,oft,n,omega},
If[Length[invec] ==  1, 
Return[invec];
];
n = Length[invec];
omega = 1;
elist = Table[invec[[j]],{j,1,n,2}];
olist = Table[invec[[j]],{j,2,n,2}];
eft = FFT[elist];
oft = FFT[olist];
netlist = Table[0,{j,0,n-1}];
For[k = 0, k ≤ n/2 - 1, k = k+1,
netlist[[k+1]] = eft[[k+1]] + omega oft[[k+1]];
netlist[[k + n/2+1]] = eft[[k+1]] - omega oft[[k+1]];
omega = omega Exp[2 Pi I/n];
];
Return[netlist];
]

If we introduce a sign parameter, we can make the inverse FFT together with the FFT just by specifing sign = -1 or sign = 1:

FFT[invec_,sign_] := Module[{k,elist,olist,netlist,eft,oft,n,omega},
If[Length[invec] ==  1, 
Return[invec];
];
n = Length[invec];
omega = 1;
elist = Table[invec[[j]],{j,1,n,2}];
olist = Table[invec[[j]],{j,2,n,2}];
eft = FFT[elist,sign];
oft = FFT[olist,sign];
netlist = Table[0,{j,0,n-1}];
For[k = 0, k ≤ n/2 - 1, k = k+1,
netlist[[k+1]] = eft[[k+1]] + omega oft[[k+1]];
netlist[[k + n/2+1]] = eft[[k+1]] - omega oft[[k+1]];
omega = omega Exp[sign 2 Pi I/n];
];
Return[netlist];
]

The inverse FFT needs to be divided by “N”, as in (7.40).

iFFT[invec_] := FFT[invec,-1]/Length[invec]

FFT[invec_] := FFT[invec,1]

Speed - Figure 7.1

This is a plot of the ratio of the FFT timing to the DFT timing, i.e. ((N Log(N))/(N^2)) = ((n (2^n))/((((2^n)))^2)) = (n/(2^n)) - we plot on a continuous range, just to see the shape, only integer values of n are allowed in practice:

Plot[n/2^n,{n,1,10},PlotStyle->GrayLevel[0.5]]

One Dimensional Examples - Figures 7.2 and 7.3

Here we send in a well-defined f = 5 example, and expect to get delta spike approximations at +/- 5:

x[t_] := Cos[2 Pi 5 t]

Make sure we have a power of 2 as the data size and make the discrete version:

NN = 2^10;
T = 1.0;
dt = T/(NN);
df = 1/(NN dt);
X = N[Table[x[t],{t,0,1-dt,dt}]];

Now take the FFT, and rotate it so that the negative frequencies are in the right spot (see (7.39)), then generate the power spectrum by attaching the frequencies and taking the magnitude squared of the FFT coefficients:

FFTx = FFT[X];
FFTx = RotateRight[FFTx,NN/2];
PSx = Table[{jj df,Abs[ FFTx[[jj+ NN/2+1]]]^2},{jj,-NN/2,NN/2-1}];

Here’s the original function,

G0 = ListLinePlot[Table[{(j+1) dt, X[[j+1]]},{j,1,Length[X]-1}], PlotStyle->GrayLevel[0.5]]

and its Fourier transform

ListLinePlot[PSx,PlotRange->{{-10,10},All},PlotStyle->GrayLevel[0.5]]

That looks good - the peaks in the power spectrum appear at the correct fequency locations (and with appropriate height).

Next, we’ll try adding a decay envelope, to see what that does to the power spectrum:

x[t_] := Cos[2 Pi 5 t] Exp[-4 t]

NN = 2^10;
T = 1.0;
dt = T/(NN);
df = 1/(NN dt);
X = N[Table[x[t],{t,0,1-dt,dt}]];

FFTx = FFT[X];
FFTx = RotateRight[FFTx,NN/2];
PSx = Table[{jj df,Abs[ FFTx[[jj+ NN/2+1]]]^2},{jj,-NN/2,NN/2-1}];

Again, the original function,

G0 = ListLinePlot[Table[{t,N[x[t]]},{t,0,1,dt}], PlotRange->{-1,1},PlotStyle-> GrayLevel[0.5]]

and the power spectrum:

G1 = ListLinePlot[PSx, PlotRange->{{-15,15},All},PlotStyle->GrayLevel[0.5]]

The peaks broaden a bit, but remain at the relevant +/- 5 Hz locations.

Low Pass FIlter - Figures 7.4 and 7.5

Below, we take invec and cut all contributions to the power spectrum that come from cutabove and above - this is Implementation 7.2.  The total time of the sample, T, must be provided to correctly define the frequency step size: df = 1/T.

LoPass[invec_, cutabove_, T_] := Module[{outvec,NN,dt,df,cutindex,FTdata,index},
NN = Length[invec];
dt = N[T/NN];
df= N[1/(NN dt)];
FTdata = FFT[N[invec]];
cutindex = Round[cutabove/df];
outvec = Table[0.0,{j,1,Length[FTdata]}];
outvec[[1]] = FTdata[[1]];
For[index = 2, index ≤ cutindex, index=index+1,
outvec[[index]] = FTdata[[index]];
outvec[[Length[FTdata]-(index-2)]] = Conjugate[FTdata[[index]]];
];
outvec = iFFT[outvec];
Return[outvec];
]

The second line of the For loop ensures that the negative frequency data is related to the positive frequency data correctly (so that the inverse FFT is real-valued).

We’ll test it out on the “GraviteaTime.wav” sound file:

data = Import["http://people.reed.edu/~jfrankli/Book/downloads/files/GraviteaTime.wav","Data"][[1]];
mdata = Table[data[[j]],{j,1,2^16}];
dt = 1.0/(41000);
df = 1/(Length[mdata] dt);
jdata= Table[{j dt,mdata[[j]]},{j,1,Length[mdata]}];

Here’s the sound wave itself:

ListLinePlot[jdata, AxesOrigin-> {0,-1.1},PlotStyle->GrayLevel[0.5]]

And we can generate the power spectrum for this clip, to see what frequencies are relevant.

NN = Length[mdata];
FFTx = FFT[mdata];
FFTx = RotateRight[FFTx,NN/2];
PSx = Table[{jj df,Abs[ FFTx[[jj+ NN/2+1]]]^2},{jj,-NN/2,NN/2-1}];

ListLinePlot[PSx, PlotRange->{{0,2000},{0,10^7}},PlotStyle->GrayLevel[0.5]]

Using the built-in ListPlay command, we can play at a variety of SampleRate(s) - here we try 41 kHz:

ListPlay[mdata,SampleRate-> 41000]

Now use the lo-pass filter to get rid of frequencies above A, 440 Hz:

lojingle = LoPass[mdata, 440, Length[mdata] dt];

lodata  = Table[{j dt,lojingle[[j]]},{j,1,Length[lojingle]}];

ListLinePlot[lodata, AxesOrigin->{0,-1.1}, PlotStyle->GrayLevel[0.5]]

and finally, check the power spectrum to ensure that the frequencies above 440 Hz are, in fact, eliminated:

NN = Length[lojingle];
FFTx = FFT[lojingle];
FFTx = RotateRight[FFTx,NN/2];
PSx = Table[{jj df,Abs[ FFTx[[jj+ NN/2+1]]]^2},{jj,-NN/2,NN/2-1}];

ListLinePlot[PSx, PlotRange->{{0,2000},{0,1 10^7}},PlotStyle->GrayLevel[0.5]]

We can play the lo-passed result:

ListPlay[Chop[lojingle],SampleRate->41000]

The top pane of the ListPlay output is a spectrogram, to be discussed below.

Two Dimensional Fourier Transform - Figure 7.6

Here is a spatially oscillating function of two variables:

u[x_,y_] := Sin[2 Pi (4)x]  Sin[2 Pi (8) y]

Set up a grid with unity spatial size:

Nx= 32;
Ny = 32;
dx = 1.0/(Nx);
dy = 1.0/(Ny);

We can plot the function to see the two different wavelengths, one in each direction:

utab = Table[u[j dx, k dy],{j,0,Nx-1},{k,0,Ny-1}];
uxytab = Table[{j dx, k dy, u[j dx,k dy]},{j,0,Nx-1},{k,0,Ny-1}];
uxytab = Partition[Flatten[uxytab],3];
Ga = ListDensityPlot[uxytab,ColorFunction->"GrayTones"]

Implementation 7.3, the two-dimensional FFT - we apply the FFT to the rows first, then the columns.

FFT2D[uxy_] := Module[{FFTy,FFTret},
FFTy = Table[FFT[uxy[[k]]],{k,1,Length[uxy]}];
FFTy = Transpose[FFTy];
FFTret = Table[FFT[FFTy[[k]]],{k,1,Length[FFTy]}];
FFTret = Transpose[FFTret];
Return[FFTret];
]

pstab = FFT2D[utab];
pstab = Abs[pstab]^2;

pstab holds the power spectrum data, but now we need to rotate that data so that the positive and negative wavenumbers are in the correct locations:

pstab = Table[RotateRight[pstab[[j]],Ny/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
pstab =  Table[RotateRight[pstab[[j]],Nx/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];

We’ll take a spacing of 1 here, since 1/(Nx dx) = 1/(Ny dy) = 1.

dkx = 1.0;
dky = 1.0;

Attach the wavenumber information, and plot:

psdat = Table[{j dkx, k dky,pstab[[j+Nx/2+1,k + Ny/2+1]]},{j,-Nx/2,Nx/2-1},{k,-Ny/2,Ny/2-1}];
psdat = Partition[Flatten[psdat],3];
Gb = ListDensityPlot[psdat, PlotRange->All,ColorFunction->"GrayTones"]

Grout = GraphicsGrid[{{Ga,Gb}}]

Check the peak placement

Select[psdat,Abs[#[[3]] - Max[psdat]] < .000001&]

Filtering

We can filter images - below, we’ll read in the Roman Broccoli jpeg, convert it to grayscale, and perform both lo and hi-pass filters on its 2D Fourier transform:

data = Import["http://people.reed.edu/~jfrankli/Book/downloads/files/rb.jpg","Data"];

utab = N[Table[data[[j,k]].{255,255,255}/Sqrt[3 (255)^2],{j,65,65+511},{k,65,65+511}]];
Ga = MatrixPlot[utab]

We’ll assume a size of 1 for the side lengths, then we have dx and dy

Nx= Length[utab];
Ny = Length[utab[[1]]];
dx = 1.0/(Nx);
dy = 1.0/(Ny);

To recover a filtered image, we must define the inverse FFT for two dimensions:

iFFT2D[uxy_] := Module[{FFTy,FFTret},
FFTy = Table[iFFT[uxy[[k]]],{k,1,Length[uxy]}];
FFTy = Transpose[FFTy];
FFTret = Table[iFFT[FFTy[[k]]],{k,1,Length[FFTy]}];
FFTret = Transpose[FFTret];
Return[FFTret];
]

Below, we compute the two dimensional Fourier transform rotate it appropriately, and then cut any element with radius greater than or equal to 10, we then perform the inverse FFT to recover an image.  Since we are removing the elements with large wavenumber (the spatial analogue of frequency), and those correspond to smaller features, we expect the filtered image to have larger blobs in it, with the small details smoothed out.

pstab = N[FFT2D[utab]];

pstab = Table[RotateRight[pstab[[j]],Ny/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
pstab =  Table[RotateRight[pstab[[j]],Nx/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
cut = 10.;
For[i = 1, i ≤ Length[pstab], i = i+1,
For[j = 1, j≤ Length[pstab], j = j+1,
If[( (i-1)-Nx/2)^2 + ((j-1) - Ny/2)^2 ≥ cut^2,
pstab[[i,j]] = 0.0;
];
];
];
pstab = Table[RotateLeft[pstab[[j]],Ny/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
pstab =  Table[RotateLeft[pstab[[j]],Nx/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
idata=  iFFT2D[pstab];
Gl = MatrixPlot[idata]

This time, we cut those elements with radius less than 10, a high-pass filter - and again show the result - here, removing the low wavenumber components means that we are left with an image with a lot of small features.

pstab = N[FFT2D[utab]];

pstab = Table[RotateRight[pstab[[j]],Ny/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
pstab =  Table[RotateRight[pstab[[j]],Nx/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
cut = 10.;
For[i = 1, i ≤ Length[pstab], i = i+1,
For[j = 1, j≤ Length[pstab], j = j+1,
If[( (i-1)-Nx/2)^2 + ((j-1) - Ny/2)^2 ≤  cut^2,
pstab[[i,j]] = 0.0;
];
];
];
pstab = Table[RotateLeft[pstab[[j]],Ny/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
pstab =  Table[RotateLeft[pstab[[j]],Nx/2],{j,1,Length[pstab]}];
pstab = Transpose[pstab];
idata=  iFFT2D[pstab];
Gh = MatrixPlot[idata]

GraphicsGrid[{{Ga,Gl,Gh}}]

Spectrogram - Figure 7.7

An implementation of the spectrogram described in Section 7.7.4 - you provide the data, the low and high frequencies fl and fu.

SpectroGram[data_, fl_, fu_] := Module[{alpha ,ldata, index,NN,Ns, T,dt,df,FFTx,scale, oput},
NN = Length[data];
Ns = Round[2 fu/fl];
alpha = Round[NN/Ns];
T = alpha/fl;
dt = T/NN;
oput = {};
Print["Using : alpha = ", alpha, " Ns = ", Ns, " T = ", T];
For[index = 1, index ≤ alpha, index=index+1,
ldata = data[[(index-1) Ns + 1 ;; (index) Ns]];
FFTx = FFT[ldata];
FFTx =  RotateRight[FFTx,Ns/2];
FFTx = Table[{index Ns dt,jj (1/(Ns dt)),Abs[ FFTx[[jj+ Ns/2+1]]]^2},{jj,0,Ns/2-1}];
oput = Join[oput,{FFTx}];
];
oput = Partition[Flatten[oput],3];
scale = Max[oput];
oput = Table[{oput[[j,1]],oput[[j,2]],oput[[j,3]]/scale},{j,1,Length[oput]}];
Print[scale];
Return[oput];
]

Here’s a test function that goes linearly from f =4 to 8 Hz - the problem is set up in Section 7.7.4:

p[t_] := Cos[2 Pi (4) t] (T-t) + Cos[2 Pi (8) t] t

Take T = 16, with N = 512 and generate data.

T = 16;
NN = 512;
dt = T/NN;
ptab =N[ Table[p[t],{t,0,T-dt,dt}]];

We’d like fl = 1 Hz, fu = 16 Hz, so we can clearly see the transition from 4 to 8 in our data:

test = SpectroGram[ptab,1.,16.];

Here’s a three dimensional plot, so we can see what is happening

ListPlot3D[test, PlotRange->{{0,16},{0,16},{0,1}}]

and the more traditional contour plot:

Grout = ListContourPlot[test,PlotRange->{0,1},ColorFunction->"GrayTones"]