Monday, 27 April 2026

The 35th Annual Running Of The University Of Oxford Live Digital Signal Processing Course, In May 2026

 The 35th annual running of the University Of Oxford Live Digital Signal Processing course will be running in Oxford, UK, from Tuesday 20th to Friday 23rd May 2026.

The courses are presented by experts from industry for Engineers in industry and over the last 30 years has trained many hundreds of Engineers, from all areas of Science and Engineering.

Here is a summary of the two courses.

Digital Signal Processing (Theory and Application) - Tuesday 26th to Thursday 28th May 2026.

https://www.conted.ox.ac.uk/courses/digital-signal-processing-theory-and-application

This course provides a good understanding of DSP principles and their implementation and equips the delegate to put the ideas into practice and/or to tackle more advanced aspects of DSP. 'Hands-on' laboratory sessions are interspersed with the lectures to illustrate the taught material and allow you to pursue your own areas of interest in DSP. The hands-on sessions use specially written software running on PCs.

Subjects include:

  • Theoretical Foundations
  • Digital Filtering
  • Fourier Transforms And Frequency Domain Processing
  • DSP Hardware And Programming
  • ASIC Implementation
  • Typical DSP Applications

Digital Signal Processing Implementation (algorithms to optimization) - Friday 29th May 2026.

A one-day supplement to the Digital Signal Processing course that takes the theory and translates it into practice.

https://www.conted.ox.ac.uk/courses/digital-signal-processing-implementation-algorithms-to-optimisation

The course will include a mixed lecture and demonstration format and has been written to be independent of target processor architecture.

The course will show how to take common DSP algorithms and map them onto common processor architectures. It will also give a guide line for how to choose a DSP device, in particular how to choose and use the correct data word length for any application.

Attendee Feedback From Previous Courses:

John is like a textbook in human form ;-)  

It was informative, enjoyable and stimulating 

Excellent content, very lively thanks to the 2 excellent presenters - Anonymous

A very good introduction to DSP theory

Excellent lecturers! Really useful information and very understandable

Great mix of theory and practice

The lecturers gave a detailed and excellent explanation of the fundamental topics of DSP with real world engineering practice.

This session closes the gap and clears up much confusion between classroom DSP theories and actual DSP implementation.

Very good session, with in-depth discussion on the math and background.


These courses will be held at the University of Oxford, UK

Copyright © 2026 Delta Numerix


Friday, 28 November 2025

Using Large Language Models To Understand Digital Signal Processing Algorithms

I recently needed to understand some detailed timing error detectors for a modem.

I have the books and have implemented a TED before but this algorithm (Gardner's) was new to me to I turned to my current favourite LLM (Chat-GPT v5.1) to explain the algorithm and it saved me hours of reading books. 😃

If you are stuck understanding an algorithm then I highly recommend this approach as it gives a detailed text based introduction and generates and runs Python code to demonstrate it.

For example: "explain the FFT algorithm" it has provided a great introduction and written a complete FFT function and then compared the results to Numpy. 😃

In my case, it took 3 attempts to get the code for the TED to compile (it's a complex algorithm) but Chat-GPT automatically detected the errors and kept going until the code ran successfully. 👍🏻

Tuesday, 26 August 2025

Numerix-DSP Digital Signal Processing And Machine Learning Videos

Here is a selection of DSP and ML related videos presented by John Edwards

DSP Online Conference

2022 - Building A Tensorflow Lite Neural Network Vibration Classifier, With A Little Help From DSP

2021 - An Introduction To High Efficiency And Multi-rate Digital Filters

2020 - Frequency Domain Signal Processing


TinyML Foundation / Edge AI Foundation

2020 - “Low MIPS & Memory Machine Learning Industrial Vibration Monitoring Solution - AKA Not All AI Applications Are Cat v Dogs on Facebook ;-)


SigLib

SigLib DSP Library Introduction

SigLib Vibration Monitoring Machine Learning Demonstration


Data Science Festival 2020 - Lunch & Learn - "The Frequency Domain And How It Can Be Used To Aid Artificial Intelligence"


The 34th Annual Running Of The University Of Oxford Digital Signal Processing Course Will Be Held Online Again, In 2025

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

The course will run from Wednesday 22 Oct 2025 to Wednesday 26 Nov 2025, 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.

Copyright © 2025 Delta Numerix

Wednesday, 22 January 2025

The 34th Annual Running Of The University Of Oxford Live Digital Signal Processing Course, In May 2025

The 34th annual running of the University Of Oxford Live Digital Signal Processing course will be running in Oxford, UK, from Tuesday 20th to Friday 23rd May 2025.

The courses are presented by experts from industry for Engineers in industry and over the last 30 years has trained many hundreds of Engineers, from all areas of Science and Engineering.

Here is a summary of the two courses.

Digital Signal Processing (Theory and Application) - Tuesday 20th to Thursday 22nd May 2025.

https://www.conted.ox.ac.uk/courses/digital-signal-processing-theory-and-application

This course provides a good understanding of DSP principles and their implementation and equips the delegate to put the ideas into practice and/or to tackle more advanced aspects of DSP. 'Hands-on' laboratory sessions are interspersed with the lectures to illustrate the taught material and allow you to pursue your own areas of interest in DSP. The hands-on sessions use specially written software running on PCs.

Subjects include:

  • Theoretical Foundations
  • Digital Filtering
  • Fourier Transforms And Frequency Domain Processing
  • DSP Hardware And Programming
  • ASIC Implementation
  • Typical DSP Applications

Digital Signal Processing Implementation (algorithms to optimization) - Friday 23rd May 2025.

A one-day supplement to the Digital Signal Processing course that takes the theory and translates it into practice.

https://www.conted.ox.ac.uk/courses/digital-signal-processing-implementation-algorithms-to-optimisation

The course will include a mixed lecture and demonstration format and has been written to be independent of target processor architecture.

The course will show how to take common DSP algorithms and map them onto common processor architectures. It will also give a guide line for how to choose a DSP device, in particular how to choose and use the correct data word length for any application.

Attendee Feedback From Previous Courses:

John is like a textbook in human form ;-)  

It was informative, enjoyable and stimulating 

Excellent content, very lively thanks to the 2 excellent presenters - Anonymous

A very good introduction to DSP theory

Excellent lecturers! Really useful information and very understandable

Great mix of theory and practice

The lecturers gave a detailed and excellent explanation of the fundamental topics of DSP with real world engineering practice.

This session closes the gap and clears up much confusion between classroom DSP theories and actual DSP implementation.

Very good session, with in-depth discussion on the math and background.


These courses will be held at the University of Oxford, UK

Copyright © 2025 Delta Numerix


Wednesday, 15 January 2025

Understanding First Order Filters

While sorting through some very old papers I came across a solution to an interesting problem that I I struggled with when I was learning DSP. I have no idea where the original problem came from so I've replicated it here, as best I can remember, along with the solution:

The following first order direct form II filter :

                 w(n)
x(n) -->+-------------------+-->y(n)
        ^         |         ^
        |       +----+      |
        |       |z^-1|      |
        |       +----+      |
        |         |         |
        |         v         |
        ----*-----------*----
           a1  w(n-1)  b1

Is defined by the following equations:

y(n) = w(n) + b1.w(n-1)     (1)

w(n) = x(n) + a1.w(n-1)     (2)

Question: Show the difference equation in terms of y and x ?

Hint: Rearranging to a direct form I filter structure will help.

Solution

Diagramatically

The original system is a Linear Time Invariant (LTI) system so the feedforward and feedback sections can be swapped without changing the system response:

x(n) -------------+-------------->y(n)
         |        ^        |
       +----+     |      +----+
       |z^-1|     |      |z^-1|
       +----+     |      +----+
         |        |        |
         v        |        v
         ----*----+----*----
            b1        a1

Hence:

y(n) = x(n) + b1.x(n-1) + a1.y(n-1)


Mathematically

From (2):

w(n-1) = x(n-1) + a1.w(n-2)     (3)

Substituting (2) and (3) into (1), to compute the output:

y(n) = x(n) + a1.w(n-1) + b1.[x(n-1) + a1.(w(n-2)]     (4)

Rearranging to combine w terms:

y(n) = x(n) + b1.x(n-1) + a1.[w(n-1) + b1.w(n-2)]     (5)

From (1):     y(n-1) = w(n-1) + b1.w(n-2)     (6)

Substituting (6) into (5) gives:

y(n) = x(n) + b1 x(n-1) + a1 y(n-1)


Copyright © 2025 Delta Numerix

Tuesday, 19 November 2024

Using Generative AI And Large Language Models (LLMs) To Write DSP Code - Autumn 2024 Update

Having previously written a couple of blog posts regarding the use of LLMs to write DSP code, I've spent the last few months working on a project that has shown the landscape has changed dramatically.

The previous blog posts are here:

Using Generative AI And Large Language Models (LLMs) To Write DSP Code


At the start of the project, Claude 3.5 Sonnet was a step up from Chat-GPT 3 and Gemini really didn't cut the mustard at all. Then Chat-GPT 4o was released and this is a game changer, it has much more knowledge about the nuances of signal processing libraries such as scipy.signal and while it still struggles sometimes to write C code, I find it is by far the best option.

While I am lucky enough to have paid access to Chat-GPT 4o, not everyone does however there is an option through GitHub Marketplace that should work for most people. The nice thing about this is that you can easily try different LLMs but I find sticking to GPT 4o is the best option for me.

Copyright © 2024 Delta Numerix


Thursday, 25 July 2024

Why Mel-frequency Cepstrum Analysis Is Not Always The Ideal Solution For Vibration Analysis

The Mel-frequency Cepstrum (MFC) and it's associated outputs, the Mel-frequency Cepstral Coefficients (MFCCs), are commonly used for speech applications such as speaker and speech recognition, using neural networks. Unfortunately, the nature of the MFC means that it is not always ideally suited to applications such as vibration analysis and predictive maintenance

The MFC uses logarithmicaly spaced frequency banks to replicate how the human ear hears sound. This approach can lead to very large savings in the number of MIPS required for the recognition part of speaker and speech recognition. Unfortunately, this logarithmic frequency space hides frequencies that are closely spaced meaning that this approach is sub-optimal for applications such as machine vibration analysis, where small variations in vibrational frequency can indicate problems with the machine, particularly the bearings.

The following diagram shows a simple Mel-spaced filterbank, with 12 separate filters:


As can be seen from the diagram, resolution of close by frequencies is a particular problem for higher frequency harmonics, where the filters have a wider bandwidth.

The problem can also be seen in the following two images, which are sampled from identical machines running with two different error modes. It can be seen that it is the higher frequency peaks (1 kHz to 2 kHz) that vary the most and this is just the region, for this Mel-spaced filterbank, where the filter bandwidths start to get exessively wide.

Vibration Mode #1

Vibration Mode #2

The solution to this problem is to use a regular Fast Fourier Transform (FFT) for the front-end processing of these types of applications and use spectral analysis of the anticpated vibration modes to observe the frequency resolution required and this will then define the FFT size required for the application.

The SigLib Digital Signal Processing and Machine Learning library includes examples for machine vibration monitoring. These can be found here.

Copyright © 2024 Delta Numerix


Thursday, 27 June 2024

Using Generative AI And Large Language Models (LLMs) To Write DSP Code

Back in March 2023 I wrote the following blog post about using Generative AI and Large Language Models (LLMs) to write code: Are Chat-GPT and Google Bard The New Frontier For Writing DSP Code?

Since then, I have used these tools in many projects and have made a number of observations. In general, the more complex the task you are setting for the LLM, the more likely the performance of each is going to diverge and also the more likely it is that, as a programmer, you are going to have to test the code extensively to find the bugs.

Using these tools is a bit like an artist generating a preliminary sketch, rather than the final polished painting, with all of the correct detail.

I have found three main uses for Generative AI in coding:

  • Writing code to meet a specification
  • Documenting / commenting existing code
  • Converting code from one language / system to another

I have tried all of the following: Gemini, Google Code Assist, Chat-GPT, Bing and Co-Pilot. I have had the best coding results with Gemini (Bard) however if I find that this is struggling then I will try them all because they all have strengths and weaknesses.

It is important that you know what you want to do because there is no guarantee you will receive a correct answer! 

A useful trick is to try the same request multiple times because, unlike a traditional search engine, an LLM with give you a different responses each time. Handily, Gemini automatically generates 3 draft solutions and you can click on the tabs provide to review each.

I have observed that LLMs are much better at writing Python than lower level languages (C/C++ etc.). In Python, it will almost certainly produce a working solution using the Numpy/Scipy library functions, that may just need some final tuning.

If you are writing code for a lower level language then the best option is often to take a two stage approach:

  • Generate Python/Numpy/Scipy code
  • Convert the Python code to C - LLMs are very good at converting Numpy/Scipy functions to C

Generative AI is very good for converting between languages and Gemini will add comments to code that does not contain original comments. This is particularly useful if you work with a colleague who is not very dilligent with their code commenting ;-). It is worth noting, however, that comments are sometimes wrong due to AI misunderstanding the intention of the code.

Sometimes the conversion process will skip complex code sections, in a program, entirely so if this happens then the next step is to copy those sections and convert them separately.

Converting code from Python to C/C++ is generally very easy because they both use 0-based array indexing. Matlab, however, is more complex because it uses 1-based array indexing and this confuses the LLM. When converting Matlab code to Python or C/C+++ then I generally use the following request, which I then follow with the code section:

convert the following matlab code, with 1 based array indexing, to Python and Numpy, with 0 based array indexing

One final example of a gotcha is that Matlab uses FIR filter order whereas Scipy uses the number of coefficients.

As well as documenting code, LLMs are very good at debugging code however it is often important to explicitly specify the language in the request, rather than leaving it to the LLM to decide what language the code is written in.

Finally, Test! Test! Test!

Copyright © 2024 Delta Numerix


Wednesday, 8 May 2024

Digital Filter Plus Filter Design Tool Now Open Sourced

Back in the 1980s one of my first tasks, as a junior engineer, was to write a very simple filter design program, in GWBasic!

In the 1990s I updated it to Borland C and added some new functionality.

In the 2000s I added a GUI front end and more functionality.

For the last 20 years this has done myself and my customers well and has been my goto filter design tool but has languished in recent times so I've finally found the time to open-source it and it's now part of the SigLib DSP library.

https://github.com/Numerix-DSP/siglib

Over the last few months, at the request of customers, I've also added commercial grade audio Automatic Gain Control and multi-dimensional Kalman Filters to the library.

Enjoy :-)

SigLib and all of it's components, including Digital Filter Plus are licensed for free for educational and personal use only. All other uses require a developer's license.

Copyright © 2024 Delta Numerix

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()

Copyright © 2024 Delta Numerix


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.

Copyright © 2024 Delta Numerix

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.

Copyright © 2024 Delta Numerix


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:])

Copyright © 2023 Delta Numerix


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.