Introduction to FIR filter design: Part II

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

Plot discrete-time signal as sequence of amplitude lines

DiscreteListPlot3[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[ .01], 
                  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) Music signal

Import wav file

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

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

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

fs = 32000;

Listen to wav file

ListPlay[wav1, SampleRate -> fs];

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

(2) General filter

General filter definition

m1  =  (a - 1)/Δ;
b1  =  (a -  θc1 m1);

    
m2  =  (a - 
              1)/(-Δ);
b2  =  (a  -  θc2 m2);

    
Hθp  =  (u[θ] -  
              u[θ - (θc1  -  Δ)])  +  
(m1 θ  +  
                b1) (u[θ  - (θc1  -  Δ)]  -  
                u[θ - θc1])  + 
a (u[θ - θc1]  -  
                u[θ  -  θc2])  + 
(m2 θ  +  b2) (u[θ  - θc2]  -  
                u[θ - (θc2  +  Δ)])  + 
(u[θ - (θc2  +  Δ)]  -  u[θ - π]);

    
Hθn  =  
        Hθp  /.  θ -> (-θ);

    
Hθ  =  Hθp  +  Hθn;

Numeric values

rules = {a -> 2,  θc1 -> 1,  θc2 -> 
            3/2,  Δ -> 1/20};

Visualize filter

Plot[Hθ /.  rules,  {θ, -π, π},  
        PlotRange -> {- .01, 2.1},  
        PlotStyle -> Blue];

Symbolic solution

\(hns1\  = \ 
      1/\((2  π)\)\ \((
∫\_\(-π\)\%\(-\((θc2 + Δ)\)\)Exp[
                  I\ n\ θ] \[DifferentialD]θ + \ ∫\_\(-\((θc2 + Δ)\)\)\%\(-θc2\)\((\((a - 1)\)/Δ\ \((θ\  + \ θc2)\)\  + \ a)\) 
                Exp[I\ n\ θ] \[DifferentialD]θ + ∫\_\(-θc2\)\%\(-θc1\)a\ Exp[
                  I\ n\ θ] \[DifferentialD]θ + ∫\_\(-θc1\)\%\(-\((θc1 - Δ)\)\)\((\((1 - 
                            a)\)/Δ\ \((θ\  + \ θc1)\)\  + \ a)\) 
                Exp[I\ n\ θ] \[DifferentialD]θ + ∫\_\(-\((θc1 - Δ)\)\)\%\((θc1 - Δ)\)Exp[
                  I\ n\ θ] \[DifferentialD]θ\  + ∫\_\((θc1 - Δ)\)\%θc1\((\((a - 
                            1)\)/Δ\ \((θ\  - \ θc1)\)\  + \ a)\) 
                Exp[I\ n\ θ] \[DifferentialD]θ + ∫\_θc1\%θc2 
                a\ Exp[I\ n\ θ] \[DifferentialD]θ + ∫\_θc2\%\((θc2 + Δ)\)\((\((1 - 
                            a)\)/Δ\ \((θ\  - \ θc2)\)\  + \ a)\) 
                Exp[I\ n\ θ] \[DifferentialD]θ + ∫\_\((θc2 + Δ)\)\%π Exp[
                  I\ n\ θ] \[DifferentialD]θ)\)\)

hns  =  hns1  //  FullSimplify;
hns  //  TraditionalForm

Limit[Sin[n π]/(n π),  n -> 0]

Corresponding approximate impulse response function (with delay)

nmax  =  50;

hn = 
        Table[{i + 
              nmax, If[i  ==  0,  
                  Limit[hns  /.  rules,  
                    n -> 0],  hns  /.  rules  /.  
                    n -> i]  //  N  //  Chop}, {i, -nmax, 
            nmax}];

ghn = 
        DiscreteListPlot2[hn,  Red,  PlotRange -> {- .21,  .21},  
          Frame -> True, Axes -> {True, False},  
          ImageSize -> 504,  AspectRatio -> 1/3];

Limit[hns  /.  rules,  n -> 0]

(3) Music FIR filter experiment #1

Composite filter

Numeric values

nmax = 50;
fs  =  32000;

First sub filter

fc1  =  375;
fc2  =  1500;
a1 =  10^(12/20);
Δ1  =  100;

rulesf1  =  {a -> 
            a1, θc1 -> (2  π fc1/
                  fs),  θc2  ->  (2  π fc2/
                  fs), Δ -> (2  π Δ1/fs)};

Hθ1  =  Hθ  /.  rulesf1;

Second sub filter

fc1  =  6000;
fc2  =  15000;
a2 =  (10)^(-12/20);
Δ2  =  100;

rulesf2  =  {a -> 
            a2,  θc1 -> (2  π fc1/
                  fs),  θc2  ->  (2  π fc2/
                  fs), Δ -> (2  π Δ2/fs)};

Hθ2  =  Hθ  /.  rulesf2;

Composite filter

Hθall  =  Hθ1*Hθ2;

gf = 
        Plot[Hθall,  {θ,  -π,  π},  
          PlotStyle -> Blue,  PlotRange -> {- .1, 4.1},  
          PlotPoints -> 200];

gfdb = 
        Plot[20  Log[10, 
              Hθall],  {θ,  -π,  π},  
          PlotStyle -> Blue,  PlotRange -> All,  
          PlotPoints -> 200];

Show[GraphicsArray[{gf, gfdb}],  
        ImageSize -> 504];

Corresponding approximate impulse response function (with delay)

First filter

nmax  =  500;
hn1 = 
        Table[If[i  ==  0,  
                Limit[hns  /.  rulesf1,  
                  n -> 0],  hns  /.  rulesf1  /.  
                  n -> i]  //  N  //  Chop, {i, -nmax, 
            nmax}];
ghn1 = 
        DiscreteListPlot2[hn1,  Red,  
          PlotRange -> {- .21,  .31},  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 504,  
          AspectRatio -> 1/3];

Second filter

nmax  =  500;
hn2 = 
        Table[If[i  ==  0,  
                Limit[hns  /.  rulesf2,  
                  n -> 0],  hns  /.  rulesf2  /.  
                  n -> i]  //  N  //  Chop, {i, -nmax, 
            nmax}];
ghn2 = 
        DiscreteListPlot2[hn2,  Red,  
          PlotRange -> {- .21,  .31},  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 504,  
          AspectRatio -> 1/3];

Show[GraphicsArray[{{ghn1},  {ghn2}}],  
        ImageSize -> 504];

Composite filter

hn = 
        ListConvolve[hn1, hn2, {1, -1}, 0];
ghn = 
        DiscreteListPlot2[hn,  Red,  
          PlotRange -> {- .055,  .055},  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 504,  
          AspectRatio -> 1/3];

hn2  =  Take[hn,  {1001 - 300, 1001 + 300}];

ghn2 = 
        DiscreteListPlot2[hn2,  Red,  
          PlotRange -> {- .055,  .055},  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 504,  
          AspectRatio -> 1/3];

Corresponding frequency response

As a function of θ

hnb  = {Table[n, {n, 0,  Length[hn2] - 1}], hn2}  //  
          Transpose;

Hθb  =  DTFT[hnb,  θ];
gfr1 = 
        Plot[Abs[Hθb], {θ, -π, π},  
          PlotStyle -> {Red},  PlotRange -> {- .1, 4.1}, 
          PlotPoints -> 200];

As a function of f

gfr2 = 
        Plot[Abs[
              Hθb]  /.  θ -> (2  π f/
                    fs), {f, -fs/2, fs/2},  PlotStyle -> {Red},  
          PlotRange -> All, PlotPoints -> 200,  
          DisplayFunction -> Identity];
gfr2b  =  
        Show[gfr2,  
          FrameTicks -> {{{-15000, "-15"},  {-10000, "-10\
"},  {-5000,  "-5"}, {0, "0"}, {5000, "5"}, {10000, "10\
"}, {15000, "15"}}, Automatic, None, Automatic},  
          PlotRange -> {- .1, 4.1},  
          DisplayFunction -> $DisplayFunction];

Show[GraphicsArray[{gfr1, gfr2b}], 
        ImageSize -> 504];

Convolve composite filter with music signal

Perform convolution

wav2 = ListConvolve[hn, wav1, {1, -1}, 0];

Listen to wav file

ListPlay[wav2, SampleRate -> fs];

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

Look at frequency content (spectrogram) of two wav files

Unfiltered sample

ns = fs;

spectrogram1  =  
        Transpose[Take[-Abs[Fourier[#]], {1, ns/2}] &  /@  
            Partition[wav1, ns, ns/2]];

gs1 = 
        ListDensityPlot[spectrogram1, Mesh -> False,  
          PlotRange -> Automatic,  
          MeshRange -> {{0, 23}, {0, 16000}},  
          AspectRatio -> 1/GoldenRatio,  
          ImageSize -> 504];

Filtered sample

ns = fs;

spectrogram2  =  
        Transpose[Take[-Abs[Fourier[#]], {1, ns/2}] &  /@  
            Partition[wav2, ns, ns/2]];

gs2 = 
        ListDensityPlot[spectrogram2, Mesh -> False,  
          PlotRange -> Automatic,  
          MeshRange -> {{0, 23}, {0, 16000}},  
          AspectRatio -> 1/GoldenRatio,  
          ImageSize -> 504];

Show both

Show[GraphicsArray[{{gs1}, {gs2}}], 
        ImageSize -> 504];

Show one column of each spectrogram

\(\(gexample1 = 
        DiscreteListPlot2[\(-\(Transpose[spectrogram1]\)[\([20]\)]\), Red, \ 
          Frame -> True, PlotRange -> {\(-10\), 10^7}, \ 
          FrameTicks -> {{{0, "0"}, {5000, "5"}, {10000, "10"}, {15000, "15"}}, Automatic, None, Automatic}, \ 
          ImageSize -> 504, AspectRatio -> 1/3];\)\)

\(\(gexample2 = 
        DiscreteListPlot2[\(-\(Transpose[spectrogram2]\)[\([20]\)]\), Red, \ 
          Frame -> True, PlotRange -> {\(-10\), 10^7}, 
          FrameTicks -> {{{0, "0"}, {5000, "5"}, {10000, "10"}, {15000, "15"}}, Automatic, None, Automatic}, 
          ImageSize -> 504, AspectRatio -> 1/3];\)\)

Show[GraphicsArray[{{gexample1}, {gexample2}}], 
        ImageSize -> 504];

(4) Music FIR filter experiment #2

Composite filter

Numeric values

nmax = 50;
fs  =  32000;

First sub filter

fc1  =  375;
fc2  =  1500;
a1 =  10^(-12/20);
Δ1  =  100;

rulesf1  =  {a -> 
            a1, θc1 -> (2  π fc1/
                  fs),  θc2  ->  (2  π fc2/
                  fs), Δ -> (2  π Δ1/fs)};

Hθ1  =  Hθ  /.  rulesf1;

Second sub filter

fc1  =  6000;
fc2  =  15000;
a2 =  (10)^(12/20);
Δ2  =  100;

rulesf2  =  {a -> 
            a2,  θc1 -> (2  π fc1/
                  fs),  θc2  ->  (2  π fc2/
                  fs), Δ -> (2  π Δ2/fs)};

Hθ2  =  Hθ  /.  rulesf2;

Composite filter

Hθall  =  Hθ1*Hθ2;

gf = 
        Plot[Hθall,  {θ,  -π,  π},  
          PlotStyle -> Blue,  PlotRange -> {- .1, 4.1},  
          PlotPoints -> 200];

gfdb = 
        Plot[20  Log[10, 
              Hθall],  {θ,  -π,  π},  
          PlotStyle -> Blue,  PlotRange -> All,  
          PlotPoints -> 200];

Show[GraphicsArray[{gf, gfdb}],  
        ImageSize -> 504];

Corresponding approximate impulse response function (with delay)

First filter

nmax  =  500;
hn1 = 
        Table[If[i  ==  0,  
                Limit[hns  /.  rulesf1,  
                  n -> 0],  hns  /.  rulesf1  /.  
                  n -> i]  //  N  //  Chop, {i, -nmax, 
            nmax}];
ghn1 = 
        DiscreteListPlot2[hn1,  Red,  
          PlotRange -> {- .21,  .31},  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 504,  
          AspectRatio -> 1/3];

Second filter

nmax  =  500;
hn2 = 
        Table[If[i  ==  0,  
                Limit[hns  /.  rulesf2,  
                  n -> 0],  hns  /.  rulesf2  /.  
                  n -> i]  //  N  //  Chop, {i, -nmax, 
            nmax}];
ghn2 = 
        DiscreteListPlot2[hn2,  Red,  
          PlotRange -> {- .21,  .31},  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 504,  
          AspectRatio -> 1/3];

Show[GraphicsArray[{{ghn1},  {ghn2}}],  
        ImageSize -> 504];

Composite filter

hn = 
        ListConvolve[hn1, hn2, {1, -1}, 0];
ghn = 
        DiscreteListPlot2[hn,  Red,  
          PlotRange -> {- .055,  .055},  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 504,  
          AspectRatio -> 1/3];

hn2  =  Take[hn,  {1001 - 300, 1001 + 300}];

ghn2 = 
        DiscreteListPlot2[hn2,  Red,  
          PlotRange -> {- .055,  .055},  Frame -> True, 
          Axes -> {True, False},  ImageSize -> 504,  
          AspectRatio -> 1/3];

Corresponding frequency response

As a function of θ

hnb  = {Table[n, {n, 0,  Length[hn2] - 1}], hn2}  //  
          Transpose;

Hθb  =  DTFT[hnb,  θ];
gfr1 = 
        Plot[Abs[Hθb], {θ, -π, π},  
          PlotStyle -> {Red},  PlotRange -> {- .1, 4.1}, 
          PlotPoints -> 200];

As a function of f

gfr2 = 
        Plot[Abs[
              Hθb]  /.  θ -> (2  π f/
                    fs), {f, -fs/2, fs/2},  
          PlotStyle -> {Red},   PlotRange -> {- .1, 4.1}, 
          PlotPoints -> 200, 
          FrameTicks -> {{{-15000, "-15"},  {-10000, "-10\
"},  {-5000,  "-5"}, {0, "0"}, {5000, "5"}, {10000, "10\
"}, {15000, "15"}}, Automatic, None, Automatic}];

Show[GraphicsArray[{gfr1, gfr2}], 
        ImageSize -> 504];

Convolve composite filter with music signal

Perform convolution

wav2 = ListConvolve[hn, wav1, {1, -1}, 0];

Listen to wav file

ListPlay[wav2, SampleRate -> fs];

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

Look at frequency content (spectrogram) of two wav files

Unfiltered sample

ns = fs;

spectrogram1  =  
        Transpose[Take[-Abs[Fourier[#]], {1, ns/2}] &  /@  
            Partition[wav1, ns, ns/2]];

gs1 = 
        ListDensityPlot[spectrogram1, Mesh -> False,  
          PlotRange -> Automatic,  
          MeshRange -> {{0, 23}, {0, 16000}},  
          AspectRatio -> 1/GoldenRatio,  
          ImageSize -> 504];

Filtered sample

ns = fs;

spectrogram2  =  
        Transpose[Take[-Abs[Fourier[#]], {1, ns/2}] &  /@  
            Partition[wav2, ns, ns/2]];

gs2 = 
        ListDensityPlot[spectrogram2, Mesh -> False,  
          PlotRange -> Automatic,  
          MeshRange -> {{0, 23}, {0, 16000}},  
          AspectRatio -> 1/GoldenRatio,  
          ImageSize -> 504];

Show both

Show[GraphicsArray[{{gs1}, {gs2}}], 
        ImageSize -> 504];

Show one column of each spectrogram

\(\(gexample1 = 
        DiscreteListPlot2[\(-\(Transpose[spectrogram1]\)[\([20]\)]\), Red, \ 
          Frame -> True, PlotRange -> {\(-10\), 10^7}, \ 
          FrameTicks -> {{{0, "0"}, {5000, "5"}, {10000, "10"}, {15000, "15"}}, Automatic, None, Automatic}, \ 
          ImageSize -> 504, AspectRatio -> 1/3];\)\)

\(\(gexample2 = 
        DiscreteListPlot2[\(-\(Transpose[spectrogram2]\)[\([20]\)]\), Red, \ 
          Frame -> True, PlotRange -> {\(-10\), 10^7}, 
          FrameTicks -> {{{0, "0"}, {5000, "5"}, {10000, "10"}, {15000, "15"}}, Automatic, None, Automatic}, 
          ImageSize -> 504, AspectRatio -> 1/3];\)\)

Show[GraphicsArray[{{gexample1}, {gexample2}}], 
        ImageSize -> 504];