Monday 27 April 2020

Python/Numpy : How Not To Generate A Sinusoidal Waveform

I was recently reviewing some Python/Numpy code that included a waveform generator. A simplified version of code looked like this :

x = np.linspace(0,2*np.pi-(2*np.pi/8),8)
np.sin(x)

This generates the following :

array([ 0.00000000e+00,  7.81831482e-01,  9.74927912e-01,  4.33883739e-01,
       -4.33883739e-01, -9.74927912e-01, -7.81831482e-01, -2.44929360e-16])

Which looks like a perfect single cycle of a sinusoid. Except it isn't !
On closer inspection, the last element in the array is, to all intents and purposes, 0, which means that this isn't a perfect single cycle of a sinusoid because that final sample is actually the first sample of the next cycle.

To generate a perfect single cycle of a sinusoid using linspace you need to account for where the last sample of the sinusoid should fall, if you were to plot it on a graph.

x = np.linspace(0,2*np.pi-(2*np.pi/8),8)
np.sin(x)

This generates the following array, which is spot on :

array([ 0.00000000e+00,  7.07106781e-01,  1.00000000e+00,  7.07106781e-01,
        1.22464680e-16, -7.07106781e-01, -1.00000000e+00, -7.07106781e-01])

In thinking about this problem, it occurred to me that this is not ideal and very likely to cause confusion becasue it is easy to forget the required modification. The main reason for the confusion is that standard Python generates and processes data from, for example, 0 to N-1 as shown in this simple Numpy example :

np.arange(8.)

Which yeilds :

array([0., 1., 2., 3., 4., 5., 6., 7.])

So returning to the original problem, a far safer way of generating the sinusoid is the following code :

x = np.arange(0., 2.*np.pi, 2.*np.pi/8.)
np.sin(x)
Which generates the following array :
array([ 0.00000000e+00,  7.07106781e-01,  1.00000000e+00,  7.07106781e-01,
        1.22464680e-16, -7.07106781e-01, -1.00000000e+00, -7.07106781e-01])

Now we have the first np.arrange() instruction to generate the time index and the second stage np.sin() to generate the sinudoid. This is clear, precise and unlikely to cause error.

Side note : of course, it would be entirely possibl to combine this into a single line instruction however I believe this causes other possibilities for error insertion.



Saturday 25 April 2020

The Difference Between FFT Spectrum and Power Spectral Density

I always teach the difference between FFT Spectrum and Power Spectral Density on my DSP courses and many students find it confusing.

This applications note from Audio Precision summarizes the subject very neatly : The Difference Between FFT Spectrum and Power Spectral Density

Functions for calculating both the FFT Spectrum and Power Spectral Density are included in the SigLib DSP Library.