Recursive difference equations

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

Linear, recursive difference equation example

y[
          n]  =  α y[n - 1]  +  (1 - α) x[n]

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}];
plr = AbsoluteOptions[g1,  PlotRange]  //  First  //  
          Last;

x[n_] := 
If[
          n ≤ 0 || 
            n > Length[
                signal], 
0, 
signal[[n]]
];

Experiment #1

Difference equation definition

Remove[Global`y];
y[
          n_] := y[
            n] = α y[n - 1]  +  (1 - α) 
              x[n];

Initial condition and constant α

y[0] = 0;
α = 0.1;

Time domain

filtered  = 
        Table[y[n], {n, 1,  Length[signal]}];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity,  
          PlotRange -> plr];
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

Remove[Global`y];
Remove[Global`n];

eq  = 
        y[n] == α y[n - 1]  +  (1 - α) 
              DiscreteDelta[n];
y[-1] = 0;

sol = 
        Solve[ZTransform[eq,  n,  z],  ZTransform[y[n], n, z]]  //  
          Flatten;
impulseTransferFunction  =  (ZTransform[y[n], n, z]  //.  
              sol~Join~{y[-3] -> 0, y[-2] -> 0,  
                  y[-1] -> 0})  //  Simplify; 
freqResp = impulseTransferFunction  /.  
            z -> Exp[I θ]  /.  θ -> 
            2  π ωn/sampleRate;

g3 = 
        Plot[Abs[freqResp], {ωn, -sampleRate/2, sampleRate/2},  
          PlotStyle -> {Red,  Thickness[ .01]},  
          PlotRange -> {- .01, 1.01}];

Experiment #2

Difference equation definition

Remove[Global`y];
y[
          n_] := y[
            n] = α y[n - 1]  +  (1 - α) 
              x[n];

Initial condition and constant α

y[0] = 0;
α = 0.5;

Time domain

filtered  = 
        Table[y[n], {n, 1,  Length[signal]}];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity,  
          PlotRange -> plr];
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

Remove[Global`y];
Remove[Global`n];

eq  = 
        y[n] == α y[n - 1]  +  (1 - α) 
              DiscreteDelta[n];
y[-1] = 0;

sol = 
        Solve[ZTransform[eq,  n,  z],  ZTransform[y[n], n, z]]  //  
          Flatten;
impulseTransferFunction  =  (ZTransform[y[n], n, z]  //.  
              sol~Join~{y[-3] -> 0, y[-2] -> 0,  
                  y[-1] -> 0})  //  Simplify; 
freqResp = impulseTransferFunction  /.  
            z -> Exp[I θ]  /.  θ -> 
            2  π ωn/sampleRate;

g3 = 
        Plot[Abs[freqResp], {ωn, -sampleRate/2, sampleRate/2},  
          PlotStyle -> {Red,  Thickness[ .01]},  
          PlotRange -> {- .01, 1.01}];

Experiment #3

Difference equation definition

Remove[Global`y];
y[
          n_] := y[
            n] = α y[n - 1]  +  (1 - α) 
              x[n];

Initial condition and constant α

y[0] = 0;
α = 0.9;

Time domain

filtered  = 
        Table[y[n], {n, 1,  Length[signal]}];
g2  =  
        ListPlot[filtered, PlotStyle -> {Red},  
          DisplayFunction -> Identity,  
          PlotRange -> plr];
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

Remove[Global`y];
Remove[Global`n];

eq  = 
        y[n] == α y[n - 1]  +  (1 - α) 
              DiscreteDelta[n];
y[-1] = 0;

sol = 
        Solve[ZTransform[eq,  n,  z],  ZTransform[y[n], n, z]]  //  
          Flatten;
impulseTransferFunction  =  (ZTransform[y[n], n, z]  //.  
              sol~Join~{y[-3] -> 0, y[-2] -> 0,  
                  y[-1] -> 0})  //  Simplify; 
freqResp = impulseTransferFunction  /.  
            z -> Exp[I θ]  /.  θ -> 
            2  π ωn/sampleRate;

g3 = 
        Plot[Abs[freqResp], {ωn, -sampleRate/2, sampleRate/2},  
          PlotStyle -> {Red,  Thickness[ .01]},  
          PlotRange -> {- .01, 1.01}];

Unstable recursive equation

Clear memory

Remove[Global`y];
Remove[Global`x];
Remove[Global`n];

Define difference equation and initial conditions

y[
          n_] := y[n] = -1/4  y[n - 1] + 1/3  y[n - 2] - 
            1/2  y[n - 3]  +  DiscreteDelta[n];

y[-1] = y[-2] = y[-3] = 0;

Simulated output for impulse input

Input

input = 
        Table[{n, DiscreteDelta[n]}, {n, 0, 
            100}];
DiscreteListPlot[input, Blue,  Frame -> True,  
        Axes -> {True, False},  
        PlotRange -> {- .1, 2.1}];

Output

output = Table[{n, y[n]}, {n, 0, 100}];
DiscreteListPlot[output, Red, PlotRange -> All,  
        Frame -> True,  Axes -> {True, False}];

Theoretical solution

Compute z-transform and transfer function

Remove[Global`y];
Remove[Global`x];
Remove[Global`n];

eq = (y[n] == -1/4  y[n - 1] + 1/3  y[n - 2] - 
              1/2  y[n - 3]  +  DiscreteDelta[n]);

y[-1] = y[-2] = y[-3] = 0;

sol = 
        Solve[ZTransform[eq,  n,  z],  ZTransform[y[n], n, z]]  //  
          Flatten;

impulseTransferFunction  =  (ZTransform[y[n], n, z]  //.  
              sol~Join~{y[-3] -> 0, y[-2] -> 0,  
                  y[-1] -> 0})  //  
          Simplify; 
impulseTransferFunction  //  TraditionalForm

Compute poles of impulse response transfer function

poles = Solve[Denominator[impulseTransferFunction]  ==  0,  
            z]  //  Flatten  //  N ; 
    poles //  TraditionalForm

Plot pole locations

Show[
        Graphics[{{Red, Circle[{0, 0}, 1]}, {Orange, 
              PointSize[ .04], (Point[{Re[z], 
                          Im[z]}]  /.  #) &  /@  poles}}],  
        AspectRatio -> Automatic, Frame -> True,  
        Axes -> True,  
        PlotLabel -> 
          StyleForm["Complex plane",  FontFamily -> Times,  
            FontSize -> 18]];

Analytic impulse response

impulseResponse = 
        InverseZTransform[impulseTransferFunction,  z,  n];
DiscreteListPlot[Table[{n, impulseResponse},  {n, 0, 100}], Red, 
        PlotRange -> All,  Frame -> True,  
        Axes -> {True, False}];

Marginally stable recursive equation

Clear memory

Remove[Global`y];
Remove[Global`x];
Remove[Global`n];

Define difference equation and initial conditions

y[
          n_] := y[n] = -1/2  y[n - 1] + 1/2  y[n - 2]  +  
            DiscreteDelta[n];

y[-1] = y[-2] = y[-3] = 0;

Simulated output for impulse input

Input

input = 
        Table[{n, DiscreteDelta[n]}, {n, 0, 
            100}];
DiscreteListPlot[input, Blue,  Frame -> True,  
        Axes -> {True, False},  
        PlotRange -> {- .1, 2.1}];

Output

output = Table[{n, y[n]}, {n, 0, 100}];
DiscreteListPlot[output, Red, PlotRange -> All,  
        Frame -> True,  Axes -> {True, False}];

Theoretical solution

Compute z-transform and transfer function

Remove[Global`y];
Remove[Global`x];
Remove[Global`n];

eq = (y[n] == -1/2  y[n - 1] + 1/2  y[n - 2] + 
              DiscreteDelta[n]);

y[-1] = y[-2] = y[-3] = 0;

sol = 
        Solve[ZTransform[eq,  n,  z],  ZTransform[y[n], n, z]]  //  
          Flatten;

impulseTransferFunction  =  (ZTransform[y[n], n, z]  //.  
              sol~Join~{y[-3] -> 0, y[-2] -> 0,  
                  y[-1] -> 0})  //  
          Simplify; 
impulseTransferFunction  //  TraditionalForm

Compute poles of impulse response transfer function

poles = Solve[Denominator[impulseTransferFunction]  ==  0,  
            z]  //  Flatten  //  N ; 
    poles //  TraditionalForm

Plot pole locations

Show[
        Graphics[{{Red, Circle[{0, 0}, 1]}, {Orange, 
              PointSize[ .04], (Point[{Re[z], 
                          Im[z]}]  /.  #) &  /@  poles}}],  
        AspectRatio -> Automatic, Frame -> True,  
        Axes -> True, 
        PlotLabel -> 
          StyleForm["Complex plane",  FontFamily -> Times,  
            FontSize -> 18]];

Analytic impulse response

impulseResponse = 
        InverseZTransform[impulseTransferFunction,  z,  n];
DiscreteListPlot[Table[{n, impulseResponse},  {n, 0, 100}], Red, 
        PlotRange -> All,  Frame -> True,  
        Axes -> {True, False}];
impulseResponse  //  TraditionalForm

Stable recursive equation

Clear memory

Remove[Global`y];
Remove[Global`x];
Remove[Global`n];

Define difference equation and initial conditions

y[
          n_] := y[n] = -1/3  y[n - 1] + 1/2  y[n - 2]  +  
            DiscreteDelta[n];

y[-1] = y[-2] = y[-3] = 0;

Simulated output for impulse input

Input

input = 
        Table[{n, DiscreteDelta[n]}, {n, 0, 
            100}];
DiscreteListPlot[input, Blue, PlotRange -> {- .1, 2.1},  
        Frame -> True,  Axes -> {True, False}];

Output

output = Table[{n, y[n]}, {n, 0, 100}];
DiscreteListPlot[output, Red, PlotRange -> All,  
        Frame -> True,  Axes -> {True, False}];

Theoretical solution

Compute z-transform and transfer function

Remove[Global`y];
Remove[Global`x];
Remove[Global`n];

eq = (y[n] == -1/3  y[n - 1] + 1/2  y[n - 2] + 
              DiscreteDelta[n]);

y[-1] = y[-2] = y[-3] = 0;

sol = 
        Solve[ZTransform[eq,  n,  z],  ZTransform[y[n], n, z]]  //  
          Flatten;

impulseTransferFunction  =  (ZTransform[y[n], n, z]  //.  
              sol~Join~{y[-3] -> 0, y[-2] -> 0,  
                  y[-1] -> 0})  //  
          Simplify; 
impulseTransferFunction  //  TraditionalForm

Compute poles of impulse response transfer function

poles = Solve[Denominator[impulseTransferFunction]  ==  0,  
            z]  //  Flatten  //  N ; 
    poles //  TraditionalForm

Plot pole locations

Show[
        Graphics[{{Red, Circle[{0, 0}, 1]}, {Orange, 
              PointSize[ .04], (Point[{Re[z], 
                          Im[z]}]  /.  #) &  /@  poles}}],  
        AspectRatio -> Automatic, Frame -> True,  
        Axes -> True, 
        PlotLabel -> 
          StyleForm["Complex plane",  FontFamily -> Times,  
            FontSize -> 18]];

Analytic impulse response

impulseResponse = 
        InverseZTransform[impulseTransferFunction,  z,  n];
DiscreteListPlot[Table[{n, impulseResponse},  {n, 0, 100}], Red, 
        PlotRange -> All,  Frame -> True,  
        Axes -> {True, False}];
impulseResponse  //  TraditionalForm

Nonlinear, recursive difference equation

Nonlinear difference equations can behave very unpredictably, and are much harder to analyze than linear difference equation. Below we illustrate this with a simple example of the following difference equation, where R is some real number greater than 0.

y[n]  =  R y[n - 1] (1 - y[n - 1])

Experiment #1

Difference equation definition

Remove[Global`y];
y[n_] := y[n] = R y[n - 1] (1 - y[n - 1]);

Initial condition and constant R

y[0] = 0.5;
R = 1. ;

Simulation result

output = Table[y[n], {n, 0, 100}];
g1 = 
        ListPlot[output,  PlotRange -> All,  Axes -> False,  
          PlotStyle -> {Red,  Thickness[ .01]}, 
          PlotLabel -> 
            StyleForm["R=" <> ToString[R],  FontFamily -> Times,  
              FontSize -> 14,  FontWeight -> Bold]];

Experiment #2

Difference equation definition

Remove[Global`y];
y[n_] := y[n] = R y[n - 1] (1 - y[n - 1]);

Initial condition and constant R

y[0] = 0.5;
R = 1.1;

Simulation result

output = Table[y[n], {n, 0, 100}];
g2 = 
        ListPlot[output,  PlotRange -> All,  Axes -> False,  
          PlotStyle -> {Magenta,  Thickness[ .01]}, 
          PlotLabel -> 
            StyleForm["R=" <> ToString[R],  FontFamily -> Times,  
              FontSize -> 14,  FontWeight -> Bold]];

Experiment #3

Difference equation definition

Remove[Global`y];
y[n_] := y[n] = R y[n - 1] (1 - y[n - 1]);

Initial condition and constant R

y[0] = 0.5;
R = 2. ;

Simulation result

output = Table[y[n], {n, 0, 100}];
g3 = 
        ListPlot[output,  PlotRange -> All,  Axes -> False,  
          PlotStyle -> {Green,  Thickness[ .01]}, 
          PlotLabel -> 
            StyleForm["R=" <> ToString[R],  FontFamily -> Times,  
              FontSize -> 14,  FontWeight -> Bold]];

Experiment #4

Difference equation definition

Remove[Global`y];
y[n_] := y[n] = R y[n - 1] (1 - y[n - 1]);

Initial condition and constant R

y[0] = 0.5;
R = 3. ;

Simulation result

output = Table[y[n], {n, 0, 100}];
g4 = 
        ListPlot[output,  PlotRange -> All,  Axes -> False,  
          PlotStyle -> {Blue,  Thickness[ .01]}, 
          PlotLabel -> 
            StyleForm["R=" <> ToString[R],  FontFamily -> Times,  
              FontSize -> 14,  FontWeight -> Bold]];

Experiment #5

Difference equation definition

Remove[Global`y];
y[n_] := y[n] = R y[n - 1] (1 - y[n - 1]);

Initial condition and constant R

y[0] = 0.5;
R = 3.1;

Simulation result

output = Table[y[n], {n, 0, 100}];
g5 = 
        ListPlot[output,  PlotRange -> All,  Axes -> False,  
          PlotStyle -> {Cyan,  Thickness[ .01]}, 
          PlotLabel -> 
            StyleForm["R=" <> ToString[R],  FontFamily -> Times,  
              FontSize -> 14,  FontWeight -> Bold]];

Experiment #6

Difference equation definition

Remove[Global`y];
y[n_] := y[n] = R y[n - 1] (1 - y[n - 1]);

Initial condition and constant R

y[0] = 0.5;
R = 3.7;

Simulation result

output = Table[y[n], {n, 0, 100}];
g6 = 
        ListPlot[output,  PlotRange -> All,  Axes -> False,  
          PlotStyle -> {Orange,  Thickness[ .01]}, 
          PlotLabel -> 
            StyleForm["R=" <> ToString[R],  FontFamily -> Times,  
              FontSize -> 14,  FontWeight -> Bold]];

All results

y[n]  =  R y[n - 1] (1 - y[n - 1])

Show[GraphicsArray[{{g1, g2}, {g3, g4}, {g5, g6}}],  
        ImageSize -> 640];