Difference equations: examples and applications

Mathematica definitions

Load needed packages for this notebook

<< Graphics`Colors`;
<< Miscellaneous`Audio`;

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}];

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]
];

Filtering definitions

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
];

Filter[n, sampleRate, {lowFreq, highFreq}] generates an ideal discrete filter for an FFT list of length 'n' (sampled at 'sampleRate'), with a bandpass from 'lowFreq' to 'highFreq' (in Hz).

Filter[n_, 
          samplingRate_, {lowF_, 
            highF_}] := 
Module[{coeff, l, 
            h}, 
coeff = 
            Table[0, {n}]; 
l = Floor[lowF*n/samplingRate]; 
          h = Ceiling[highF*n/samplingRate]; 
If[
            l  ==  0,  l  =  1;  
            coeff = ReplacePart[coeff, 1, 1]]; 
coeff = 
            ReplacePart[coeff, 1, Table[{i}, {i, l + 1, h + 1}]]; 
          coeff = ReplacePart[coeff, 1, 
              Table[{i}, {i, n + 1 - h, 
                  n + 1 - l}]]
];

FrequencyResponse[h_,  ωn_] := 
      h . Table[Exp[-I ωn k], {k, 0, Length[h] - 1}]

Non-recursive difference equation example: simple low-pass filter

y[n]  =  
        1/L (x[n]  +  x[n - 1]  +  x[n - 2]  +   ...   +  
              x[n - L + 1])
 y[n]  =  
        1/L (∑_(k = 0)^(L - 1))x[n - k]

Noisy, discrete-time waveform

sampleRate = 500;

f1[t_] := Sin[2  π 4  t]  +  Random[] - 1/2;
signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
g1 = ListPlot[signal,  PlotStyle -> {Blue}];

Filter #1

Time domain

L = 2;

h = Table[1/L  //  N, {L}];
filtered  =  
        ListCorrelate[h,  signal,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Thickness[ .01], Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Frequency domain

plrng  = {0, 40};
fftSignal  =  Fourier[signal];
gf1 = 
        DiscreteListPlot[FourierRearrange[Abs[fftSignal], sampleRate], Blue, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
fftFiltered  =  Fourier[filtered];
gf2 = 
        DiscreteListPlot[FourierRearrange[Abs[fftFiltered], sampleRate], Red, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{gf1, gf2}], ImageSize -> 640];

Filter frequency response

Plot[
        Abs[FrequencyResponse[h, 
            2  π ωn/sampleRate]], {ωn, -sampleRate/2, 
          sampleRate/2},  
        PlotStyle -> {Red,  Thickness[ .01]}];

Filter #2

Time domain

L = 10;

h = Table[1/L  //  N, {L}];
filtered  =  
        ListCorrelate[h,  signal,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Thickness[ .01], Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Frequency domain

plrng  = {0, 40};
fftSignal  =  Fourier[signal];
gf1 = 
        DiscreteListPlot[FourierRearrange[Abs[fftSignal], sampleRate], Blue, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
fftFiltered  =  Fourier[filtered];
gf2 = 
        DiscreteListPlot[FourierRearrange[Abs[fftFiltered], sampleRate], Red, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{gf1, gf2}], ImageSize -> 640];

Filter frequency response

Plot[
        Abs[FrequencyResponse[h, 
            2  π ωn/sampleRate]], {ωn, -sampleRate/2, 
          sampleRate/2},  PlotStyle -> {Red,  Thickness[ .01]},  
        PlotRange -> All];

Filter #3

Time domain

L = 100;

h = Table[1/L  //  N, {L}];
filtered  =  
        ListCorrelate[h,  signal,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Thickness[ .01], Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Show[g1, g2];

Frequency domain

plrng  = {0, 40};
fftSignal  =  Fourier[signal];
gf1 = 
        DiscreteListPlot[FourierRearrange[Abs[fftSignal], sampleRate], Blue, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
fftFiltered  =  Fourier[filtered];
gf2 = 
        DiscreteListPlot[FourierRearrange[Abs[fftFiltered], sampleRate], Red, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{gf1, gf2}], ImageSize -> 640];

Filter frequency response

Plot[
        Abs[FrequencyResponse[h, 
            2  π ωn/sampleRate]], {ωn, -sampleRate/2, 
          sampleRate/2},  PlotStyle -> {Red},  PlotRange -> All, 
        PlotPoints -> 100];

Pass through single frequency

sampleRate = 500;

f1[t_] := Sin[2  π 20  t];
signal2  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
g1 = ListPlot[signal2,  PlotStyle -> {Blue}];

L = 100;

h = Table[1/L  //  N, {L}];
filtered  =  
        ListCorrelate[h,  signal2,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Non-recursive difference equation example: simple high-pass filter #1

y[n]  =  
        1/2  x[n]  -  1/2 x[n - 1]

Noisy, discrete-time waveform

sampleRate = 500;

f1[t_] := Sin[2  π 4  t]  +  Random[] - 1/2;
signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
g1 = ListPlot[signal,  PlotStyle -> {Blue}];

Filter

Time domain

h = {1/2, -1/2};
filtered  =  
        ListCorrelate[h,  signal,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Frequency domain

plrng  = {0, 40};
fftSignal  =  Fourier[signal];
gf1 = 
        DiscreteListPlot[FourierRearrange[Abs[fftSignal], sampleRate], Blue, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
fftFiltered  =  Fourier[filtered];
gf2 = 
        DiscreteListPlot[FourierRearrange[Abs[fftFiltered], sampleRate], Red, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{gf1, gf2}], ImageSize -> 640];

Filter frequency response

Plot[
        Abs[FrequencyResponse[h, 
            2  π ωn/sampleRate]], {ωn, -sampleRate/2, 
          sampleRate/2},  
        PlotStyle -> {Red,  Thickness[ .01]}];

Non-recursive difference equation example: simple high-pass filter #2

y[n]  =  
        1/4  x[n]  -  
          1/2  x[n - 1]  +  1/4 x[n - 2]

Noisy, discrete-time waveform

sampleRate = 500;

f1[t_] := Sin[2  π 4  t]  +  Random[] - 1/2;
signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
g1 = ListPlot[signal,  PlotStyle -> {Blue}];

Filter

Time domain

h = {1/4, -1/2, 1/4};
filtered  =  
        ListCorrelate[h,  signal,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Frequency domain

plrng  = {0, 40};
fftSignal  =  Fourier[signal];
gf1 = 
        DiscreteListPlot[FourierRearrange[Abs[fftSignal], sampleRate], Blue, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
fftFiltered  =  Fourier[filtered];
gf2 = 
        DiscreteListPlot[FourierRearrange[Abs[fftFiltered], sampleRate], Red, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{gf1, gf2}], ImageSize -> 640];

Filter frequency response

Plot[
        Abs[FrequencyResponse[h, 
            2  π ωn/sampleRate]], {ωn, -sampleRate/2, 
          sampleRate/2},  
        PlotStyle -> {Red,  Thickness[ .01]}];

Non-recursive difference equation example: simple high-pass filter #3

y[n]  =  
        1/8  x[n]  -  3/8  x[n - 1]  +  
          3/8  x[n - 2] - 1/8 x[n - 3]

Noisy, discrete-time waveform

sampleRate = 500;

f1[t_] := Sin[2  π 4  t]  +  Random[] - 1/2;
signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
g1 = ListPlot[signal,  PlotStyle -> {Blue}];

Filter

Time domain

h = {1/8, -3/8, 3/8, -1/8};
filtered  =  
        ListCorrelate[h,  signal,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Frequency domain

plrng  = {0, 40};
fftSignal  =  Fourier[signal];
gf1 = 
        DiscreteListPlot[FourierRearrange[Abs[fftSignal], sampleRate], Blue, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
fftFiltered  =  Fourier[filtered];
gf2 = 
        DiscreteListPlot[FourierRearrange[Abs[fftFiltered], sampleRate], Red, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{gf1, gf2}], ImageSize -> 640];

Filter frequency response

Plot[
        Abs[FrequencyResponse[h, 
            2  π ωn/sampleRate]], {ωn, -sampleRate/2, 
          sampleRate/2},  
        PlotStyle -> {Red,  Thickness[ .01]}];

Non-recursive difference equation example: simple high-pass filter #4

y[n]  =  
        1/16  x[n]  -  1/4  x[n - 1]  +  3/8  x[n - 2] - 
          1/4  x[n - 3] + 1/16 x[n - 4]

Noisy, discrete-time waveform

sampleRate = 500;

f1[t_] := Sin[2  π 4  t]  +  Random[] - 1/2;
signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
g1 = ListPlot[signal,  PlotStyle -> {Blue}];

Filter

Time domain

h = {1/16, -1/4, 3/8, -1/4, 
          1/16};
filtered  =  
        ListCorrelate[h,  signal,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Frequency domain

plrng  = {0, 40};
fftSignal  =  Fourier[signal];
gf1 = 
        DiscreteListPlot[FourierRearrange[Abs[fftSignal], sampleRate], Blue, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
fftFiltered  =  Fourier[filtered];
gf2 = 
        DiscreteListPlot[FourierRearrange[Abs[fftFiltered], sampleRate], Red, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{gf1, gf2}], ImageSize -> 640];

Filter frequency response

Plot[
        Abs[FrequencyResponse[h, 
            2  π ωn/sampleRate]], {ωn, -sampleRate/2, 
          sampleRate/2},  
        PlotStyle -> {Red,  Thickness[ .01]}];

Non-recursive difference equation example: well-designed band-pass filter

y[n]  =  -0.0193 x[n]  +  
          0.0099 x[n - 1]  -  0.0003 x[n - 2]  +  0.0276 x[n - 3] + 
          0.0000 x[n - 4] - 0.0649 x[n - 5]  +  0.0264 x[n - 6] - 
          0.0126  x[n - 7] + 0.1188  x[n - 8] + 0.0000  x[n - 9] - 
          0.3331  x[n - 10] + 0.2422  x[n - 11] + 0.2422  x[n - 12] - 
          0.3331  x[n - 13] + 0.0000  x[n - 14] + 0.1188  x[n - 15] - 
          0.0126  x[n - 16] + 0.0264  x[n - 17] - 0.0649  x[n - 18] + 
          0.0000  x[n - 19] + 0.0276  x[n - 20] - 0.0003  x[n - 21] + 
          0.0099  x[n - 22] - 0.0193  x[n - 23]
y[n] = (∑_(k = 0)^23) (Subscript[b, k]) x[n - k]

Noisy, discrete-time waveform

sampleRate = 500;

f1[t_] := Sin[2  π 4  t]  +  Random[] - 1/2;
signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
g1 = ListPlot[signal,  PlotStyle -> {Blue}];

Filter

Time domain

h1 = {-0.0193, 0.0099, -0.0003, 0.0276, 0.0000, -0.0649, 
          0.0264, -0.0126, 0.1188, 0.0000, -0.3331, 
          0.2422};
h  =  h1~Join~Reverse[h1];
filtered  =  
        ListCorrelate[h,  signal,  {-1, -1}, 
          0];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

Frequency domain

plrng  = {0, 40};
fftSignal  =  Fourier[signal];
gf1 = 
        DiscreteListPlot[FourierRearrange[Abs[fftSignal], sampleRate], Blue, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
fftFiltered  =  Fourier[filtered];
gf2 = 
        DiscreteListPlot[FourierRearrange[Abs[fftFiltered], sampleRate], Red, 
          PlotRange -> plrng,  Frame -> True, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{gf1, gf2}], ImageSize -> 640];

Filter frequency response

Plot[
        Abs[FrequencyResponse[h, 
            2  π ωn/sampleRate]], {ωn, -sampleRate/2, 
          sampleRate/2},  
        PlotStyle -> {Red,  Thickness[ .01]}];

Simple non-recursive filters on speech example

Import wav file

Adjust the directory setting and input file below to your system configuration.

dir = "/Users/nechyba/public_html/eel3135/data";
file = "Speech.wav";

SetDirectory[dir];
wav1 = ReadSoundFile[file, PrintHeader -> True];

sampleRate1 = 11127;

Original wav file

ListPlay[wav1, SampleRate -> sampleRate1];

Click on the sound plot above to hear the speech signal.

Simple low-pass filter

y[n]  =  
        1/L (x[n]  +  x[n - 1]  +  x[n - 2]  +   ...   +  
              x[n - L + 1])
 y[n]  =  
        1/L (∑_(k = 0)^(L - 1))x[n - k]

L = 100;

filtered  =  
        ListCorrelate[Table[1/L  //  N,  {L}],  wav1,  {-1, -1}, 
          0];
ListPlay[filtered,  SampleRate -> sampleRate1];

Simple high-pass filter

y[n]  =  
        1/4  x[n]  -  
          1/2  x[n - 1]  +  1/4 x[n - 2]

filtered  =  
        ListCorrelate[{1/4, -1/2, 1/4},  wav1,  {-1, -1}, 
          0];
ListPlay[filtered,  SampleRate -> sampleRate1];

Simple non-recursive filters in two dimensions (image)

Import gray-scale image

Adjust the directory setting and input file below to your system configuration.

dir = "/Users/nechyba/public_html/eel3135/data";
file = "sample_image.jpg";

SetDirectory[dir];
im1 = Import[file]; 
g1 = Show[im1,  ImageSize -> 500];

Low-pass running-average filter

L = 20;

{imtst, tmp1, tmp2, tmp3}  =  List  @@  im1[[1]];
mask  =  Table[1, {L}, {L}];
imout1 = ListConvolve[mask, imtst];
range  =  {Min[Flatten[imout1]], Max[Flatten[imout1]]};
g2 = 
        Show[Graphics[Raster[imout1, tmp1, range, tmp3]], 
          AspectRatio -> Automatic, 
          ImageSize -> 500];

Show[GraphicsArray[{g1, g2}], ImageSize -> 640];

High-pass filter

{imtst, tmp1, tmp2, tmp3}  =  List  @@  im1[[1]];
mask  = {{1, 1, 1}, {1, -8, 1}, {1, 1, 1}};
imout1 = Abs[ListConvolve[mask, imtst]];
range  =  {Min[Flatten[imout1]], Max[Flatten[imout1]]};
g2 = 
        Show[Graphics[Raster[imout1, tmp1, range, tmp3]], 
          AspectRatio -> Automatic, 
          ImageSize -> 500];

Show[GraphicsArray[{g1, g2}], ImageSize -> 640];