Discrete Fourier Transform (DFT) examples

Note: This notebook assumes that all definitions in ipcv_definitions.nb have been evaluated first. It also assumes that the directory in the Set image directory points to an image directory with the following images: cat256x256.ppm, part128x128.ppm and bakery256x256.ppm.

1D DFT: Example #1: single cosine

Sampled sequence: cos(2π t)

x[t_] := Cos[2  π t];
xd =  Sampler[x, fs = 10, {0, 4}];
If[OddQ[Length[xd]],  xd  =  Drop[xd, -1]];
gx  =  
        DiscreteListPlot[xd,  Blue,  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 500,  
          AspectRatio -> 1/3];

View as "image"

ShowImage[{Last  /@  xd}, Round[500/Length[xd]], 
        1.2, -1.2];

Magnitude plot of DFT as a function of index k

fft  =  DFT[xd];
g3 = 
        DiscreteListPlot[{#[[1]],  Abs[#[[2]]]} &  /@  fft, 
          Red,  PlotRange -> All,  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

View as "image"

tmp  =  Log[1 + Abs[Last  /@  fft]];
ShowImage[{tmp}, Round[500/Length[tmp]], 
        1.2  Max[tmp]];

Phase plot of DFT as a function of index k

fft  =  DFT[xd];
g4 = 
        DiscreteListPlot[{#[[1]],  
                Arg[#[[2]] + 0.0000001]} &  /@  fft, Red,  
          PlotRange -> {-π,  π},  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

Rearrange magnitude plot of DFT as a function of frequency f

g5 = 
        DiscreteListPlot[FourierRearrange[Abs[Last  /@  fft], fs], Red, 
          PlotRange -> All,  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

View as "image"

tmp  =  
        Log[1 + (Last  /@  
                FourierRearrange[Abs[Last  /@  fft], 
                  fs])];
ShowImage[{tmp}, Round[500/Length[tmp]], Max[tmp]];

1D DFT: Example #2: DC + sum of cosines

Sampled sequence: 2 + cos(2π t) + 2 cos(4π t)

x[t_] := 2 + Cos[2  π t] + 2 Cos[4  π t];
xd =  Sampler[x, fs = 10, {0, 4}];
If[OddQ[Length[xd]],  xd  =  Drop[xd, -1]];
gx  =  
        DiscreteListPlot[xd,  Blue,  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 500,  
          AspectRatio -> 1/3];

View as "image"

ShowImage[{Last  /@  xd}, Round[500/Length[xd]], 6, 
        0];

Magnitude plot of DFT as a function of index k

fft  =  DFT[xd];
g3 = 
        DiscreteListPlot[{#[[1]],  Abs[#[[2]]]} &  /@  fft, 
          Red,  PlotRange -> All,  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

View as "image"

tmp  =  Log[1 + Abs[Last  /@  fft]];
ShowImage[{tmp}, Round[500/Length[tmp]], 
        1.2  Max[tmp]];

Phase plot of DFT as a function of index k

fft  =  DFT[xd];
g4 = 
        DiscreteListPlot[{#[[1]],  
                Arg[#[[2]] + 0.0000001]} &  /@  fft, Red,  
          PlotRange -> {-π,  π},  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

Rearrange magnitude plot of DFT as a function of frequency f

g5 = 
        DiscreteListPlot[FourierRearrange[Abs[Last  /@  fft], fs], Red, 
          PlotRange -> All,  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

View as "image"

tmp  =  
        Log[1 + (Last  /@  
                FourierRearrange[Abs[Last  /@  fft], 
                  fs])];
ShowImage[{tmp}, Round[500/Length[tmp]], Max[tmp]];

1D DFT: Example #3: step (edge)

Sampled sequence

fs = 10;
xd =  
        Table[-1, {20}]~Join~Table[1, {20}];
If[OddQ[Length[xd]],  xd  =  Drop[xd, -1]];
gx  =  
        DiscreteListPlot2[xd,  Blue,  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 500,  
          AspectRatio -> 1/3];

View as "image"

ShowImage[{xd}, Round[500/Length[xd]], 2, -2];

Magnitude plot of DFT as a function of index k

fft  =  DFT[xd];
g3 = 
        DiscreteListPlot[{#[[1]],  Abs[#[[2]]]} &  /@  fft, 
          Red,  PlotRange -> All,  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

View as "image"

tmp  =  Log[1 + Abs[Last  /@  fft]];
ShowImage[{tmp}, Round[500/Length[tmp]], 
        1.2  Max[tmp]];

Phase plot of DFT as a function of index k

fft  =  DFT[xd];
g4 = 
        DiscreteListPlot[{#[[1]],  
                Arg[#[[2]] + 0.0000001]} &  /@  fft, Red,  
          PlotRange -> {-π,  π},  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

View as "image"

tmp  =  
        Arg[(Last  /@  
                fft) +  .0000001 ];
ShowImage[{tmp}, Round[500/Length[tmp]], 
        1.2  π, -1.2 π];

Rearrange magnitude plot of DFT as a function of frequency f

g5 = 
        DiscreteListPlot[FourierRearrange[Abs[Last  /@  fft], fs], Red, 
          PlotRange -> All,  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

View as "image"

tmp  =  
        Log[1 + (Last  /@  
                FourierRearrange[Abs[Last  /@  fft], 
                  fs])];
ShowImage[{tmp}, Round[500/Length[tmp]], Max[tmp]];

Rearrange phase plot of DFT as a function of frequency f

g5 = 
        DiscreteListPlot[
          FourierRearrange[Arg[Last  /@  fft  +   .0000001], fs], Red, 
          PlotRange -> All,  Frame -> True,  
          Axes -> False,  AspectRatio -> 1/3,  
          ImageSize -> 500];

View as "image"

tmp  =  
        Last /@  FourierRearrange[Arg[Last  /@  fft  +   .0000001], 
            fs];
ShowImage[{tmp}, Round[500/Length[tmp]], 
        1.2  π, -1.2 π];

Filtering in the frequency domain

lpf = Filter[Length[xd], fs, {0, 2}];

DiscreteListPlot[FourierRearrange[lpf,  fs], Magenta, 
        PlotRange -> All,   ImageSize -> 500,  
        Frame -> True, AspectRatio -> 1/3, 
        ImageSize -> 500];

DiscreteListPlot[
        InverseFourier[lpf*(Last  /@  fft)]  //  Chop,  Blue,  
        Frame -> True, Axes -> {True, False},  
        ImageSize -> 500,  AspectRatio -> 1/3];

View as "image"

tmp  =  
        InverseFourier[lpf*(Last  /@  fft)]  //  
          Chop;
ShowImage[{tmp}, Round[500/Length[tmp]], 2, -2];

Set image directory

If[StringQ[$ImagePath], $path = $ImagePath,  $path = "."]; 
    SetDirectory[$path];

2D DFT: Image example #1

Load gray-scale image

img = ReadImage["cat256x256.ppm"];
scale  =  1.5;

Display image

ShowImage[img, scale];

DFT of image

DFT of image

fft = Fourier[img];

Magnitude DFT (frequencies radiate from the center from highest to lowest)

tmp = Log[1 + Abs[fft]];
ShowImage[tmp, scale,  Max[tmp]];

Phase DFT

ShowImage[Arg[fft], scale, π, -π];

DFT of image with frequencies rearranged

DFT

shift = 
        Table[(-1)^(x + y), {x, 
            0, (Dimensions[img] // First) - 1}, {y, 
            0, (Dimensions[img] // Last) - 1 }];
fft = Fourier[img*shift];

Magnitude DFT (frequencies radiate from the center from lowest to highest)

tmp = Log[1 + Abs[fft]];
ShowImage[tmp, scale,  Max[tmp]];

Phase DFT

ShowImage[Arg[fft], scale, π, -π];

How important is phase?

Take DFT (frequencies radiate from the center from lowest to highest)

fft = Fourier[img*shift];

Eliminate magnitude information

fftphase  =  Map[(#/Abs[#]) &,  fft,  {2}];

Magnitude DFT (constant)

ShowImage[Abs[fftphase], scale, 2];

Take inverse DFT of constant-magnitude image DFT

ShowImage[Abs[InverseFourier[fftphase]] //  Chop, 
        scale,  .01];

Filtering in the frequency domain

Define 2D low-pass filter

LowPassFilter2D[dim_,  
          cf_] := 
Module[{m, c, 
            r}, 
m  = 
            ZeroMatrix[dim]; 
c = (dim)/2  //  
              N; 
r = cf*dim/2; 
Do[
            m[[x, y]] = 
              If[(x - c)^2  +  (y - c)^2  <  r^2, 1, 0], {x, 1, 
              dim}, {y, 1, 
              dim}]; 
m
];

Take DFT (frequencies radiate from the center from lowest to highest)

fft = Fourier[img*shift];

Define low-pass filter

cf =  .2;
lpf  =  
        LowPassFilter2D[Dimensions[img]  //  First, cf];

Visualize low-pass filter

ShowImage[lpf, scale, 2];

Apply filter and view resulting filtered image

tmp = Abs[InverseFourier[fft*lpf]];
ShowImage[tmp, scale];

Define high-pass filter

cf =  .2;
hpf = 
        1 - LowPassFilter2D[Dimensions[img]  //  First, 
            1 - cf];

Visualize high-pass filter

ShowImage[hpf, scale, 2];

tmp = Abs[InverseFourier[fft*hpf]];
ShowImage[tmp, scale, 30];

2D DFT: Image example #2

Load gray-scale image

img = 
        ColorToGrayScale[
          ReadImage["part128x128.ppm"]];
scale  =  3;

Display image

ShowImage[img, scale];

DFT of image

DFT

fft = Fourier[img];

Magnitude DFT (frequencies radiate from the center from highest to lowest)

tmp = Log[1 + Abs[fft]];
ShowImage[tmp, scale,  Max[tmp]];

Phase DFT

ShowImage[Arg[fft], scale, π, -π];

DFT of image with frequencies rearranged

DFT

shift = 
        Table[(-1)^(x + y), {x, 
            0, (Dimensions[img] // First) - 1}, {y, 
            0, (Dimensions[img] // Last) - 1 }];
fft = Fourier[img*shift];

Magnitude DFT (frequencies radiate from the center from lowest to highest)

tmp = Log[1 + Abs[fft]];
ShowImage[tmp, scale,  Max[tmp]];

Phase DFT

ShowImage[Arg[fft], scale, π, -π];

How important is phase?

Take DFT (frequencies radiate from the center from lowest to highest)

fft = Fourier[img*shift];

Eliminate magnitude information

fftphase  =  Map[(#/Abs[#]) &,  fft,  {2}];

Magnitude DFT (constant)

ShowImage[Abs[fftphase], scale, 2];

Take inverse DFT of constant-magnitude image DFT

ShowImage[Abs[InverseFourier[fftphase]] //  Chop, 
        scale,  .01];

Filtering in the frequency domain

Define 2D low-pass filter

LowPassFilter2D[dim_,  
          cf_] := 
Module[{m, c, 
            r}, 
m  = 
            ZeroMatrix[dim]; 
c = (dim)/2  //  
              N; 
r = cf*dim/2; 
Do[
            m[[x, y]] = 
              If[(x - c)^2  +  (y - c)^2  <  r^2, 1, 0], {x, 1, 
              dim}, {y, 1, 
              dim}]; 
m
];

Take DFT (frequencies radiate from the center from lowest to highest)

fft = Fourier[img*shift];

Define low-pass filter

cf =  .2;
lpf  =  
        LowPassFilter2D[Dimensions[img] // First, cf];

Visualize low-pass filter

ShowImage[lpf, scale, 2];

Apply filter and view resulting filtered image

tmp = Abs[InverseFourier[fft*lpf]];
ShowImage[tmp, scale];

Define high-pass filter

cf =  .2;
hpf = 
        1 - LowPassFilter2D[Dimensions[img] // First, 1 - cf];

Visualize high-pass filter

ShowImage[hpf, scale, 2];

tmp = Abs[InverseFourier[fft*hpf]];
ShowImage[tmp, scale, 30];

2D DFT: Image example #3

Load gray-scale image

img = 
        ReadImage["bakery256x256.ppm"];
scale  =  1.5;

Display image

ShowImage[img, scale];

DFT of image

fft = Fourier[img];

Magnitude DFT (frequencies radiate from the center from highest to lowest)

tmp = Log[1 + Abs[fft]];
ShowImage[tmp, scale,  Max[tmp]];

Phase DFT

ShowImage[Arg[fft], scale, π, -π];

DFT of image with frequencies rearranged

shift = 
        Table[(-1)^(x + y), {x, 
            0, (Dimensions[img] // First) - 1}, {y, 
            0, (Dimensions[img] // Last) - 1 }];
fft = Fourier[img*shift];

Magnitude DFT (frequencies radiate from the center from lowest to highest)

tmp = Log[1 + Abs[fft]];
ShowImage[tmp, scale,  Max[tmp]];

Phase DFT

ShowImage[Arg[fft], scale, π, -π];

How important is phase?

Take DFT (frequencies radiate from the center from lowest to highest)

fft = Fourier[img*shift];

Eliminate magnitude information

fftphase  =  Map[(#/Abs[#]) &,  fft,  {2}];

Magnitude DFT (constant)

ShowImage[Abs[fftphase], scale, 2];

Take inverse DFT of constant-magnitude image DFT

ShowImage[Abs[InverseFourier[fftphase]] //  Chop, 
        scale,  .01];

Filtering in the frequency domain

Define 2D low-pass filter

LowPassFilter2D[dim_,  
          cf_] := 
Module[{m, c, 
            r}, 
m  = 
            ZeroMatrix[dim]; 
c = (dim)/2  //  
              N; 
r = cf*dim/2; 
Do[
            m[[x, y]] = 
              If[(x - c)^2  +  (y - c)^2  <  r^2, 1, 0], {x, 1, 
              dim}, {y, 1, 
              dim}]; 
m
];

Take DFT (frequencies radiate from the center from lowest to highest)

fft = Fourier[img*shift];

Define low-pass filter

cf =  .2;
lpf  =  
        LowPassFilter2D[Dimensions[img] // First, cf];

Visualize low-pass filter

ShowImage[lpf, scale, 2];

Apply filter and view resulting filtered image

tmp = Abs[InverseFourier[fft*lpf]];
ShowImage[tmp, scale];

Define high-pass filter

cf =  .2;
hpf = 
        1 - LowPassFilter2D[Dimensions[img] // First, 1 - cf];

Visualize high-pass filter

ShowImage[hpf, scale, 2];

tmp = Abs[InverseFourier[fft*hpf]];
ShowImage[tmp, scale, 30];

Key shortcoming Fourier analysis: poor localization of frequencies

Example #1: sum of sinusoids

Continuous-time waveform

x1[t_] := 
        Cos[2  π 10  t]  +  Cos[2  π 25 t]  +  
          Cos[2  π 50 t]  +  
          Cos[2  π 100 t];
gx1 = 
        Plot[x1[t], {t, 0, 1},  PlotPoints -> 100, 
          PlotStyle -> {Blue},  
          ImageSize -> 500];
x1[t]  //  TraditionalForm

Discrete Fourier Transform (DFT)

sampleRate = 1000;

signal  =  
        Table[x1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  
          N;

fft  =  Fourier[signal];
gf1 = 
        ListPlot[FourierRearrange[Abs[fft], sampleRate], 
          PlotStyle -> {Red}, 
          PlotRange -> {{-125, 125}, All},  Frame -> True,  
          ImageSize -> 500];

Example #2: time-varying sinusoid

Continuous-time waveform

u[t_] := UnitStep[t];

x2[
          t_] := (u[t] - u[t -  .3]) 
            Cos[2  π 10  t]  +  (u[t -  .3] - u[t -  .6]) 
            Cos[2  π 25 t]  + (u[t -  .6] - u[t -  .8]) Cos[
              2  π 50 t]  +  (u[t -  .8] - u[t - 1. ]) 
            Cos[2  π 100 t];
gx2 = 
        Plot[x2[t], {t, 0, 1},  PlotPoints -> 100, 
          PlotStyle -> {Blue},  
          ImageSize -> 500];
x2[t]  //  TraditionalForm

Fast Fourier Transform (FFT)

sampleRate = 1000;

signal  =  
        Table[x2[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  
          N;

fft  =  Fourier[signal];
gf2 = 
        ListPlot[FourierRearrange[Abs[fft], sampleRate], 
          PlotStyle -> {Red}, 
          PlotRange -> {{-125, 125}, All},  Frame -> True,  
          ImageSize -> 500];

Display all

Show[GraphicsArray[{{gx1, gx2}, {gf1, gf2}}], 
        ImageSize -> 640];