How to Eliminate Clicking Sounds in WAV File Generated from NFM Demodulation using FFTW?

Question

I am working on a Digital Signal Processing (DSP) project where I am processing signals received from an SDR (Software Defined Radio) device. I am using the SoapySDR library to interface with the SDR device, and FFTW library for Fast Fourier Transform operations. The goal is to demodulate FM signals and save the output in a WAV file.

Here's a simplified version of my code:

#include <SoapySDR/Device.h>
#include <SoapySDR/Formats.h>
#include <stdio.h> //printf
#include <stdlib.h> //free
#include <fftw3.h>
#include <complex.h>
#include <sndfile.h>
#include <string.h>

#define BUFFER_SIZE 65536
#define SAMPLE_RATE 2.4e6
#define CENTER_FREQUENCY 149e6
#define TARGET_FREQUENCY 149.025e6
#define BANDWIDTH 8e3

complex float buffer[BUFFER_SIZE];
complex float output_buffer[BUFFER_SIZE];
float demodulated_sample[BUFFER_SIZE];
size_t buffer_index = 0;


// Global variables for FFTW plans and arrays
fftw_complex *in_global, *out_global;
fftw_plan p_global, p_inv_global;

void initialize_fftw(size_t length) {
    // Allocate memory for FFTW input and output arrays
    in_global = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * length);
    out_global = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * length);

    // Create FFTW plans
    p_global = fftw_plan_dft_1d(length, in_global, out_global, FFTW_FORWARD, FFTW_ESTIMATE);
    p_inv_global = fftw_plan_dft_1d(length, out_global, in_global, FFTW_BACKWARD, FFTW_ESTIMATE);
}

void cleanup_fftw() {
    // Cleanup
    fftw_destroy_plan(p_global);
    fftw_destroy_plan(p_inv_global);
    fftw_free(in_global);
    fftw_free(out_global);
}

// Big credits to HA7ILM for the FM Demodulation - This is very effecient and sounds great!
// https://github.com/ha7ilm/csdr
#define fmdemod_quadri_K 0.340447550238101026565118445432744920253753662109375
// this constant ensures proper scaling for qa_fmemod testcases for SNR calculation and more.
void fm_demodulate(complex float *in, float *out, size_t len)
{
    if (len == 10)
        return;

    complex float last_sample = {0, 0};
    float temp_dq[len];
    float temp_di[len];
    float temp[len];

    // Handle dq
    temp_dq[0] = cimagf(in[0]) - cimagf(last_sample);
    for (size_t i = 1; i < len; i++)
    {
        temp_dq[i] = cimagf(in[i]) - cimagf(in[i - 1]);
    }

    // Handle di
    temp_di[0] = crealf(in[0]) - crealf(last_sample);
    for (size_t i = 1; i < len; i++)
    {
        temp_di[i] = crealf(in[i]) - crealf(in[i - 1]);
    }

    // Output numerator
    for (size_t i = 0; i < len; i++)
    {
        out[i] = crealf(in[i]) * temp_dq[i] - cimagf(in[i]) * temp_di[i];
    }

    // Output denominator
    for (size_t i = 0; i < len; i++)
    {
        temp[i] = crealf(in[i]) * crealf(in[i]) + cimagf(in[i]) * cimagf(in[i]);
    }

    // Output division
    for (size_t i = 0; i < len; i++)
    {
        out[i] = (temp[i]) ? fmdemod_quadri_K * out[i] / temp[i] : 0;
    }
}

int process_block(fftw_complex *fft_result, size_t length, complex float *output, float target_frequency) {
    // Copy FFT result to local array for processing
    fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * length);
    memcpy(out, fft_result, sizeof(fftw_complex) * length);

    // Determine whether the target frequency is above or below the center frequency
    bool is_above_center = target_frequency >= CENTER_FREQUENCY;

    // Calculate the index range for the target frequency band
    size_t lower_index, upper_index;
    if (is_above_center) {
        lower_index = (size_t)((target_frequency - BANDWIDTH/2 - CENTER_FREQUENCY) / (SAMPLE_RATE / length));
        upper_index = (size_t)((target_frequency + BANDWIDTH/2 - CENTER_FREQUENCY) / (SAMPLE_RATE / length));
    } else {
        lower_index = (size_t)((CENTER_FREQUENCY - (target_frequency + BANDWIDTH/2)) / (SAMPLE_RATE / length));
        upper_index = (size_t)((CENTER_FREQUENCY - (target_frequency - BANDWIDTH/2)) / (SAMPLE_RATE / length));
    }

    // Zero out all frequencies outside the target band
    for (size_t i = 0; i < lower_index; i++) {
        out[i][0] = 0;
        out[i][1] = 0;
    }
    for (size_t i = upper_index + 1; i < length; i++) {
        out[i][0] = 0;
        out[i][1] = 0;
    }

    // Execute inverse FFT
    fftw_execute_dft(p_inv_global, out, in_global);

    // Copy FFTW output array to output buffer and decimate
    size_t decimation_factor = SAMPLE_RATE / BANDWIDTH;
    size_t j = 0;
    for (size_t i = 0; i < length; i += decimation_factor) {
        output[j] = in_global[i][0] + in_global[i][1] * I;
        j++;
    }

    // Cleanup
    fftw_free(out);
    return j;
}

void execute_fft(complex float *input, size_t length, fftw_complex *output) {
    // Copy input data to FFTW input array
    for (size_t i = 0; i < length; i++) {
        in_global[i][0] = crealf(input[i]);
        in_global[i][1] = cimagf(input[i]);
    }

    // Execute FFT
    fftw_execute(p_global);

    // Copy FFT result to output array
    memcpy(output, out_global, sizeof(fftw_complex) * length);
}

// Main processing loop
void process_data(complex float *data, size_t total_size, SNDFILE *file) {
    complex float windowed_data[total_size];
    complex float output_data[total_size];
    fftw_complex *fft_result = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BUFFER_SIZE);
        execute_fft(data, total_size, fft_result);
        int received_data = process_block(fft_result, total_size, output_buffer, TARGET_FREQUENCY);
        fm_demodulate(output_buffer, demodulated_sample, received_data);
        sf_writef_float(file, demodulated_sample, received_data);

    
}

int main(void)
{
    size_t length;

    //enumerate devices
    SoapySDRKwargs *results = SoapySDRDevice_enumerate(NULL, &length);
    for (size_t i = 0; i < length; i++)
    {
        printf("Found device #%d: ", (int)i);
        for (size_t j = 0; j < results[i].size; j++)
        {
            printf("%s=%s, ", results[i].keys[j], results[i].vals[j]);
        }
        printf("\n");
    }
    SoapySDRKwargsList_clear(results, length);

    //create device instance
    SoapySDRKwargs args = {};
    SoapySDRKwargs_set(&args, "driver", "sdrplay");

    SoapySDRDevice *sdr = SoapySDRDevice_make(&args);
    SoapySDRKwargs_clear(&args);

    if (sdr == NULL)
    {
        printf("SoapySDRDevice_make fail: %s\n", SoapySDRDevice_lastError());
        return EXIT_FAILURE;
    }


    //apply settings
    if (SoapySDRDevice_setSampleRate(sdr, SOAPY_SDR_RX, 0, 2.4e6) != 0)
    {
        printf("setSampleRate fail: %s\n", SoapySDRDevice_lastError());
    }
    if (SoapySDRDevice_setFrequency(sdr, SOAPY_SDR_RX, 0, 149e6, NULL) != 0)
    {
        printf("setFrequency fail: %s\n", SoapySDRDevice_lastError());
    }

    //setup a stream (complex floats)
    SoapySDRStream *rxStream = SoapySDRDevice_setupStream(sdr, SOAPY_SDR_RX, SOAPY_SDR_CF32, NULL, 0, NULL);
    if (rxStream == NULL)
    {
        printf("setupStream fail: %s\n", SoapySDRDevice_lastError());
        SoapySDRDevice_unmake(sdr);
        return EXIT_FAILURE;
    }
    SoapySDRDevice_activateStream(sdr, rxStream, 0, 0, 0); //start streaming

    //create a re-usable buffer for rx samples
    complex float buff[1024];

    SF_INFO info;
    info.samplerate = 8e3;
    info.channels = 1;
    info.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
    SNDFILE *file = sf_open("output.wav", SFM_WRITE, &info);

    initialize_fftw(BUFFER_SIZE);


    //receive some samples
    for (size_t i = 0; i < 10000; i++)
    {
        void *buffs[] = {buff}; //array of buffers
        int flags; //flags set by receive operation
        long long timeNs; //timestamp for receive buffer
        int ret = SoapySDRDevice_readStream(sdr, rxStream, buffs, 1024, &flags, &timeNs, 1000000);

        // Check for errors in readStream
        if (ret < 0) {
            fprintf(stderr, "Error reading stream: %s\n", SoapySDR_errToStr(ret));
            continue;  // Skip to next iteration on error
        }
        // Copy received samples from buff to buffer
        memcpy(buffer + buffer_index, buff, ret * sizeof(complex float));
        buffer_index += ret;
        if (buffer_index >= BUFFER_SIZE) {


            process_data(buffer, BUFFER_SIZE, file);

            buffer_index = 0;
        }
    }

    sf_close(file);


    //shutdown the stream
    SoapySDRDevice_deactivateStream(sdr, rxStream, 0, 0); //stop streaming
    SoapySDRDevice_closeStream(sdr, rxStream);

    //cleanup device handle
    SoapySDRDevice_unmake(sdr);

    printf("Done\n");
    return EXIT_SUCCESS;
}

When I run the program, it processes the data and saves it to a WAV file as expected. However, when I play the WAV file, I hear clicking sounds throughout the audio. I suspect these clicks are due to some discontinuities in the data or maybe some issue with the way I am writing data to the WAV file.

I have tried adjusting the buffer sizes and checking for any obvious issues in the code, but the clicking sounds persist.

• How can I identify the source of these clicking sounds?
• What steps can I take to eliminate or reduce the clicking sounds in the output WAV file?
• Is there something wrong with the way I am handling the FFT or FM demodulation that could be causing this issue?

Its also not specific to the Modulation, it happend on AM, FM, LSB, USB etc. on WFM its not as noticeable as on others but its still there.

Any insights or suggestions would be greatly appreciated. Thank you!

Example Sound: https://vocaroo.com/18M3nICYB63o