[PATCH i-g-t 3/3] lib/igt_audio: Replace GSL FFT usage with meow_fft

[Peter Senna Tschudin]

On 13.09.2024 12:44, Kamil Konieczny wrote:
> From: Ryszard Knop <ryszard.knop at intel.com>
> 
> Tested by running kms_chamelium tests related with audio and comparing
> GSL and meow_fft outputs in a separate test program. It appears meow_fft
> is slightly less accurate than GSL (result differs after the 6th decimal
> place), but the error is low enough that it does not matter here.
> 
> v2: sort headers, remove blank spaces, fixed codestyle (Kamil)
> 

Reviewed-by: Peter Senna Tschudin <peter.senna at linux.intel.com>
> Signed-off-by: Ryszard Knop <ryszard.knop at intel.com>
> Signed-off-by: Kamil Konieczny <kamil.konieczny at linux.intel.com>
> ---
>  lib/igt_audio.c | 74 +++++++++++++++++++++++++++++++++++++++----------
>  1 file changed, 60 insertions(+), 14 deletions(-)
> 
> diff --git a/lib/igt_audio.c b/lib/igt_audio.c
> index e0b1bafe1..193f9aa05 100644
> --- a/lib/igt_audio.c
> +++ b/lib/igt_audio.c
> @@ -28,12 +28,14 @@
>  
>  #include <errno.h>
>  #include <fcntl.h>
> -#include <gsl/gsl_fft_real.h>
> +#include <limits.h>
>  #include <math.h>
>  #include <unistd.h>
>  
>  #include "igt_audio.h"
> +#include "igt_aux.h"
>  #include "igt_core.h"
> +#include "meow_fft/meow_fft.h"
>  
>  #define FREQS_MAX 64
>  #define CHANNELS_MAX 8
> @@ -322,6 +324,53 @@ static double hann_window(double v, size_t i, size_t N)
>  	return v * 0.5 * (1 - cos(2.0 * M_PI * (double) i / (double) N));
>  }
>  
> +/**
> + * run_fft:
> + * @array: The signal to run FFT on
> + * @length: The signal buffer length (must be a power of 2)
> + *
> + * Run the Fast Fourier Transform on the provided signal, whose
> + * length must be a power of 2. The return value points to an
> + * array of N/2 structs with real (.r) and imaginary (.j) parts.
> + *
> + * For the 0-th FFT element, only the real part is valid - its
> + * imaginary part is always 0 and is not saved anywhere.
> + *
> + * For the (N/2)-th element, again, only the real part is valid.
> + * The array does not have enough space to save it as a separate
> + * element, so its real part is saved in the 0-th element's
> + * imaginary part.
> + */
> +static Meow_FFT_Complex *run_fft(double *array, size_t length)
> +{
> +	size_t workset_bytes;
> +	Meow_FFT_Complex *fft_data;
> +	struct Meow_FFT_Workset_Real *fft_workset;
> +
> +	igt_assert(length >= 2 && is_power_of_two(length));
> +
> +	// Get size for a N point fft working on non-complex (real) data.
> +	workset_bytes = meow_fft_generate_workset_real(length, NULL);
> +	if (workset_bytes == 0)
> +		return NULL;
> +
> +	fft_data = malloc(sizeof(Meow_FFT_Complex) * length);
> +	if (!fft_data)
> +		return NULL;
> +
> +	fft_workset = (struct Meow_FFT_Workset_Real *)malloc(workset_bytes);
> +	if (!fft_workset) {
> +		free(fft_data);
> +		return NULL;
> +	}
> +
> +	meow_fft_generate_workset_real(length, fft_workset);
> +	meow_fft_real(fft_workset, array, fft_data);
> +	free(fft_workset);
> +
> +	return fft_data;
> +}
> +
>  /**
>   * Checks that frequencies specified in signal, and only those, are included
>   * in the input data.
> @@ -333,11 +382,12 @@ bool audio_signal_detect(struct audio_signal *signal, int sampling_rate,
>  			 int channel, const double *samples, size_t samples_len)
>  {
>  	double *data;
> +	Meow_FFT_Complex *fft_data;
>  	size_t data_len = samples_len;
>  	size_t bin_power_len = data_len / 2 + 1;
>  	double bin_power[bin_power_len];
>  	bool detected[FREQS_MAX];
> -	int ret, freq_accuracy, freq, local_max_freq;
> +	int freq_accuracy, freq, local_max_freq;
>  	double max, local_max, threshold;
>  	size_t i, j;
>  	bool above, success;
> @@ -360,27 +410,22 @@ bool audio_signal_detect(struct audio_signal *signal, int sampling_rate,
>  	freq_accuracy = sampling_rate / data_len;
>  	igt_debug("Allowed freq. error: %d Hz\n", freq_accuracy);
>  
> -	ret = gsl_fft_real_radix2_transform(data, 1, data_len);
> -	if (ret != 0) {
> +	fft_data = run_fft(data, data_len);
> +	if (!fft_data) {
>  		free(data);
>  		igt_assert(0);
>  	}
>  
> -	/* Compute the power received by every bin of the FFT.
> -	 *
> -	 * For i < data_len / 2, the real part of the i-th term is stored at
> -	 * data[i] and its imaginary part is stored at data[data_len - i].
> -	 * i = 0 and i = data_len / 2 are special cases, they are purely real
> -	 * so their imaginary part isn't stored.
> -	 *
> +	/* Compute the power received by every bin of the FFT. The run_fft
> +	 * function docs explain the special cases outside of the for loop.
>  	 * The power is encoded as the magnitude of the complex number and the
>  	 * phase is encoded as its angle.
>  	 */
> -	bin_power[0] = data[0];
> +	bin_power[0] = fft_data[0].r;
>  	for (i = 1; i < bin_power_len - 1; i++) {
> -		bin_power[i] = hypot(data[i], data[data_len - i]);
> +		bin_power[i] = hypot(fft_data[i].r, fft_data[i].j);
>  	}
> -	bin_power[bin_power_len - 1] = data[data_len / 2];
> +	bin_power[bin_power_len - 1] = fft_data[0].j;
>  
>  	/* Normalize the power */
>  	for (i = 0; i < bin_power_len; i++)
> @@ -487,6 +532,7 @@ bool audio_signal_detect(struct audio_signal *signal, int sampling_rate,
>  		}
>  	}
>  
> +	free(fft_data);
>  	free(data);
>  
>  	return success;