Monday, 20 March 2023

Is Chat-GPT 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.

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


Conclusion

For simple functions such as FIR and IIR filters, the code generated does exactly what was asked however the code is not really useable for embedded DSP applications. The main reasons are that when asking simple requests, Chat-GPT doesn'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, 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. And similarly compare the results when Google Bard is published.

Note: I understand that Bing uses Chat-GPT V4.0. I'll give it a whirl and let you know.

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.