Spectral leackage and windowing

Mathematica definitions

Load needed packages for this notebook

<< Graphics`Colors`;
<< Calculus`FourierTransform`;
<< Graphics`Graphics`;

Set up defaults

Turn off annoying warning messages

Off[General::spell];
Off[General::spell1];

Set default options for specific functions

SetOptions[Plot, Frame -> True];
SetOptions[ListPlot, Frame -> True,  PlotJoined -> True,  
        Axes -> {True, False}];
SetOptions[Fourier, 
        FourierParameters -> {1, -1}];
SetOptions[InverseFourier, 
        FourierParameters -> {1, -1}];
SetOptions[FourierTransform,  
        FourierParameters -> {0, -2 π}];
SetOptions[FourierTransform,  
        FourierParameters -> {0, -2 π}];
SetOptions[InverseFourierTransform,  
        FourierParameters -> {0, -2 π}];

Visualization functions

Fix time variable for sampled signals

FixTime[signal_, 
          sampleRate_] := 
Module[{time, 
            data}, 
time = 
            Table[i/sampleRate, {i, 1, 
                Length[signal]}]; 
data  =  {time, 
                signal} // Transpose
];

Generate periodic waveform from one period

PeriodicSignal[f, {dmin, dmax}, {t, min, max}] generates a periodic waveform for the interval [min, max]. It is assumed that one period of the function f(t) is defined over the interval [dmin, dmax].

PeriodicSignal[
          f_, {dmin_, dmax_}, {t_, min_, 
            max_}] := 
Module[{T, imin, 
            imax}, 
T  =  
            dmax  -  dmin; 
imin = 
            Floor[(min - dmin)/T]; 
imax = 
            Ceiling[(max - dmax)/T];  
Sum[
            f[t - i T] (UnitStep[t - dmin - i T]  -  
                  UnitStep[t - dmax - i T]), {i, imin, 
              imax}]
];

Plot periodic waveforms

PeriodicPlot[f, {dmin, dmax}, {t, min, max}, opts] plots a a periodic waveform for the interval [min, max]. It is assumed that one period of the function f(t) is defined over the interval [dmin, dmax]. The 'opts' argument passed plotting options to the Mathematica function Plot[ ].

PeriodicPlot[f_, {dmin_, dmax_}, {t_, min_, max_}, 
          opts___] := 
Module[{fnc}, 
fnc = PeriodicSignal[f, {dmin, dmax}, {t, min, max}]; 
Plot[
            fnc, {t, min, max}, opts]
];

Plot discrete-time signal as sequence of amplitude lines

DiscreteListPlot[signal_, color_:  Black, 
          opts___] := 
Module[{data, g1, 
            g2}, 
If[
            Length[signal[[1]]]  ==  2, 
data = 
              signal, 
data = 
              FixTime[signal, 
                1]
]; 
g1 = Line[{{#[[1]], 0}, {#[[1]], #[[2]]}}] &  /@  
              data; 
g2 = {PointSize[ .015], 
                  Point[#]} &  /@  data; 
Show[
            Graphics[{color, g1, g2}], 
            opts]
];

Plot discrete-time signal as sequence of amplitude lines

DiscreteListPlot2[signal_, color_:  Black, 
          opts___] := 
Module[{data, g1, 
            g2}, 
If[
            Length[signal[[1]]]  ==  2, 
data = 
              signal, 
data = 
              FixTime[signal, 
                1]
]; 
g1 = Line[{{#[[1]], 0}, {#[[1]], #[[2]]}}] &  /@  
              data; 
g2 = {PointSize[ .005], 
                  Point[#]} &  /@  data; 
Show[
            Graphics[{color, g1, g2}], 
            opts]
];

Fourier rearrange

FourierRearrange[fft] rearranges the result of an FFT (Mathematica function Fourier[ ]), to a easier to visualize the spectrum represented by 'fft'. The length of the list 'fft' is assumed to be even; if it is not, the last term in 'fft' is dropped.

FourierRearrange[fft_,  
          samplingRate_] := 
Module[{fft2, n, 
            n2}, 
If[EvenQ[Length[fft]],  
            fft2  =  fft,  
            fft2  =  Drop[fft, -1]];  
n = 
            Length[fft]; 
          n2 = n/2; 
{Table[ 
                samplingRate i/n, {i, -n2 + 1, n2}], 
              Take[fft2, {n2 + 2, n2*2}]~Join~Take[fft2, {1, n2 + 1}]} // 
            Transpose
];

Short-cut functions

Dirac delta function

δ[t_] := DiracDelta[t];

Unit step function

u[t_] := UnitStep[t];

Pulse of width w

p[t_, w_] := u[t + w/2] - u[t - w/2];

Hack to plot δ(t)

PlotConvert[X_, 
        w_] := 
X  /.  {δ[a_] -> 
            p[a, w]}

Sampler (CT to DT converter)

sampler[xc_, fs_, 
          range_: {0, 5}] := 
Table[{n, xc[n/fs]}, {n, 
            range[[1]]*fs, range[[2]]*fs}];

Sampled frequency spectrum

SampledSpectrum[X_, f_, fs_, n_:  5] := 
Sum[
          X[f - k fs], {k, -n, n}];

Definition of DTFT

DTFT[
          xd_, θ_] := 
Module[{n, xn, n1, 
            n2}, 
{n, xn}  =  
            xd  //  Transpose; 
n1  =  First[n];  
          n2  =  Last[n]; 
Sum[
            xn[[k - n1 + 1]] Exp[-I k θ], {k, n1, 
              n2}]
];

Definition of DFT

DFT[
          x_] := 
Module[{x2,  k}, 
If[
            Length[First[x]]  >  1,  x2  =  Last  /@  x,  
            x2  =  x]; 
{Table[
                k, {k, 0, Length[x2] - 1}], Fourier[x2]}  //  
            Transpose];

Definition of spectral leakage window

Window[t_,  
        interval_] := (0.54  -   .46* 
            Cos[2  π/interval t])

Windowing example #1 (single cosine wave  (2+1/6) Hz)

Unwindowed example w/zero padding

x[t_] := Cos[(4  π  +  2  π (1/6)) t];
xda =  sampler[x, fs = 20, {0, 3}];
xd2 =  
        Table[{n, 0}, {n, -fs, -1}]~Join~xda~Join~ 
          Table[{n, 0}, {n, 3*fs + 1, 4*fs}];
 If[OddQ[Length[xda]],  xda  =  Drop[xda, -1]];
gwa  =  
        DiscreteListPlot2[xd2,  Blue,  Frame -> True, 
          Axes -> {True, False},  AspectRatio -> 1/3,  
          ImageSize -> 504,  PlotRange -> All];

Plot windowing function

x[t_] := Window[t,  3];
xdb =  sampler[x, fs, {0, 3}];
xd2 =  
        Table[{n, 0}, {n, -fs, -1}]~Join~xdb~Join~ 
          Table[{n, 0}, {n, 3*fs + 1, 4*fs}];
 If[OddQ[Length[xdb]],  xdb  =  Drop[xdb, -1]];
gwb  =  
        DiscreteListPlot2[xd2,  Blue,  Frame -> True, 
          Axes -> {True, False},  AspectRatio -> 1/3,  
          ImageSize -> 504,  PlotRange -> All];

Sampled and windowed sequence

x[t_] :=  
        Cos[(4  π  +  2  π (1/6)) t] 
          Window[t, 3];
xdc =  sampler[x, fs, {0, 3}];
xd2 =  
        Table[{n, 0}, {n, -fs, -1}]~Join~xdc~Join~ 
          Table[{n, 0}, {n, 3*fs + 1, 4*fs}];
 If[OddQ[Length[xdc]],  xdc  =  Drop[xdc, -1]];
gwc  =  
        DiscreteListPlot2[xd2,  Blue,  Frame -> True, 
          Axes -> {True, False},  AspectRatio -> 1/3,  
          ImageSize -> 504,  PlotRange -> All];

Show together

Show[GraphicsArray[{{gwa}, {gwb}, {gwc}}],  
        ImageSize -> 504];

Magnitude plot of unwindowed DTFT as a function of frequency

Xθ  = DTFT[xda, θ]  //  Simplify; 
    g1 = Plot[
        Abs[Xθ  /.  θ ->  2  π f/fs], {f, -fs/
            2, fs/2},  PlotStyle -> {Blue,  Dashing[{ .01}]},  
        PlotRange -> All, PlotPoints -> 200];

Magnitude plot of DFT as a function of frequency

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

Rearrange and scale DFT result

g5 = 
        DiscreteListPlot[FourierRearrange[Abs[Last  /@  fft], fs], Red, 
          PlotRange -> All, ImageSize -> 504,  
          Frame -> True,  Axes -> False,  
          DisplayFunction -> Identity];
g6  =  Show[g1, g5];

Show both

gw1 = 
        Show[GraphicsArray[{{g1, g5,  g6}}],  
          ImageSize -> 504*1.5];

Magnitude plot of windowed DTFT as a function of frequency

Xθ  = DTFT[xdc, θ]  //  Simplify; 
    g1 = Plot[
        Abs[Xθ  /.  θ ->  2  π f/fs], {f, -fs/
            2, fs/2},  PlotStyle -> {Blue,  Dashing[{ .01}]},  
        PlotRange -> All, PlotPoints -> 200];

Magnitude plot of DFT as a function of frequency

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

Rearrange and scale DFT result

g5 = 
        DiscreteListPlot[FourierRearrange[Abs[Last  /@  fft], fs], Red, 
          PlotRange -> All, ImageSize -> 504,  
          Frame -> True,  Axes -> False,  
          DisplayFunction -> Identity];
g6  =  Show[g1, g5];

Show both

gw2 = 
        Show[GraphicsArray[{{g1, g5,  g6}}],  
          ImageSize -> 504*1.5];

Show all

Show[gw1];
Show[gw2];

Windowing example #2 (example #6 waveform, DC, (1+1/6) Hz, (2+1/6)Hz)

Unwindowed example w/zero padding

x[
          t_] := (1  +  2 Cos[2 (1 + 1/6) π t]  +  
            4 Cos[(4  π  +  2  π (1/6)) 
                  t]);
xda =  sampler[x, fs = 20, {0, 3}];
xd2 =  
        Table[{n, 0}, {n, -fs, -1}]~Join~xda~Join~ 
          Table[{n, 0}, {n, 3*fs + 1, 4*fs}];
 If[OddQ[Length[xda]],  xda  =  Drop[xda, -1]];
gwa  =  
        DiscreteListPlot2[xd2,  Blue,  Frame -> True, 
          Axes -> {True, False},  AspectRatio -> 1/3,  
          ImageSize -> 504,  PlotRange -> All];

Plot windowing function

x[t_] := Window[t,  3];
xdb =  sampler[x, fs, {0, 3}];
xd2 =  
        Table[{n, 0}, {n, -fs, -1}]~Join~xdb~Join~ 
          Table[{n, 0}, {n, 3*fs + 1, 4*fs}];
 If[OddQ[Length[xdb]],  xdb  =  Drop[xdb, -1]];
gwb  =  
        DiscreteListPlot2[xd2,  Blue,  Frame -> True, 
          Axes -> {True, False},  AspectRatio -> 1/3,  
          ImageSize -> 504,  PlotRange -> All];

Sampled and windowed sequence

x[
          t_] := (1  +  2 Cos[2 (1 + 1/6) π t]  +  
              4 Cos[(4  π  +  2  π (1/6)) t]) 
          Window[t, 3];
xdc =  sampler[x, fs, {0, 3}];
xd2 =  
        Table[{n, 0}, {n, -fs, -1}]~Join~xdc~Join~ 
          Table[{n, 0}, {n, 3*fs + 1, 4*fs}];
 If[OddQ[Length[xdc]],  xdc  =  Drop[xdc, -1]];
gwc  =  
        DiscreteListPlot2[xd2,  Blue,  Frame -> True, 
          Axes -> {True, False},  AspectRatio -> 1/3,  
          ImageSize -> 504,  PlotRange -> All];

Show together

Show[GraphicsArray[{{gwa}, {gwb}, {gwc}}],  
        ImageSize -> 504];

Magnitude plot of unwindowed DTFT as a function of frequency

Xθ  = DTFT[xda, θ]  //  Simplify; 
    g1 = Plot[
        Abs[Xθ  /.  θ ->  2  π f/fs], {f, -fs/
            2, fs/2},  PlotStyle -> {Blue,  Dashing[{ .01}]},  
        PlotRange -> All, PlotPoints -> 200];

Magnitude plot of DFT as a function of frequency

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

Rearrange and scale DFT result

g5 = 
        DiscreteListPlot[FourierRearrange[Abs[Last  /@  fft], fs], Red, 
          PlotRange -> All, ImageSize -> 504,  
          Frame -> True,  Axes -> False,  
          DisplayFunction -> Identity];
g6  =  Show[g1, g5];

Show both

gw1 = 
        Show[GraphicsArray[{{g1, g5,  g6}}],  
          ImageSize -> 504*1.5];

Magnitude plot of windowed DTFT as a function of frequency

Xθ  = DTFT[xdc, θ]  //  Simplify; 
    g1 = Plot[
        Abs[Xθ  /.  θ ->  2  π f/fs], {f, -fs/
            2, fs/2},  PlotStyle -> {Blue,  Dashing[{ .01}]},  
        PlotRange -> All, PlotPoints -> 200];

Magnitude plot of DFT as a function of frequency

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

Rearrange and scale DFT result

g5 = 
        DiscreteListPlot[FourierRearrange[Abs[Last  /@  fft], fs], Red, 
          PlotRange -> All, ImageSize -> 504,  
          Frame -> True,  Axes -> False,  
          DisplayFunction -> Identity];
g6  =  Show[g1, g5];

Show both

gw2 = 
        Show[GraphicsArray[{{g1, g5,  g6}}],  
          ImageSize -> 504*1.5];

Show all

Show[gw1];
Show[gw2];