QAM modulation in C without complex library

demodulation.c

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

typedef struct _complex {
    double imag;
    double real;
} complex_t;

const int SAMPLE_RATE = 160E3;
const float FREQUENCY = 40E3;
const int N_OSCILATIONS = 4;
const int N_SAMPLES_PER_SYMBOL = (N_OSCILATIONS * SAMPLE_RATE) / FREQUENCY;
const float OMEGA_PART = -2 * M_PI * FREQUENCY / SAMPLE_RATE;
const int SIGNAL_LENGTH = 2000;

int main() {
    int signal[SIGNAL_LENGTH] = ...;

    // The filter (`filterx`) is calculated only once, so you can store it
    complex_t filterx[N_SAMPLES_PER_SYMBOL];
    for (int i = 0; i < N_SAMPLES_PER_SYMBOL; i++) {
        filterx[i].real = cos(OMEGA_PART * i) / (N_SAMPLES_PER_SYMBOL / 2);
        filterx[i].imag = sin(OMEGA_PART * i) / (N_SAMPLES_PER_SYMBOL / 2);
    }

    int symbol_index = 0;
    for (int i = 0; i < SIGNAL_LENGTH - N_SAMPLES_PER_SYMBOL + 1; i += N_SAMPLES_PER_SYMBOL) {
        complex_t sum = {
            .real = 0,
            .imag = 0
        };

        // This part is computationaly the most expensive, however it can be accelerated
        // by utilising built in DSP convolution acceleration. Note that in case of using the convolution
        // `filterx` has to be inversed.
        for (int j = 0; j < N_SAMPLES_PER_SYMBOL; j++) {
            sum.real += filterx[j].real * signal[i + j];
            sum.imag += filterx[j].imag * signal[i + j];
        }

        // Here you will probably need some PLL and equalisation to get rid of real-world distortions
        // Please check fantastically explained algorithms by Joseph D. Gaeddert:
        // - PLL: https://liquidsdr.org/blog/pll-howto/
        // - LMS Equalizer: https://liquidsdr.org/blog/lms-equalizer/
        int symbol = ((sum.imag > 0) << 1) | (sum.real < 0);
        printf("%d ", symbol);
    }
    printf("\n");
}

modulation.c

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

typedef struct _complex {
    double imag;
    double real;
} complex_t;

const double SYMBOL_DURATION = 50E-6;
const int SAMPLE_RATE = 100E3;
const int FREQUENCY = 40E3;

const float SAMPLE_DURATION = 1.0 / SAMPLE_RATE;
const int N_SAMPLES_PER_SYMBOL = SYMBOL_DURATION / SAMPLE_DURATION;

int main() {
    complex_t symbol = {
        .real = 1,
        .imag = 0
    };
    complex_t carrier[N_SAMPLES_PER_SYMBOL];
    double signal[N_SAMPLES_PER_SYMBOL];

    // Generate carrier
    for (size_t i = 0; i < N_SAMPLES_PER_SYMBOL; i++) {
        double omega = 2 * M_PI * FREQUENCY * i * SAMPLE_DURATION;
        carrier[i].real = cos(omega);
        carrier[i].imag = sin(omega);
    }
    
    // Generate signal. We can ignore imaginary part
    for (size_t i = 0; i < N_SAMPLES_PER_SYMBOL; i++) {
        signal[i] = carrier[i].real * symbol.real - carrier[i].imag * symbol.imag;
    }

    for (size_t i = 0; i < N_SAMPLES_PER_SYMBOL; i++) {
        printf("%f ", signal[i]);
    }
    printf("\n");
}

modulation_complex.c

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

const double SYMBOL_DURATION = 50E-6;
const int SAMPLE_RATE = 100E3;
const int FREQUENCY = 40E3;

const float SAMPLE_DURATION = 1.0 / SAMPLE_RATE;
const int N_SAMPLES_PER_SYMBOL = SYMBOL_DURATION / SAMPLE_DURATION;

int main() {
    double carrier[N_SAMPLES_PER_SYMBOL];
    double signal[N_SAMPLES_PER_SYMBOL];

    for (size_t i = 0; i < N_SAMPLES_PER_SYMBOL; i++) {
        carrier[i] = cpow(M_E, (2*I) * M_PI * FREQUENCY * i * SAMPLE_DURATION);
    }

    double complex symbol = 1 + 0 * I;
    for (size_t i = 0; i < N_SAMPLES_PER_SYMBOL; i++) {
        signal[i] = carrier[i] * symbol;
    }

    for (size_t i = 0; i < N_SAMPLES_PER_SYMBOL; i++) {
        printf("%f ", creal(signal[i]));
    }
    printf("\n");
}