Friday 8 March 2024

Plotting a Spectrogram In Python, Using Numpy and Matplotlib

When performing frequency domain (FFT) based processing it is often useful to display a spectrogram of the frequency domain results. While there is a very good SciPy spectrogram function, this takes time domain data and does all of the clever stuff. However if you are processing data in the frequency domain you often just want to build the spectrogram dataset and keep appending FFT results to it.

A spectrogram is a 3D plot, with the following configuration:

  • Time is on the X axis
  • Frequency is on the Y axis
  • Frequency magnitude is shown in colour

This program will use the Matplotlib function imshow() to display the spectrogram.

There are two key tricks to using imshow() for this purpose:

  • Rotate the dataset so that time is on the x-axis and frequency is on the y-axis
  • Scale the x and y axes labels to correctly show time and frequency
  • Remove the second half of the FFT results - once the magnitude of the FFT result has been calculated, the two halves of the result are mirror images so we can discard the upper half.

Here's the code:

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

# Global Configuration
plotXAxisSecondsFlag = True                             # Set to True to Plot time in seconds, False to plot time in samples
plotYAxisHzFlag = True                                  # Set to True to Plot frequency in Hz, False to plot frequency in bins

Fs = 10000                                              # Sampling Frequency (Hz)
timePeriod = 10                                         # Time period in seconds
sampleLength = Fs*timePeriod

sinusoidFrequency = 1000                                # Frequency of sine wave (Hz)

fftLength = 256                                         # Length of the FFT
halfFftLength = fftLength >> 1

window = np.hanning(fftLength)

time = np.arange(sampleLength) / float(Fs)              # Generate sinusoid + harmonic with half the magnitude
x = np.sin(2*np.pi*sinusoidFrequency*time) * ((time[-1] - time)/1000) # Decreases in amplitude over time
x += 0.5 * np.sin(2*2*np.pi*sinusoidFrequency*time) * (time/1000)     # Increases in amplitude over time

# Add FFT frames to the spectrogram list - Note, we use a Python list here becasue it is very easy to append to
spectrogramDataset = []

i = 0
while i < (len(x) - fftLength):                         # Step through whole dataset
  x_discrete = x[i:i + fftLength]                       # Extract time domain frame
  x_discrete = x_discrete * window                      # Apply window function
  x_frequency = np.abs(np.fft.fft(x_discrete))          # Perform FFT
  x_frequency = x_frequency[:halfFftLength]             # Remove the redundant second half of the FFT result
  spectrogramDataset.append(x_frequency)                # Append frequency response to spectrogram dataset
  i = i + fftLength

# Plot the spectrogram
spectrogramDataset = np.asarray(spectrogramDataset)     # Convert to Numpy array then rotate and flip the dataset
spectrogramDataset = np.rot90(spectrogramDataset)
z_min = np.min(spectrogramDataset)
z_max = np.max(spectrogramDataset)
plt.figure()
plt.imshow(spectrogramDataset, cmap='gnuplot2', vmin = z_min, vmax = z_max, interpolation='nearest', aspect='auto')
plt.title('Spectrogram')
freqbins, timebins = np.shape(spectrogramDataset)
xlocs = np.float32(np.linspace(0, timebins-1, 8))
if plotXAxisSecondsFlag == True:
  plt.xticks(xlocs, ["%.02f" % (i*spectrogramDataset.shape[1]*fftLength/(timebins*Fs)) for i in xlocs]) # X axis is time (seconds)
  plt.xlabel('Time (s)')
else:
  plt.xticks(xlocs, ["%.02f" % (i*spectrogramDataset.shape[1]*fftLength/timebins) for i in xlocs])      # X axis is samples
  plt.xlabel('Time (Samples)')

ylocs = np.int16(np.round(np.linspace(0, freqbins-1, 11)))
if (plotYAxisHzFlag == True):
  plt.yticks(ylocs, ["%.02f" % (((halfFftLength-i-1)*Fs)/fftLength) for i in ylocs])  # Y axis is Hz
  plt.ylabel('Frequency (Hz)')
else:
  plt.yticks(ylocs, ["%d" % int(halfFftLength-i-1) for i in ylocs])                   # Y axis is Bins
  plt.ylabel('Frequency (Bins)')
plt.show()

 

Friday 12 January 2024

The 33rd Annual Running Of The University Of Oxford Digital Signal Processing Course Will Be Held Online Again, In 2024

The course first moved online in 2020 and has received excellent reviews from the attendees.

The course will run from Monday 8th April to Friday 17th May 2024, with live online classes one afternoon per week.

Based on the classroom course, Digital Signal Processing (Theory and Application), this online course consists of weekly live online tutorials and also includes a software lab that can be run remotely. We'll include all the same material, many of the existing labs and all the interaction of the regular course.

Online tutorials are delivered via Microsoft Teams once each week and practical exercises are set to allow you to practice the theory during the week. 

You will also have access to the course VLE (virtual learning environment) to communicate with other students, view and download course materials and tutor support is available throughout.

Code examples will be provided although no specific coding experience is required. 

The live tutorials will be on Wednesday each week from 13:00 - 14:30 and 15:00 - 16:30 (GMT) with a 30-minute break in between.

You should allow for 10 - 15 hours study time per week in addition to the weekly lessons and tutorials.

After completing the course, you should be able to understand the workings of the algorithms we explore in the course and how they can solve specific signal processing problems.

Full details are available here: https://www.conted.ox.ac.uk/courses/digital-signal-processing-online.


Friday 29 December 2023

SigLib Now Includes Kalman Filtering Functions

Over the holiday period I decided to refresh my knowledge of Kalman filters by watching the excellent video series here: Kalman Filter YouTube Lessons.

Here is a diagram to show the architecture of the Kalman Filter:


As a result, I have now added 1D and 2D Kalman filters to the SigLib DSP Library, which can be found here: https://github.com/Numerix-DSP/siglib.

In the comments here https://www.youtube.com/watch?v=Fuy73n6_bBc&list=PLX2gX-ftPVXU3oUFNATxGXY90AULiqnWT&index=27, Gergely Takács mentioned his excellent Matlab and Python examples that can be found here: https://github.com/gergelytakacs/planeKF. Inspired by Gergely I have added the same example to the SigLib DSP Library examples.


Wednesday 8 November 2023

Deploying Matlab On ARM Using Codegen

Several years ago I worked on two projects that required code to be deployed on ARM devices, where the original algorithms had been developed in Matlab. I researched Matlab's Codegen capabilities and realized that it has some very useful features so I've generated a generic example, with no customer code, to demonstrate the capabilities.

Among the many great features of Codegen is Whole Project code optimization. For example, if a value is stored in a global variable then passed to the underlying generated C code function, Codegen will bypass passing the value on the stack and merely read it from the global memory pool, in the sub-function.

The project can be found here: https://github.com/Numerix-DSP/Matlab_To_C.

There are several tricks to achieve good quality code from Matlab/Codegen and I've covered a few of them in the documentation. There are also several areas that Codegen is not able to support directly like interrupt service routines on the embedded devices - The project shows how to integrate the code with embedded I/O routines.

The project is Windows based but with a few minor tweeks can run perfectly well under Linux or OSX. For testing the generated C code on a host, the project includes an example host program that shows how to handle file I/O in a Matlab compatible format and also uses SigLib for .csv file I/O and Gnuplot to replace the Matlab graphics plotting functions.

This project also includes a batch file to convert the Matlab code to C to allow it to run on an STM32 device, with the simulation files stored on a USM memory stick. The File I/O funcationality can easily be replaced to run the code from and Interrupt Service Routine instead.

For my NXP loving friends I have tried to run the program on an LPC55S69 EVK but unfortunately, I'm having a few problems with USB file I/O. Once that is resolved, I'll port to that device.

Following that I will add scripts to swap out the Matlab DSP function calls (e.g. FFT) and replace them with calls to the ARM CMSIS-DSP library.


Tuesday 18 July 2023

The 31st Annual Running Of The University Of Oxford Digital Signal Processing Course Will Be Held Online Again, In 2023

 The 31st annual running of the University Of Oxford Digital Signal Processing course will include an online edition that will be running over a six week period, from Monday 25th September to Friday 3rd November 2023.

The course first moved online in 2020 and has received excellent reviews from the attendees.

Based on the classroom course, Digital Signal Processing (Theory and Application), this online course consists of weekly live online tutorials and also includes a software lab that can be run remotely. We'll include all the same material, many of the existing labs and all the interaction of the regular course.

Online tutorials are delivered via Microsoft Teams once each week and practical exercises are set to allow you to practice the theory during the week. 

You will also have access to the course VLE (virtual learning environment) to communicate with other students, view and download course materials and tutor support is available throughout.

Code examples will be provided although no specific coding experience is required. 

The live tutorials will be on Wednesday each week from 13:00 - 14:30 and 15:00 - 16:30 (GMT) with a 30-minute break in between.

You should allow for 10 - 15 hours study time per week in addition to the weekly lessons and tutorials.

After completing the course, you should be able to understand the workings of the algorithms we explore in the course and how they can solve specific signal processing problems.

Full details are available here: https://www.conted.ox.ac.uk/courses/digital-signal-processing-online

Monday 19 June 2023

Simple Python Data Plotter For DSP Log Files

I often log DSP data to .log files, which are typically text files with comma separated columns. In order to plot the data I use a variant of the following Python file, which extracts the data columns and plots the results.

For easier human reading, I typically use ", " for separating the columns, rather than just a single ",". This requires the use of the separator specifier and the use of the 'python' engine because the 'c' engine does not support regex separators.

The nice thing about this is that human readable text can be inserted between the columns of data.

The main thing to be aware of is that all rows in the file must include the same number of columns otherwise read_csv() fails to read the file correctly.

To use this program:

python plot.py log_file.log


# Script to plot the data in the extracted columns of the log file
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys

def main(argv):
  print(argv[0])

  df = pd.read_csv(argv[0], header=None, sep=', ', engine='python')
  df = df.head(50000)                                       # Truncate samples for better plotting

  f, ax = plt.subplots()

  inputSamples = df.iloc[:,1]
  plt.plot(inputSamples, label="inputSamples")

  gainedSamples = df.iloc[:,3]
  plt.plot(gainedSamples, label="gainedSamples")

  meanSquaredSum = df.iloc[:,7]
  plt.plot(meanSquaredSum, label="meanSquaredSum")

  ymin, ymax = ax.get_ylim()
  ax.set_yticks(np.round(np.linspace(ymin, ymax, 8), 2))    # Print 8 y-axis ticks

  plt.legend(loc="upper left")
  plt.show()

if __name__ == "__main__":
   main(sys.argv[1:])




Monday 20 March 2023

Are Chat-GPT and Google Bard The New Frontier For Writing DSP Code?

I needed to take a break from doing some consultancy work and writing code for customers but the weather is very British today so rather than brave the outdoors I thought I'd have a play with Chat-GPT and see how good it is at writing DSP code.

I asked Chat-GPT the following: "write C code for FIR filter".

Here are the results, along with the supporting comments provided by Chat-GPT:

Here is an example code for a Finite Impulse Response (FIR) filter implementation in C language. This implementation assumes that the input signal samples and filter coefficients are stored in arrays.

#include <stdio.h>
#include <stdlib.h>

#define
N 100  // length of input signal and filter coefficients
#define M 10   // length of filter

int main()
{
    float input_signal[N];    // input signal
    float filter_coeff[M];    // filter coefficients
    float output_signal[N];   // output signal after filtering
    int i, j;

    // Generate input signal and filter coefficients
    for (i = 0; i < N; i++) {
        input_signal[i] = (float)rand() / RAND_MAX;   // generate random input signal samples
    }
    for (i = 0; i < M; i++) {
        filter_coeff[i] = (float)rand() / RAND_MAX;   // generate random filter coefficients
    }

    
// Apply FIR filter
    for (i = 0; i < N; i++) {
        output_signal[i] = 0.0;
        for (j = 0; j < M; j++) {
            if (i - j >= 0) {
                output_signal[i] += input_signal[i - j] * filter_coeff[j];
            }
        }
    }

    // Print output signal
    for (i = 0; i < N; i++) {
        printf("%f ", output_signal[i]);
    }

    return
0;
}

In this implementation, we first generate random input signal samples and filter coefficients using the rand() function. Then, we apply the FIR filter by iterating over each sample of the input signal and computing the output sample as a weighted sum of the current and previous input samples, using the filter coefficients. Finally, we print the output signal.

Note that this implementation assumes that the length of the input signal and filter coefficients are fixed, and that the input signal is buffered in memory. In practice, more sophisticated algorithms and data structures may be required to handle larger input signals or streaming data.

Chat-GPT was similarly impressive for "write C code for IIR filter".

Writing DSP code always gets more fun when trying to write a more complex function, such as the Fast Fourier Transform.

It is very interesting to note that each Chat-GPT response is different and, in the case of the FFT example, I had to ask it to generate 3 responses before having a solution that was semi-optimal and included the complete code for the FFT.

  • Solution #1 was a very interesting yet very impractical recursive solution!
  • Solution #2 gave a solution that used the FFFTW library, which is more like it but I was hoping for raw C code
  • Solution 3 gave the following solution:

#include <stdio.h>
#include <math.h>

#define PI 3.14159265358979323846

void fft(double *x_re, double *x_im, int n) {
    int i, j, k;
    double c, s, t1, t2;
    double tr, ti;
    int m = log2(n);

    // bit reversal permutation
    j = 0;
    for (i = 0; i < n - 1; i++) {
        if (i < j) {
              tr = x_re[i];
              ti = x_im[i];
              x_re[i] = x_re[j];
              x_im[i] = x_im[j];
              x_re[j] = tr;
              x_im[j] = ti;
        }
        k = n / 2;
          while (k <= j) {
              j -= k;
              k /= 2;
        }
        j += k;
    }

   
 // Cooley-Tukey radix-2 algorithm
    for (i = 0; i < m; i++) {
        int len = pow(2, i);
        for (j = 0; j < n; j += 2 * len) {
              for (k = 0; k < len; k++) {
                  c = cos(-2 * PI * k / (2 * len));
                  s = sin(-2 * PI * k / (2 * len));
                  t1 = x_re[j + k + len] * c - x_im[j + k + len] * s;
                  t2 = x_im[j + k + len] * c + x_re[j + k + len] * s;
                  x_re[j + k + len] = x_re[j + k] - t1;
                  x_im[j + k + len] = x_im[j + k] - t2;
                  x_re[j + k] += t1;
                  x_im[j + k] += t2;
              }
        }
    }
}


While solution 3 includes the complete code, it does call functions from the standard math library so is far from optimal.

The following request solved the math library issue "write C code for applying the FFT to a 64 point sine wave using a look-up table for the twiddle factors". Try it for yourself, the code is very good but if we want a real FFT then we need to tweak the request slightly.

The following request "write C code for applying the real FFT to a 64 point sine wave using a single look-up table for the twiddle factors and return the complex result" generated a very useable piece of code although it still uses the standard math pow function:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

#define PI 3.14159265358979323846
#define N 64

void rfft(double *x, int n, double complex *twiddle_factors, double complex *y) {
    int i, j, k;
    double c, s, t1, t2;
    double complex tmp;

    
// bit reversal permutation and complex conversion
    j = 0;
    for (i = 0; i < n - 1; i++) {
        if (i < j) {
            tmp = x[i];
            x[i] = x[j];
            x[j] = tmp;
        }
        y[i] = x[i] + I * 0.0;
    }
    y[n-1] = x[n-1] + I * 0.0;

    // Cooley-Tukey radix-2 algorithm
    for (i = 0; i < log2(n); i++) {
        int len = pow(2, i);
        int half_len = len / 2;
        int table_step = N / len;
        for (j = 0; j < n; j += 2 * len) {
            for (k = 0; k < half_len; k++) {
                c = creal(twiddle_factors[k * table_step]);
                s = cimag(twiddle_factors[k * table_step]);
                tmp = c * y[j + k + half_len] - s * conj(y[j + k + half_len]);
                y[j + k + half_len] = y[j + k] - tmp;
                y[j + k] += tmp;
            }
        }
    }

    // output conversion
    for (i = 1; i < n/2; i++) {
        y[n-i] = conj(y[i]);
    }
    y[n/2] = creal(y[n/2]);
}

int main() {
    double x[N];
    double complex y[N/2+1];
    double complex twiddle_factors[N/2];
    int i;

    
// generate input signal (64-point sine wave with frequency 4 Hz)
    for (i = 0; i < N; i++) {
        x[i] = sin(2 * PI * i * 4 / N);
    }

    
// generate twiddle factors
    for (i = 0; i < N/2; i++) {
        twiddle_factors[i] = cexp(-2 * PI * I * i / N);
    }

    
// perform RFFT
    rfft(x, N, twiddle_factors, y);

    // print output
    for (i = 0; i < N/2+1; i++) {
        printf("%d: %f + %fi\n", i, creal(y[i]), cimag(y[i]));
    }

    return 0;
}

While the use of the pow function is easily resolved, this does require a manual tweak to the code.

My final C test used the following request: "Write some C code to apply a 20 tap FIR filter to a sine wave and then analyze the results with a fast fourier transform". This resulted in some very neat code that also included a Hamming window on the input data and calculated the magnitude of the complex FFT result, all be it on the 4th attempt.

I did wonder how much knowledge Chat-GPT has of the SigLib DSP library so I tried "write C code for applying the real FFT to a 64 point sine wave using the SigLib DSP library". There were a few minor errors in the code output but generally, it was found that Chat-GPT had a pretty good knowledge about the library. One of the outputs actually included a full description of each step in the generated code :-).

When I'm writing completely new code, that is not based on anything I have written previously, I typically head over to Python (or Matlab, if that is what the customer prefers). So I asked Chat-GPT "Write some Python code to apply a 20 tap FIR filter to a sine wave and then analyze the results with a fast fourier transform". This code was very neat because it basically generated a very small but useable piece of code that used scipy.

I also tried the same request with Matlab instead of Python and I'm glad to report that Chat-GPT understands both languages, in addition to C :-)

Finally, I thought I could fool the beast by replacing the "Python code" request with "Analog Devices assembly code" and "Texas Instruments assembly code". It is interesting that the code generated is, indeed, assembly language for the appropriate DSPs but generates calls to the more complex functions as shown in this snippet:

_fft:
        ; perform FFT using TI's DSP library
        call    #_do_fft


Google Bard

Since writing the original blog, I've had some time to test Google Bard. The results are almost identical to Chat-GPT. The code generated looks like it has come from the same sources, which I guess is not surprising really. I also found that Bard would generated very different results for each request and often required several attempts to get the perfect solution.

Bing Chat

Microsoft have their own version of Chat-GPT that can be used via Bing Chat however there are some differences to Chat-GPT and Bard.

Bing Chat doesn't seem to allow you to re-submit the same question because it seems to be "smart enough" to know that it has already answered the question and will tell you that it has already provided an answer. This is unfortunate if the provided solution is sub-optimal.

Bing Chat also only works in Edge. When using Bing in another browser it just acts like a traditional search engine - not a major limitation but one to remember.

VSCode Integration

For programming, the best option I've found is to install CodeGPT in VSCode, along with the API Key for Chat-GPT. Now all that is required is to write the comment for your code and hit <ctrl>+<shift>+i and start configuring the code that is generated :-).

Conclusion

For simple functions such as FIR and IIR filters, the code generated by Chat-GPT and Google Bard does exactly what was asked however the code is not really useable for embedded DSP applications. The main reasons are that the bots don't necessarily separate initialization and run-time code and also the code is not always the most optimized for the majority of CPU/DSP architectures.

When using Chat-GPT or Bard, it is often necessary to initiate several requests before the most optimum solution is provided.

For more complex functions such as the FFT, it is necessary to have a lot more knowledge of the application requirements and what optimizations can be performed. For example, when performing the FFT of a sine wave the user needs to understand that it is not necessary to use the complex FFT.

This code was generated using Chat-GPT v3.5, it will be very interesting to perform the same requests with Chat-GPT V4.0, when that is publicly available.

Final Thought

I look forward to the day when I can ask one of these new bots to optimize my new super-duper DSP function and, British weather permitting, I can spend the rest of the day on my mountain bike or in my canoe :-)

Copyright © 2023 Sigma Numerix Ltd.