Introduction to signals and systems

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

Sample signals

Simple signals

Example #1

s1  =  Sin[2  π t];

Plot[s1, {t, 0, 2},  
        PlotStyle -> Blue];
ƒ[t]  == s1  //  TraditionalForm

Example #2

s2  = Sin[2  π 500 t];

Plot[s2, {t, 0, 2},  PlotPoints -> 100,  
        PlotStyle -> Blue]; 
ƒ[t] == s2  //  TraditionalForm

Play[s2, {t, 0, 2}];

Example #3

s3 = Sin[π t] Sin[2  π 500 t];

Plot[s3, {t, 0, 2},  
        PlotStyle -> Blue];
ƒ[t] == s3  //  TraditionalForm

Play[s3, {t, 0, 2}];

Speech signal

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;

Listen to wav file

ListPlay[wav1, SampleRate -> sampleRate1];

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

Visualize wav file

Normalize wav sound file.

normalized1  =  wav1/(2. ^15);
sound1 = FixTime[normalized1, sampleRate1];

Plot entire wav sound file.

ListPlot[sound1, PlotRange -> All,  
        PlotStyle -> Blue];

Zoom in a little

ListPlot[Take[sound1, {11000, 12000}],  PlotRange -> All,  
        PlotStyle -> Blue];

Zooming in some more

g1  = 
        ListPlot[Take[sound1, {11000, 11100}],  PlotRange -> All,  
          PlotStyle -> Blue];

Plot above waveform as a discrete-time sequence

g2  = 
        DiscreteListPlot[Take[sound1, {11000, 11100}],  Blue, 
          PlotRange -> All,  Frame -> True,  
          Axes -> {True, False}];

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];
wav2 = ReadSoundFile[file, PrintHeader -> True];

sampleRate2 = 32000;

Listen to wav file

ListPlay[wav2, 
        SampleRate -> sampleRate2];

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

Sample images

Gray-scale image

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

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

SetDirectory[dir];
im1 = Import[file]; 
Show[im1];

Color image

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

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

SetDirectory[dir];
im2 = Import[file]; 
Show[im2];

Continuous-time vs. discrete-time signals

Continuous-time signal #1

f1[t_] := 
        Exp[-(t - 4)^2] 
          Sin[3  π t];
ƒ[t] == f1[t]  //  TraditionalForm

g1 = 
        Plot[f1[t], {t, 0, 8},  PlotRange -> All,  
          PlotStyle -> {Blue}];

Sample discrete-time signal (continuous-time to discrete-time conversion)

Adjust sampling rate below (in Hz) to see different sampling results.

sampleRate = 10;

sampled  = 
        Table[{t, f1[t]}, {t, 0, 8, 
            1/sampleRate}];
g2 = 
        DiscreteListPlot[sampled, Red, PlotRange -> All,  
          Frame -> True,  
          Axes -> {True, False}];
Show[g2, g1];

Sample discrete-time signal (continuous-time to discrete-time conversion, very undersampled)

Adjust sampling rate below (in Hz) to see different sampling results.

sampleRate = 2;

sampled  = 
        Table[{t, f1[t]}, {t, 0, 8, 
            1/sampleRate}];
g2 = 
        DiscreteListPlot[sampled, Red, PlotRange -> All,  
          Frame -> True,  
          Axes -> {True, False}];
Show[g2, g1];

Continuous-time signal #2

Pulse1[t_] := 
        UnitStep[t]  -  UnitStep[t - 1/2]  -  
          1/2;
fnc1 = PeriodicSignal[Pulse1, {0, 1}, {t, 0, 4}];

g1 = 
        Plot[fnc1, {t, 0, 4},  PlotStyle -> {Blue}, 
          PlotRange -> All];

Sample signal (continuous-time to discrete-time conversion)

Adjust sampling rate below (in Hz) to see different sampling results.

sampleRate = 6;

sampled  = 
        Table[{ti, fnc1  /.  t -> ti}, {ti, 0, 4, 
            1/sampleRate}];
g2 = 
        DiscreteListPlot[sampled, Red, PlotRange -> All,  
          Frame -> True,  
          Axes -> {True, False}];
Show[g2, g1];

Listening to frequencies

Simple sine waves

400Hz sine wave

Play[Sin[2  π 400 t], {t, 0, 2}];

800Hz sine wave

Play[Sin[2  π 800 t], {t, 0, 2}];

2KHz sine wave

Play[Sin[2  π 2000 t], {t, 0, 2}];

More complex sine waves

Sum of two sine waves

f1[t_] := 
        Sin[2  π 400 t]  +  
          Sin[2  π 600 t];
Plot[f1[t], {t, 0, 0.02}];
Play[f1[t], {t, 0, 2}];
ƒ[t] == f1[t]  //  TraditionalForm

Product of two sine waves

f2[t_] := 
        Sin[π t] Sin[2  π 400 t];
Play[f2[t], {t, 0, 2}];
ƒ[t] == f2[t]  //  TraditionalForm

Time domain ⟺ Frequency domain

Example #1(a): simple sinusoid

Continuous-time waveform

f1[t_] := Sin[2  π 10  t];
Plot[f1[t], {t, 0, 1},  
        PlotStyle -> Blue];
ƒ[t] == f1[t]  //  TraditionalForm

Sampled continuous-time waveform

sampleRate = 1000;

signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
ListPlot[signal,  PlotStyle -> Blue];

Discrete spectrum: Fast Fourier Transform (FFT)

fft  =  Fourier[signal];
g1 = 
        DiscreteListPlot[FourierRearrange[Abs[fft], sampleRate], Red, 
          PlotRange -> All,  Frame -> True];

Show[g1,  PlotRange -> {{-20, 20}, All}];

Example #1(b): simple sinusoid with DC offset

Continuous-time waveform

f1[t_] := Sin[2  π 10  t] + 1;
Plot[f1[t], {t, 0, 1},  
        PlotStyle -> Blue];
ƒ[t] == f1[t]  //  TraditionalForm

Sampled continuous-time waveform

sampleRate = 1000;

signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
ListPlot[signal,  PlotStyle -> Blue];

Discrete spectrum: Fast Fourier Transform (FFT)

fft  =  Fourier[signal];
g1 = 
        DiscreteListPlot[FourierRearrange[Abs[fft], sampleRate], Red, 
          PlotRange -> All,  Frame -> True];

Show[g1,  PlotRange -> {{-20, 20}, All}];

Example #2: sum of two sinusoids

Continuous-time waveform

f1[t_] := 
        Sin[2  π 10  t]  +  
          2  Sin[2  π 30 t];
Plot[f1[t], {t, 0, 1},  
        PlotStyle -> Blue];
ƒ[t] == f1[t]  //  TraditionalForm

Sampled continuous-time waveform

sampleRate = 1000;

signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
ListPlot[signal,  PlotStyle -> Blue];

Discrete spectrum: Fast Fourier Transform (FFT)

fft  =  Fourier[signal];
g1 = 
        DiscreteListPlot[FourierRearrange[Abs[fft], sampleRate], Red, 
          PlotRange -> All,  Frame -> True];

Show[g1,  PlotRange -> {{-40, 40}, All}];

Example #3: product of two sinusoids

Continuous-time waveform

f1[t_] := 
        Sin[2  π t] Sin[2  π 30 t];
Plot[f1[t], {t, 0, 1},  
        PlotStyle -> Blue];
ƒ[t] == f1[t]  //  TraditionalForm

Sampled continuous-time waveform

sampleRate = 1000;

signal  =  
        Table[f1[t], {t, 0, 1, 1/(sampleRate - 1)}]  //  N;
ListPlot[signal,  PlotStyle -> Blue];

Discrete spectrum: Fast Fourier Transform (FFT)

fft  =  Fourier[signal];
g1 = 
        DiscreteListPlot[FourierRearrange[Abs[fft], sampleRate], Red, 
          PlotRange -> All,  Frame -> True];

Show[g1,  PlotRange -> {{0, 40}, All}];

Does this make sense?

ƒ[t]  == (f1[t]  // TrigReduce)  //  
      TraditionalForm

Note that the expression 'sin(2π t) sin(60π t)' can also be expressed as '1/2 [(cos(58π t)-cos(62π t)]'. Thus two frequencies are present in the signal 29Hz and 31Hz.

Example #4: square wave

Continuous-time waveform

Pulse1[t_] := 
        UnitStep[t] - UnitStep[t - 1/20] - 1/2;
fnc1  =  
        PeriodicSignal[
          Pulse1, {0, 1/10}, {t, 0, 1}];
Plot[fnc1, {t, 0, 1}, PlotPoints -> 100,  
        PlotStyle -> Blue];

Sampled continuous-time waveform

sampleRate = 1000;

signal  =  
        Table[fnc1  /.  t -> ti, {ti, 0, 1, 
              1/(sampleRate - 1)}]  //  N;
ListPlot[signal, PlotStyle -> Blue];

Discrete spectrum: Fast Fourier Transform (FFT)

fft  =  Fourier[signal];
g1 = 
        DiscreteListPlot[FourierRearrange[Abs[fft], sampleRate], Red, 
          PlotRange -> All,  Frame -> True];

Example #5: single pulse

Continuous-time waveform

f1[t_] := 
        UnitStep[t - 9/20] - 
          UnitStep[t - 11/20];
Plot[f1[t], {t, 0, 1}, PlotPoints -> 100,  
        PlotStyle -> Blue];

Sampled continuous-time waveform

sampleRate = 1000;

signal  =  
        Table[f1[t] , {t, 0, 1, 1/(sampleRate - 1)}]  //  
          N;
ListPlot[signal,  PlotStyle -> Blue];

Discrete spectrum: Fast Fourier Transform (FFT)

fft  =  Fourier[signal];
g1 = 
        DiscreteListPlot[FourierRearrange[Abs[fft], sampleRate], Red, 
          PlotRange -> All,  Frame -> True];

Show[g1, PlotRange -> {{-50, 50}, All}];