Murphy’s Laws in the DSP World (Part 1)

First Six Laws Explained

A Pandora’s box of unexpected issues gets opened the moment you move from the real world of analog signals and enter the world of digital signal processing (DSP). In Part 1 of this new article series, the authors explain the first six “Murphy’s Laws of DSP” and provide you with methods and techniques to navigate around them.

To analyze these signals requires that we must use a sensor to convert the original, long, continuous-in-time signals into short bursts of digital numbers that a computer program can handle. However, a Pandora’s box of unexpected issues gets opened the moment we move from the real world of analog signals and enter the world of digital signal processing (DSP). “For how long and how frequently should we sample the incoming signal?” become important things to know. Unfortunately, “Twice as much as half as much” is often a valid response—meaning there is not a “correct” answer, just a choice of less than perfect methods.

Take for example the MRI data shown in Photo 1. Imaging starts with an object or person placed in a strong magnetic field. Applying a special series of electromagnetic pulses and alternating magnetic fields cause the object to emit what are called k-space signals, which can be sampled—as shown in false colors in Photo 1a. The emitted signals contain information about the position of object in the main magnetic field and its physical characteristics. Applying an inverse discrete Fourier transform (IDFT) to the k-space data generates a digital image of the MRI phantom as displayed in Photo 1b. Note the fine detail in the comb-like feature at the center bottom of the picture, and the smoothness of the block features at the top right and bottom left.

Waiting a while between exciting and capturing the emitted signals means that a reconstructed brain image will emphasize cancer cells whose k-space signal strength decays more slowly than normal cells. Spending the same amount of time capturing individual heart images leads to blurring as the heart muscles are moving faster than the MRI data can be collected.

Two solutions have been offered for this problem. In Photo 1c the measurement time is reduced by capturing only a center portion of the k-space data, 1/9th of the data. This leads to the IDFT reconstructed image Photo 1d. There are a number of distortions introduced. For example, the bottom comb-like feature now looks blurred. The block features are no longer smooth but are full of distortions. These distortions are called edge ringing and are produced because the amount of k-space data has been reduced, or truncated. Application of various DSP algorithms, and prior knowledge of what the image “should” look like, CAN ALLOW the recovery of close to the original image quality.

Sparse imaging changes the characteristics of the alternating magnetic fields in order to collect a random selection of 1/9th of the k-space data outside a completely captured small central portion of k-space data, Photo 1e.

The capture of a wider range of random k-space values has allowed an increase in image resolution, clear comb, smoother blocks surfaces in Photo 1f but with an increased image noise level. Clever use of DSP algorithms, and more image prior knowledge, CAN ALLOW, the reduction of the increased noise level without a total loss of the improved resolution.

Notice the capitalized words: CAN ALLOW. We have to first identify the algorithms, then show, using representative data, that they improve things before finally implementing the algorithms on actual hardware. The installation of mathematical scripting language tools such as MATLAB or Octave, (see end of the article) can be used to handle the DSP simulations needed to tackle the first two parts. Then the generated input and output results can be used to validate the correctness, and/or limitations, of the hardware real-time implementation.

The trouble is any work done in a digital world has some very interesting and non-obvious differences compared to what the developer thinks is “common sense”. Just when the developer thinks that things are under control, problems jump out of the woodwork and the simulations just don’t give the results that you would expect.

Even worse, sometimes the simulation evaluations give great results, but later work with real data shows the algorithms actually failed! In this series of articles, we look at how variations of Murphy’s Law lurk to trap the unwary DSP developer.

MURPHY’S DSP LAWS I, II AND III

You don’t need to go very far into digital signal processing before finding the impact of Murphy. With the script Investigate_DSP_ML_I_II_III.m, Listing 1, we generate a 2 Hz signal which we have captured over a time duration of 1 second. Listing 2 contains the MATLAB / Octave class information to generate the following signal and display it:

However, we must get past the first Murphy DSP Law in order to properly investigate the results for a cosinusoidal signal phase φ=0°, and a sinusoidal signal, φ=-90°.

DSP_ML-I: Things will get misleading when you don’t remember that DSP signals consist of a finite number of sampled values.

As discussed early, the first thing to tackle is “How fast should you sample?”—a problem explored in Listing 1. A quick web search on “DSP sampling rates” will lead to the Nyquist criterion: “You must sample at least twice as fast as the fastest frequency in the signal being examined”. While this is true, such advice leads to the signals shown in Figure 1a with a 2 Hz signal sampled at 4.02 Hz, slightly faster than the Nyquist sampling criterion, 2×2 Hz. The generated signal is very difficult to visualize, even though all the signal information is mathematically present. Allowing the MATLAB/Octave packages to connect the sampled points, small circles, with lines is also misleading as these circles represent the sampled values of a smoothly varying signal and not the sampled edges of a jagged triangular signal.

DSP_ML-II: Any lines drawn to connect the DSP sample values are merely a visual convenience and are probably going to lead you to misinterpret the signal’s true characteristics and behavior.

Figure 1b shows that better visualization occurs with signals sampled at 16 Hz, which is 4 times the Nyquist rate. The fact the black and red sample points are related to cosinusoidal and sinusoidal signals becomes more obvious. However, we must remember the second DSP Murphy’s Law, which continually appears in many different forms during DSP analysis: Those straight edges and pointy peaks are not real!

Listing 2 also shows that we have had to do a not-so-obvious something special to avoid falling foul of the third DSP Murphy’s Law, which we will soon see provides interesting surprises on many occasions that we need to avoid! The DSP_ML-III issue to remember here is that the sample times and data at time t=m∆T starting at time
t =0 will be stored in array element m+1 and not in element m because of the way the arrays timeSampled[] and timeData[] are stored by MATLAB / Octave compared to other languages such as C/C++.

DSP_ML-III: To ensure your results pan out, don’t calculate something one sample point out!

We must watch out for many variants of DSP_ML-III during DSP programming.

MURPHY’S DSP LAWS IV AND V

In Listing 3, we explore how to visualize the frequency characteristics of time domain signals. In this script, we investigate the frequency characteristics of the time domain signal:

a)
%  Place this script in file “C:\DSP_ML\Investigate_DSP_ML_I_II_III.m” 

function Investigate_DSP_ML_I_II_III ( )
   clear classes;   %  Force MATLAB / OCTAVE to see class changes
   amplitude = 1;  duration  = 1;   
   phase = 0; %  Using phase 0 degrees generates a cosinusoid signal 
   DisplaySimpleSignal(amplitude, duration, phase);
   phase = -90;  % Phase -90 degrees generates a sinusoid signal      
   DisplaySimpleSignal(amplitude, duration, phase);

__________________________________________________________________________

b)
%  Place this script in file “C:\DSP_ML\DisplaySimpleSignal.m” 
  
function DisplaySimpleSignal(amplitude, duration, phase)
  frequency = 2;   aboveNyquistFreq = 2.01 * frequency;  
  fastSampling = 4 * 2 * frequency;    
  nyquistSampled = TimeDomainSignal(amplitude, frequency, phase, duration, aboveNyquistFreq);
  DisplayTimeSignal(nyquistSampled, ...
           sprintf('NYQUIST SAMPLED 2 Hz PHASE %d', phase), '-o');
  fastSampledSignal = TimeDomainSignal(amplitude,frequency, phase, duration, fastSampling);
  DisplayTimeSignal(fastSampledSignal, ...
                  sprintf('FAST SAMPLED 2 Hz PHASE %d', phase), '-o');

% Place in file ”C:\DSP_ML\@TimeDomainSignal\TimeDomainSignal.m”

classdef TimeDomainSignal   
  properties
    timeData = [ ];          numPts = [ ];
    timeSampled = [ ];    samplingFrequency = [ ];  
  end
 methods                                 % Time domain signal constructor
    function obj = TimeDomainSignal(Amp, freq, phase, duration, samplingFreq)
      obj.samplingFrequency = samplingFreq;
      deltaT = 1  / samplingFreq;
               % First sample at time 0 stored starting at array element 1
      obj.timeSampled = 0 : deltaT : duration;  
% Change  obj.timeSampled when fixing DSP_ML IV (after Figure 2b)   
      % obj.timeSampled = 0 : deltaT : (duration - deltaT);
      obj.numPts = length(obj.timeSampled);
      obj.timeData = Amp *  ...
                   cos(2*pi * (freq * obj.timeSampled + phase / 360));        
    end
                                            % Time domain display method
    function DisplayTimeSignal(signal, figureName, lineStyle)
       figure('Name', figureName);     hold on;
       lineStyle = strcat(lineStyle, 'k');
       plot(signal.timeSampled, signal.timeData, lineStyle);
       ylabel('\fontsize{14}AMPLITUDE');   
       xlabel('\fontsize{14}TIME WHEN SAMPLED (s)'); 
       pause(3);         % Short pause between displays=
    end
                                            % Overload plus (+) operator
   function result = plus(signal1,signal2)  
      result = signal1; 
      result.timeData = result.timeData + signal2.timeData;
   end
  end  % methods
end

% Place in file “C:\DSP_ML\Investigate_DSP_ML_IV_V.m”

function Investigate_DSP_ML_IV_V( ) 
  clear classes;   %  Force MATLAB / OCTAVE to see class changes
  fastSampling = 32; duration = 1; phase0 = 0; phase90 = -90;
 
% Dual Signal -- DC + 4 Hz  -- Phase 90
  signalDC = TimeDomainSignal(1.0, 0, phase0,  ...
                                                  duration, fastSampling);
  signal4HzPhase0 = TimeDomainSignal(1.0, 4, phase0, duration, fastSampling);
  dualSignalPhase0 = signalDC + signal4HzPhase0;

  FrequencySpectrum_V1( dualSignalPhase0,  ...
                          'Frequency Spectrum DC + 4 Hz Phase 0', '-o');  
  FrequencySpectrum_V2( dualSignalPhase0, ...
                          'Frequency Spectrum DC + 4 Hz Phase 0', '-o');
  
%% Dual Signal -- DC + 4 Hz  -- Phase -90 -- Freq all real
  signal4HzPhase90 = TimeDomainSignal(1.0, 4, phase90, ...
                                                 duration, fastSampling);
  dualSignalPhase90 = signalDC + signal4HzPhase90;
  FrequencySpectrum_V2(dualSignalPhase90, ...
                          'Frequency Spectrum DC + 4 Hz Phase -90', '-o');

%   NEW CLASS METHOD  -- PLACE IN  
%  “C:\DSP_ML\@TimeDomainSignal\FrequencySpectrum_V1.m”
function [frequency, freqData] ...
                           = FrequencySpectrum_V1(signal, name, lineStyle)
  figure('Name', name);     hold on;
  freqData = fft(signal.timeData);
  deltaF = signal.samplingFrequency / signal.numPts;
  frequency = deltaF *(0 : (signal.numPts - 1));
 
  plot(frequency, real(freqData), strcat(lineStyle, 'k')); 
  plot(frequency, imag(freqData), strcat(lineStyle, 'r')); 
  ylabel('\fontsize{14}AMPLITUDE'); 
  xlabel('\fontsize{14}FREQUENCY SPECTRUM (Hz)'); 
  legend('REAL()', 'IMAG()', 'location', 'NorthEast'); 
  pause(3);     %    Short pause between displays
end

%   NEW CLASS METHOD  -- PLACE IN  
%  “C:\DSP_ML\@TimeDomainSignal\FrequencySpectrum_V2.m”
function [frequency, freqData] ...  
                         = FrequencySpectrum_V2(signal, name, lineStyle)
  figure('Name', name);     hold on;
  freqData = fft(signal.timeData);
   %!!!! Use inbult fftshift( )  method to swap the frequency data halves 
  freqData = fftshift(freqData);   
%       numPts / 2 negative frequencies
%       DC + numPts / 2 - 1 positive frequencies
  numPtsBy2 = signal.numPts / 2;
  deltaF = signal.samplingFrequency / signal.numPts;
  frequency = deltaF * (-numPtsBy2 : (numPtsBy2 - 1));
% lineStyle allows display of sampled values and /or connection lines
  plot(frequency, real(freqData), strcat(lineStyle, 'k')); 
  plot(frequency, imag(freqData), strcat(lineStyle, 'r')); 
  ylabel('\fontsize{14}AMPLITUDE'); 
  xlabel('\fontsize{14} FREQUENCY SPECTRUM (Hz)'); 
  legend('REAL()', 'IMAG()', 'location', 'NorthEast');     
  pause(3);     %    Short pause between displays
end

This is the combination of a DC signal (0 Hz) and a 4 Hz cosinusoidal signal with equal amplitudes. In version 1 of the frequency display code FrequencySpectrum_V1(), we directly generate the discrete Fourier transform, DFT, of the 1 second long DualSignal.timeData[] array using the built-in fast Fourier transform function, fft().

Figure 2a shows we generate something very different than the expected frequency spectrum with two equal-amplitude components at 0 Hz and 4 Hz. Yes, there is a signal component around 0 Hz, and various, rather than one, frequency components around 4 Hz, but where did the frequency components around 28 Hz come from? We have now experienced our first of many sneak attacks from Murphy DSP-Law IV!

DSP_ML-IV: During DSP, sometimes nothing is as real as it seems!

The false signals around 28 Hz in Figure 2a arises from a misinterpretation of the results from the discrete Fourier transform operation (DFT). At a sampling rate of 32 Hz, the displayed output frequencies don’t cover the range from 0 Hz to 32 Hz. Instead they cover the range from –(32/2) = -16 Hz to nearly +(32/2) = +16 Hz.

In addition, the way the efficient implementation of the DFT using the fft() places the negative frequency components above the positive frequency components means that the signal components around 28 Hz are actually the results from frequencies around -4 Hz. Figure 2b shows the corrected display using the FrequencySpectrum_V2() script. This first uses the inbuilt fftshift() function to flip the two parts of the frequency spectrum, and then renumbers the frequency axis.

When correcting this display, we have invoked the fifth DSP Murphy’s Law for the first of many DSP occasions

DSP_ML-V: Applying fftshift() at least once every day, will help keep frequency confusion at bay!

The splitting of the DualSignal(t) signals into components around 0 Hz, 4 Hz and
–4 Hz can be understood as follows. DSP theory states mathematically the cosinusoidal signal:

can be expressed as the sum of two equal half-amplitude complex exponentials:

This means that we should expect that the time domain term:

will appear displayed at frequencies +4 Hz and – 4 Hz with amplitudes A/2. The DC term splits into two components at +0 Hz and -0 Hz., which “recombine” to one component at full amplitude.

SNEAKS ATTACKS FROM LAWS I AND III

The third problem still to solve is the following: Why do we get a spread of frequencies in the regions around +4 Hz and -4 Hz rather than the expected, very sharp frequency spikes from the following signal?

The answer is we have forgotten the DSP_ML I message about “dealing with discrete signals during DSP calculations”, which then allows a third DSP Murphy’s Law variant to remind us we had also forgotten to “avoid being one point out” during DSP calculations.

We sampled a signal for one second. A sampling rate of 32 Hz implies a time between samples of ∆T=1/32s. If we start sampling at time t =0 we get up to 33 values stored in timeSampled[] and timeData[] arrays. However, ML-I reminds us that DSP signals consist of a finite number of sampled values! Each of those 33 values are providing DualSignal characteristics over a time ∆T. This means we have not determined the characteristics of a signal of duration 1 second but the characteristics of a signal of duration 1/32 s longer, 33∆T=1.03125s!

Figure 2c shows that we get the expected sharp peaked frequency information for cosine signal 1+cos(2π4t) if we change the sampling time information obj.timeSampled, line 13, in the class method TimeDomainSignal, Listing 2, from 33-points to 32-points as follows to represent a signal sampled ever 1/32 s for exactly 1 second in duration:

function obj = TimeDomainSignal(Amp, 
 freq, phase, duration, samplingFreq)
  obj.samplingFrequency = samplingFreq;
  deltaT = 1 / samplingFreq;

% Change obj.timeSampled to use 32 points not 33
% obj.timeSampled = 0 : deltaT : duration – Line 13
 obj.timeSampled = 0 : deltaT : (duration - deltaT);
 obj.numPts = length(obj.timeSampled);
 obj.timeData = Amp * ...
   cos(2*pi * (freq * obj.timeSampled +
   phase / 360));
end

Figure 2d shows similar sharp peaked frequency information from a 1+sin(2π4t) signal. This signal is the difference of two complex exponentials:

The division by the 2j gives a sin() spectrum composed only of imaginary frequency components, whereas the cos() spectrum, Figure 2c, is made up of only real components.

While generating the expected sharp-peaked frequency displays for the 1+cos(2π4t) and 1+sin(2π4t) signals we have actually experienced several more variations of the one point-out DSP_ML-III. After applying fftshift() the DC component is at frequency[] array location numPts / 2 + 1, not at the halfway location numPts / 2. This means that there are numPts / 2 negative frequencies (-16 to -1 Hz), a DC (0 Hz) signal, and numPts / 2 – 1 positive frequencies
(1 to +15 Hz).

We must take a look at one of the underlying characteristics of the DFT if we want to understand the big change between the appearance of the 33-point frequency spectrum, Figure 2b, and that of the 32-point spectrum, Figure 2c.

Figure 3a shows 3 copies of the 32-point, 1+sin(2π4t), signal placed end to end along the time axis. The black line is the original 32-point 1+sin(2π4t) signal. The red lines show copies of this signal placed one data length, 32∆T seconds, up and down in time. Figure 3b shows 3 copies of the 33 point, 1+sin(2π4t) signal. The black line is the original 33-point 1+sin(2π4t) signal with the red lines representing copies of this signal shifted by the length of this data set,
33∆T seconds, in time.

In Figure 3a we have drawn a blue dotted line from the right end of the lower time copy to the left end of the original signal, and from right end of the original signal to the upper copy’s left end. When we draw these lines, we are emphasizing the cyclic properties inherent (meaning you can’t get around this property) in performing DFT frequency analysis. With these lines drawn, the 3 copies of the 32-point 1+sin(2π4t) signal seem to join together smoothly. By comparison, the copies of the 33-point 1+sin(2π4t) signal has a bump in the transitions between them in Figure 3b.

If you perform this multiple copy display analysis on any signal you have captured then you will find that “smooth merging of time shifted copies” will lead to a sharped peaked spectra, for example Figure 2d. If the merging is discontinuous, then the spectra will be distorted, for example Figure 2b.

Another way of looking at this problem is to think about the 3 copies of the original signal as one longer signal to be transformed into the frequency domain. The longer smooth signal, Figure 3a, looks like a single sinusoid—so it will lead to sharp peaks in the original 32-point spectrum, Figure 2d. The longer Figure 3b signal, will need extra frequency components to characterize those bumps so it will lead to a wide range of additional frequency components around the expected spectral peaks, such as in
Figure 2b.

Here’s the lesson to take away: If you have a signal with obvious (large) sinusoidal components, then measuring EXACTLY a whole number of cycles for each component will lead to sharp frequency peaks. Hint: You can always measure for a longer time and then discard the extra signals after you have decided how to do the analysis. But what if you have many components so it is impossible to measure a whole number of cycles for each component? For the answer to that question you’ll need to wait for the discussion on DSP data windowing in a later article.

MURPHY’S DSP LAW VI

Let’s move onto the frequency spectrum of a signal that is a little more realistic than the simple sum of pure single-tone sinusoids used in Figures 2 and 3. Figure 4a shows a short burst of a sinusoidal signal situated in the quarter middle section of the 1 second sampling period generated by the code given in Listing 4.

The DSP literature says that such a burst signal should be considered as a
1 second long sinusoidal signal multiplied by a ¼ second long square window (shown in blue). This means that the spectrum of this combination is theoretically the DFT of the first signal convolved with the DFT of the second signal.

The first spectrum is two spikes along the imaginary frequency axis as shown in
Figure 2d. Convolving this with the spectrum of the second signal, a sin(kf) / kf or sinc() function, should lead to the combined spectrum shown in blue in Figure 4b. However, that is not what you get—you have just encountered the sixth DSP Murphy’s Law.

DSP_ML-VI: By symmetry, anything that causes a problem during DSP frequency domain analysis is bound to be introducing a totally equivalent, but probably not immediately obvious, problem during DSP time domain analysis—and vice-versa!

If you examine the spectrum in Figure 4b, and apply a little bit of Sherlock Holmes’ intuition, you can see that the red signal is actually the expected blue signal multiplied by the alternating sequence +1, -1, +1, -1, +1, -1 …. Somehow we have managed to generate a convolution with a high frequency version of the sinc() following application of the fft() function during Listing 3’s FrequencySpectrum_V2().

To solve frequency domain problems generated by the DSP_ML-IV and V Laws shown in Figure 2a, we had to take into account special fft() properties. This required flipping the calculated two halves of the frequency data series by applying the fftshift() operator AFTER applying the fft() in Listing 3. We now need to also flip the initial time domain series end-over-end by applying a second fftshift() operation to solve the new phase Figure 4 problems introduced during DSP_ML-VI, this time BEFORE applying the fft() rather than after. This operation is implemented in the (hopefully) final FrequencySprectrum() method given in Listing 5.

Figure 5a shows the expected result. The final spectrum is two spikes entirely along the imaginary frequency axis from the sin() wave convolved with the spectrum of the second signal, a sinc() function with no real components.

% Place in file “C:\DSP_ML\Investigate_DSP_ML_VI.m”

function Investigate_DSP_ML_VI( )
  clear classes;   %  Force MATLAB / OCTAVE to see class changes   
    fastSampling = 128; duration = 1; phase0 = 0; phase90 = -90;
  %% Generate burst of high frequency signal with phase -90
  signal32HzPhase90 = TimeDomainSignal(1.0, 32, phase90, ...
                                               duration, fastSampling);
  startFraction = 3/8;    endFraction = 5/8;  % Keep centre quarter
  [burstSignal32HzPhase90, burstEnvelope] = ...
        GenerateBurstSignal(signal32HzPhase90,  startFraction, endFraction);                                                                             
  DisplayTimeSignal(burstSignal32HzPhase90, ...
                              'Burst 32 Hz, Phase -90', '-o');
  fprintf( ' Adding implied data window\n');
  plot(burstSignal32HzPhase90.timeSampled, burstEnvelope,  ...
                                                    '-b', 'linewidth', 3);
  legend('BURST', 'IMPLIED DATA WINDOW', 'location', 'NorthEast');
  pause (3);

%%  Time shifted frequency display -- generate first to get envelope
  [frequency, timeShiftedFreqDataPhase90] = ...
      FrequencySpectrum(burstSignal32HzPhase90, ...
        'DUMMY SHAPE NEEDED FOR LATER', '-');
%% Burst spectrum should involve sinc functions, but does not
  FrequencySpectrum_V2(burstSignal32HzPhase90, ...
          'Frequency Spectrum Burst 32 Hz Phase 90 - No time shift', '-');
%%  Overlay with true sinc functions  
   plot(frequency, imag(timeShiftedFreqDataPhase90), '-b', 'linewidth', 3);
   pause (3);
   fprintf( 'Added expected frequency shape\n');
%%  Time shifted input shows correct frequency display 
  [frequency, timeShiftedFreqDataPhase90] = ...
           FrequencySpectrum(burstSignal32HzPhase90, ...
           'Frequency Spectrum Burst 32 Hz Phase 90 - time shifted', '-');
%% 16 Hz signal burst with phase 0 has distorted spectrum
  signal32HzPhase0 = TimeDomainSignal(1.0, 32, phase0, ...
                                                duration, fastSampling);
  startFraction = 3/8;    endFraction = 5/8;  % Keep centre quarter
  [burstSignal32HzPhase0, burstEnvelope] = ...
        GenerateBurstSignal(signal32HzPhase0, startFraction, endFraction); 
  [frequency, timeShiftedFreqDataPhase0] = ...
           FrequencySpectrum(burstSignal32HzPhase0, ...
           'Frequency Spectrum Burst 32 Hz Phase 0 - With time shift', '-');

% Place in file  “C:\DSP_ML\@TimeDomainSignal\GenerateBurstSignal.m”

function [burstSignal, burstEnvelope] =   ...
                 GenerateBurstSignal(signal, startFraction, endFraction) 
% During FFT fraction 0.5 (1/2) point was located at numPts / 2 + 1
% So be consistent when locating the starting point of the signal burst 
  burstEnvelope = zeros(signal.numPts, 1);
  startPoint = floor(signal.numPts * startFraction + 1); % Keep start point  % DSP_ML III - Avoid one point out
  endPoint = floor(signal.numPts * endFraction + 1) - 1; %Exclude endpoint  burstEnvelope(startPoint: endPoint) = 1;

  burstSignal = signal;      % Keep portion [startFraction, endFraction)
  burstSignal.timeData =  burstSignal.timeData .* burstEnvelope';
end

% Place in file “C:\DSP_ML\@TimeDomainSignal\FrequencySpectrum.m” 

function [frequency, freqData] ...                     %  NEW CLASS METHOD
                            = FrequencySpectrum(signal, name, lineStyle)
  figure('Name', name);     hold on;
  % DSP_ML VI -- Frequency phase incorrect unless fold time data
  signal.timeData = fftshift(signal.timeData);  %%% Swap time data halves
  freqData = fft(signal.timeData);
 
% Frequency spectrum incorrect unless you fold final frequency data
% There are        numPts / 2 negative frequencies
%              DC + numPts / 2 - 1 positive frequencies
  freqData = fftshift(freqData);            %%% Swap frequency data halves
  numPtsBy2 = signal.numPts / 2;
  deltaF = signal.samplingFrequency / signal.numPts;
  frequency = deltaF *(-numPtsBy2:(numPtsBy2 - 1));

  plot(frequency, real(freqData), strcat(lineStyle, 'k'));  hold on;
  plot(frequency, imag(freqData), strcat(lineStyle, 'r'),  'linewidth', 2); 
  ylabel('\fontsize{14}AMPLITUDE'); 
  xlabel('\fontsize{14}FREQUENCY SPECTRUM (Hz)');  
  legend('REAL()', 'IMAG()', 'location', 'NorthEast'); 
  pause(3);   
end

MORE SNEAK ATTACKS

Figure 5b shows the frequency spectrum of a short burst of a cos() signal placed in the middle quarter of the 1 second sampling period—a mere 90-degree phase shift from the sin() situation in Figure 5a. Yes, there is the expected convolution of a sinc function with the two peaks placed along the real frequency axis, shown in black. However, the frequency spectrum is not entirely real. There is now an additional, strong, sinusoidal imaginary frequency component, shown in red.

Misquoting the palace guard in Shakespeare’s Hamlet, we can say: “There is something still wrong in the state of our DSP analysis.” But the problem to solve is not simple. Where did we go wrong…again? We can do our first sleuthing by rewriting DSP_ML-VI as: “Whatever happens (or goes-wrong) in the frequency domain can happen (or go wrong) in the time domain!”

If we have a pure cosinusoidal signal in the time domain, then that will introduce two spikes in the frequency domain because of the DSP property:

Applying DSP_ML-VI, means that the appearance of an EXTRA cosinusoidal pattern in the frequency spectrum implies that there are two EXTRA spikes, isolated samples, in the time domain. We have suffered two further sneak attacks of DSP_ML III one-outness and gathered two points more data than we intended.

We can go further in our sleuthing. The signal in the frequency domain has 16 full cycles indicating that the extra points we did not intend to sample were collected at times ±16∆T.

LAW IV REVISITED

If we unintentionally are using 2 extra points to produce Figure 5b, then we must have ALSO unintentionally captured 2 extra points when producing Figure 5a! So where is the extra sinusoidal component in Figure 5a, which should be present if we apply the equal opportunity for doing things wrong DSP_ML-VI! The answer is hidden in the rewritten:

DSP_ ML-IV: During DSP, sometimes nothing is as real or correct as it seems!

The extra component is present in Figure 5a, but by luck, the extra captured points happened to have zero values, so they generate extra frequency components of zero amplitude! Actually, we should say “unfortunately”, rather than “by luck”, the signal components are not seen. We could have released our DSP code out into the real world to cause total confusion to our customers when they applied this script to their own signals and got unexpected results. It’s always better to recognize a problem is present than have the problem hidden from you.

We will leave it to our readers to fix the script GenerateBurstSignal.m AND all the other sneak attacks from the DSP Murphy’s Law that perhaps will not be uncovered in our code until later articles in this series, or maybe never!

TOOLS SMOOTH THE WAY

All the instances of Murphy’s Law we have examined in this article are not appearing because we are using MATLAB or Octave scripting languages. These tools make it easier to visualize the problems, but the issues are associated with inherent properties of the DSP data points themselves. They can reappear whether we implement the DSP algorithms in C, C++, C#, python, Java and so on.

DSP_ML-III is particularly important when you move between algorithms programmed using a simulation tool language to real-time implementation on hardware. When your results don’t pan out, you’ve probably calculated something just one sample point out! Remember, C++ arrays start at index 0 while MATLAB and Octave arrays start at index 1.

GETTING STARTED

In this section, we outline how to set up the directories needed for the DSP_ML project, install either Octave or MATLAB, and check the installation. For your convenience, these instructions are repeated on the Circuit Cellar article materials webpage along with all the needed web links. By the end of the article you will have a number of main scripts, such as Investigate_DSP_I_II_III.m, in the directory “C:\DSP_ML” and class methods such as “TimeDomainSignal.m” and variants of the frequency analysis scripts, “FrequencySpectrum_V1.m”, in the class directory, “C:\DSP_ML\@TimeDomainSignal”, as shown in center left windows in Photo 2 and Photo 3.

• Photo 2
• Photo 3

Generate a folder “C:\DSP_ML” to hold the main DSP_ML scripts and a second folder, “C:\DSP_ML\@TimeDomainSignal to hold the TimeDomainSignal project class methods (functions) used to generate and display example DSP signals and procedures.

Install Octave or MATLAB following the instructions here. Type the main directory name “C:\DSP_ML” into the current directory window, top window in Photo 2 and
Photo 3, once you activate your installed program.

This series of articles is designed as a series of small projects that can each be done in an hour on different days. You can get all the code from the Circuit Cellar Article Code & Files webpage. Alternately cut-and-paste scripts directly from the pdf. If you take the pdf cut-and-paste approach watch out as some syntax may not carry over properly into the MATLAB / Octave script, especially the single quotes and any solid lines separating listing parts.

To create a script, start by clicking on “File | New | New Script” on the top left corner of the main page of the Octave screen, or “File | New” in MATLAB.

For many reasons and under circumstances that would fill several web-pages—for example, storing DSP_ML directory on a USB stick with Octave—we recommend these steps:

1. Right click on Octave / MATLAB icon and select properties. Change Start in option from %USERPROFILE% to “C:\DSP_ML”
2. Add script startup.m (Listing 6a) to directory “DSP_ML”. Then, close and restart the program.

To check the install, cut-and-paste the code for Listing 6b from the code listing from the Circuit Cellar Article Code & Files webpage into the editor window, and save as the script “C:\DSP_ML\HelloWorld.m”. You can run the test code, and generate a graphic plot, by typing HelloWorld in the command window; bottom left in Octave, bottom center in MATLAB.

a)
function startup( )
   addpath('C:\DSP_ML');
   clear classes;

__________________________________________________________________

b)
function HelloWorld( )
   fprintf('Hello, World!\n');
 
   close all;   % Close any previous figures
   figure('Name', '2  Lines');
   plot(1:10, '-ok'); hold on;
   plot(2:11, '-r');
 
   legend('FIRST LINE',  'SECOND LINE');
end  

Listing 6 (a) Recommend script “C:\DSP_ML\startup.m” (b) This script “C:\DSP_ML\HelloWorld.m” can be used to check the Octave or MATLAB installation.

Installing Octave: Create a folder “C:\DownLoads” and go to the link provided on the Circuit Cellar article materials webpage to download a version of Octave, for example. “octave-4.2.1-w64-installer.exe” for a 64-bit system. Run your downloaded .exe file to get the installation started. Click on various “Nexts“ to accept the license agreement. At one stage, you will be asked to provide a directory to install the program, for example “C:\Octave”. Be patient when various command windows appear, they will disappear.

To test the installation, click on the “Octave GUI” icon on your desktop. With the Start In option change, the current directory will be “C:\DSP_ML”. Cut-and-paste in the HelloWorld() code and run.

Restarting Octave solves an inconsistent problem with Octave not recognizing newly added class methods.

Installing MATLAB: Create a folder “C:\DownLoads” and go to the link provided on the Circuit Cellar article materials webpage to click on the “Test drive MATLAB” link to get the 30-day MATLAB free trial. You will need to supply basic personal information such as your email, location and planned usage. Another link on the Circuit Cellar article materials webpage provides full details of the installation steps. To test the installation, click on the “MATLAB” icon on your desktop. With the Start In option change, the current directory will be “C:\DSP_ ML”. Cut-and-paste in the HelloWorld() code and run. 

Authors’ Note: Thanks to summer undergraduate researchers Maaz Khurram and Jacob Ridgway for their help in finalizing the article.

DSP_ML Project Files can be cut and pasted from the .pdf file directly into your script file .m or can be found here:  github.com/SMILE-Projects/SMILE_Projects

Mathworks | www.mathworks.com
GNU Octave | www.gnu.org/software/octave

PUBLISHED IN CIRCUIT CELLAR MAGAZINE • JUNE 2018 #335 – Get a PDF of the issue

Michael Smith, Mai Tanaka and Ehsan Shahrabi Farahani

+ posts

Mike Smith (Mike.Smith@ucalgary.ca) is a professor in Computer Engineering at the University of Calgary, Canada. When not singing with Calgary’s adult recreational choir Up2Something, Mike’s main interests are in developing new biomedical engineering algorithms and moving them onto multi-core and multiple-processor embedded systems in a systematic and reliable fashion. Mike has recently become a convert to the application of agile methodologies in the embedded environment. Mike was one of the notable contributors interviewed in Circuit Cellar’s 25th Anniversary issue.

Mai Tanaka (tanakam@ucalgary.ca.) is a graduate student in Electrical and Computer Engineering at the University of Calgary, Canada. Her research interests include the transfer and translation of information at the human machine interface such as speech processing.

Ehsan Shahrabi Farahani (ehsan.shahrabi.f@gmail.com) has completed his Masters in Electrical and Computer Engineering at the University of Calgary, Canada. His research interests include digital signal processing, image processing, biomedical signal processing, medical data analysis, medical imaging and speech processing.

• Michael Smith, Mai Tanaka and Ehsan Shahrabi Farahani

  https://circuitcellar.com/author/michael-smith-mai-tanaka-and-ehsan-shahrabi-farahani/

  Murphy's Laws in the DSP World (Part 3)
• Michael Smith, Mai Tanaka and Ehsan Shahrabi Farahani

  https://circuitcellar.com/author/michael-smith-mai-tanaka-and-ehsan-shahrabi-farahani/

  Murphy's Laws in the DSP World (Part 2)