Frequency response of FIR filters

Mathematica definitions

Load needed packages for this notebook

<< Graphics`Colors`;
<< Calculus`FourierTransform`;
<< 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}];
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}]
];

(1) Simplest low-pass FIR filter example

y[n]  =  1/2 (x[n]  +  x[n - 1])

Define impulse response sequence h(n)

L = 2;
hn  =  Table[{i, 1/L},  {i, 0, L - 1}]

Frequency response function H((e^j)θ)

Frequency response function H((e^j)θ)

Hθ  = DTFT[hn, θ]  //  Chop;  
    Hθ  //  TraditionalForm

Plots of frequency response function H((e^j)θ)as a function of θ

g1 = 
        Plot[Abs[Hθ], {θ, -π, π},  
          PlotStyle -> {Red},  PlotRange -> All, 
          PlotPoints -> 200,  
          DisplayFunction -> Identity];
g2 = 
        Plot[Arg[Hθ], {θ, -π, π},  
          PlotStyle -> {Red},  PlotRange -> All, 
          PlotPoints -> 200,  
          DisplayFunction -> Identity];
Show[GraphicsArray[{{g1, g2}}], ImageSize -> 504];

Output of system for simple cosine input: time-domain convolution

Sampling frequency

fs = 500;

Part of input sequence

f1[t_] := Cos[2  π 50  t] ;
xn  =  Table[f1[t], {t, 0, 1 - 1/fs, 1/fs}]  //  N;
g1 = 
        DiscreteListPlot[Take[xn, {100, 150}],  Blue,  
          DisplayFunction -> Identity,  
          Frame -> True];

y[n] =  x[n] * h[n]

yn = 
        Drop[ListConvolve[Last  /@  hn, xn, {1, -1}, 
            0], -1];

g2 = 
        DiscreteListPlot[Take[yn, {100, 150}],  Red,  
          DisplayFunction -> Identity,  Frame -> True,  
          Evaluate[AbsoluteOptions[g1, 
                PlotRange][[1]]]];
Show[GraphicsArray[{{g1, g2}}],  
        ImageSize -> 504];

Output of system for simple cosine input: use frequency-response function

Value H((e^j)θ)at 50Hz (θ = 2π 50/(Subscript[f, s])= 0.2π)

H50  =  
      Hθ  /.  θ -> (2  π 50/fs)

{Abs[H50],  Arg[H50]}

\(TraditionalForm\`y[
        n] = \(\(|\)\(H(
          e\^jθ)\)\(|\)\(cos[\(\(0.2  π\ n\)\(+\)\)\  ↼ H(e\^jθ)]\)\)\)

yn = 
        Table[{n - 99, Abs[H50] Cos[2  π n 50/fs  +  Arg[H50]]}, {n, 
            0, 499}];

g3 = 
        DiscreteListPlot[Take[yn, {100, 150}],  Green,  
          DisplayFunction -> Identity,  Frame -> True,  
          Evaluate[AbsoluteOptions[g1, 
                PlotRange][[1]]]];
Show[GraphicsArray[{{g1, g2}, {g1, g3}}],  
        ImageSize -> 504];

(2) Example of phase distortion

Linear phase (no distortion)

x1[n_] := 
        Cos[2  π n/50]  +  
          2 Cos[4  π n/50];
y1[n_] := 
        Cos[2  π n/50  -  5 (2  π/50)]  +  
          2 Cos[4  π n/50  -  5 (4  π/50)];

xn  =  Table[{n, x1[n]}, {n, 0, 100}];
xn2  =  
        Table[{n, y1[n]}, {n, 0, 100}];
g1a = 
        DiscreteListPlot[xn,  Blue,  DisplayFunction -> Identity,  
          Frame -> True];
g2a = 
        DiscreteListPlot[xn2,  Red,  DisplayFunction -> Identity,  
          Frame -> True];

ga = 
        Show[GraphicsArray[{{g1a, g2a}}],  
          ImageSize -> 504];

Nonlinear phase (distortion)

x1[n_] := 
        Cos[2  π n/50]  +  
          2 Cos[4  π n/50];
y1[n_] := 
        Cos[2  π/50  n  - 5]  +  
          2 Cos[4  π /50  n  -  5];

xn  =  Table[{n, x1[n]}, {n, 0, 100}];
xn2  =  
        Table[{n, y1[n]}, {n, 0, 100}];
g1b = 
        DiscreteListPlot[xn,  Blue,  DisplayFunction -> Identity,  
          Frame -> True];
g2b = 
        DiscreteListPlot[xn2,  Red,  DisplayFunction -> Identity,  
          Frame -> True];

gb = 
        Show[GraphicsArray[{{g1b, g2b}}],  
          ImageSize -> 504];

(3) Time-domain convolution ⟺ frequency-domain multiplication: an example

Random x[n] and h[n]

xn  =  
        Table[{i,  Random[Integer, {-10, 10}]}, {i, 0, 
            9}];
hn  =  
        Table[{i, Random[Integer,  {-10, 10}]}, {i, 0, 
            9}];

Plot x[n] and h[n]

Input sequence x[n]

DiscreteListPlot[xn,  Blue,  Frame -> True, 
        PlotRange -> All,  Axes -> {True, False}];

Impulse response function h[n]

DiscreteListPlot[hn, Blue,  Frame -> True, 
        PlotRange -> All,  Axes -> {True, False}];

Frequency-domain solution

Y(
        (e^j)θ) = X((e^j)θ) H(
          (e^j)θ)

Xθ  =  DTFT[xn, θ];
Hθ  =  DTFT[hn, θ];
Yθ  =  Xθ*Hθ;

Plot frequency representation

g1 = 
        Plot[Abs[Xθ*Hθ], {θ, -π, π},  
          PlotStyle -> {Red},  PlotRange -> All, 
          PlotPoints -> 200,  
          DisplayFunction -> Identity];
g2 = 
        Plot[Arg[Xθ* Hθ], {θ, -π, π},  
          PlotStyle -> {Red},  PlotRange -> All, 
          PlotPoints -> 200,  
          DisplayFunction -> Identity];
Show[GraphicsArray[{{g1, g2}}], ImageSize -> 504];

Compute inverse DTFT y[
          n] = (∫_(-π)^π) Y(
              (e^j)θ) (e^jn)θ dθ

yn =  
        1/(2  π) (∫_(-π)^π) 
              Yθ Exp[
                I n θ] \[DifferentialD]θ;

Plot resulting time-domain sequence y[n]

yn1  =  
      Table[{i, Limit[yn, n -> i]},  {i, 0,  
          Length[xn] + Length[hn] - 2}]; 
    g3  =  DiscreteListPlot[yn1,  Red,  Frame -> True, 
        PlotRange -> All,  Axes -> {True, False}];

Time-domain solution

y[n] = x[n] * h[n] (convolution)

yn  = 
        ListConvolve[Last  /@  hn, Last  /@  xn, {1, -1}, 
          0];

Plot resulting time-domain sequence y[n]

yn2  =  {Table[n,  {n,  0,  Length[yn] - 1}],  yn}  //  
        Transpose; 
    g6  =  DiscreteListPlot[yn2,  Green,  Frame -> True, 
        PlotRange -> All,  Axes -> {True, False}];

Compute Y((e^j)θ) from  y[n]

Yθ  =  DTFT[yn2, θ];

Plot frequency representation

g4 = 
        Plot[Abs[Yθ], {θ, -π, π},  
          PlotStyle -> {Green},  PlotRange -> All, 
          PlotPoints -> 200,  
          DisplayFunction -> Identity];
g5 = 
        Plot[Arg[Yθ], {θ, -π, π},  
          PlotStyle -> {Green},  PlotRange -> All, 
          PlotPoints -> 200,  
          DisplayFunction -> Identity];
Show[GraphicsArray[{{g4, g5}}], ImageSize -> 504];

Show both solutions

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