<commit_before>// 
// mskmodem_test.c
//

#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <math.h>
#include "liquid.h"

#define OUTPUT_FILENAME "mskmodem_test.m"

// print usage/help message
void usage()
{
    printf("mskmodem_test -- minimum-shift keying modem example\n");
    printf("options (default values in <>):\n");
    printf("  h     : print help\n");
    printf("  k     : samples/symbol,         default:  8\n");
    printf("  n     : number of data symbols, default: 10\n");
    printf("  s     : SNR [dB] <30>\n");
}

int main(int argc, char*argv[]) {
    // options
    unsigned int k=8;                   // filter samples/symbol
    unsigned int num_data_symbols = 20; // number of data symbols
    float SNRdB = 30.0f;                // signal-to-noise ratio [dB]

    int dopt;
    while ((dopt = getopt(argc,argv,"hk:n:s:")) != EOF) {
        switch (dopt) {
        case 'h': usage();                         return 0;
        case 'k': k = atoi(optarg);                break;
        case 'n': num_data_symbols = atoi(optarg); break;
        case 's': SNRdB = atof(optarg);            break;
        default:
            exit(1);
        }
    }

    unsigned int i;

    // derived values
    unsigned int num_symbols = num_data_symbols;
    unsigned int num_samples = k*num_symbols;
    float nstd = powf(10.0f, -SNRdB/20.0f);

    // arrays
    unsigned char sym_in[num_symbols];      // input symbols
    float phi[num_samples];                 // transmitted phase
    float complex x[num_samples];           // transmitted signal
    float complex y[num_samples];           // received signal
    float complex z[num_samples];           // output...
    unsigned char sym_out[num_symbols];     // output symbols

    // create transmit/receive interpolator/decimator
    enum {
        TXFILT_SQUARE=0,
        TXFILT_RCOS_FULL,
        TXFILT_RCOS_HALF,
        TXFILT_GMSK,
    } tx_filter_type = TXFILT_RCOS_HALF;

    float theta = 0.0f;
    unsigned int ht_len = 0;
    unsigned int tx_delay = 0;
    float * ht = NULL;
    switch (tx_filter_type) {
    case TXFILT_SQUARE:
        // regular MSK
        theta = M_PI / 4.0f;
        ht_len = k;
        tx_delay = 1;
        ht = (float*) malloc(ht_len *sizeof(float));
        for (i=0; i<ht_len; i++)
            ht[i] = M_PI / (2.0f * k);
        break;
    case TXFILT_RCOS_FULL:
        // full-response raised-cosine pulse
        theta = M_PI / 4.0f;
        ht_len = k;
        tx_delay = 1;
        ht = (float*) malloc(ht_len *sizeof(float));
        for (i=0; i<ht_len; i++)
            ht[i] = M_PI / (2.0f*k) * (1.0f - cosf(2.0f*M_PI*i/(float)ht_len));
        break;
    case TXFILT_RCOS_HALF:
        // partial-response raised-cosine pulse
        theta = M_PI / 4.0f;
        ht_len = 3*k;
        tx_delay = 2;
        ht = (float*) malloc(ht_len *sizeof(float));
        for (i=0; i<ht_len; i++)
            ht[i] = 0.0f;
        for (i=0; i<2*k; i++)
            ht[i+k/2] = 0.5f * M_PI / (2.0f*k) * (1.0f - cosf(2.0f*M_PI*i/(float)(2*k)));
        break;
    case TXFILT_GMSK:
        ht_len = 2*k*3+1;
        tx_delay = 3;
        ht = (float*) malloc(ht_len *sizeof(float));
        liquid_firdes_gmsktx(k,3,0.35f,0.0f,ht);
        for (i=0; i<ht_len; i++)
            ht[i] *= 1.0f / (float)k;
        break;
    default:
        fprintf(stderr,"error: %s, invalid tx filter type\n", argv[0]);
        exit(1);
    }
    for (i=0; i<ht_len; i++)
        printf("ht(%3u) = %12.8f;\n", i+1, ht[i]);
    interp_rrrf interp_tx = interp_rrrf_create(k, ht, ht_len);
    free(ht);

    // generate symbols and interpolate
    // phase-accumulating filter (trapezoidal integrator)
    float b[2] = {0.5f,  0.5f};
    if (tx_filter_type == TXFILT_SQUARE) {
        b[0] = 0.0f;
        b[1] = 1.0f;
    }
    float a[2] = {1.0f, -1.0f};
    iirfilt_rrrf integrator = iirfilt_rrrf_create(b,2,a,2);
    iirfilt_rrrf_execute(integrator, theta, &theta);
    for (i=0; i<num_symbols; i++) {
        sym_in[i] = rand() % 2;
        interp_rrrf_execute(interp_tx, sym_in[i] ? 1.0f : -1.0f, &phi[k*i]);

        // accumulate phase
        unsigned int j;
        for (j=0; j<k; j++) {
            iirfilt_rrrf_execute(integrator, phi[i*k+j], &theta);
            x[i*k+j] = cexpf(_Complex_I*theta);
        }
    }
    iirfilt_rrrf_destroy(integrator);

    // push through channel
    for (i=0; i<num_samples; i++)
        y[i] = x[i] + nstd*(randnf() + _Complex_I*randnf())*M_SQRT1_2;
    
    // create decimator
    unsigned int m = 3;
    float bw = 0.9f / (float)k;
    firfilt_crcf decim_rx = firfilt_crcf_create_kaiser(2*k*m+1, bw, 60.0f, 0.0f);
    printf("bw = %f\n", bw);

    // run receiver
    unsigned int n=0;
    float complex z_prime = 0.0f;
    for (i=0; i<num_samples; i++) {
        // push through filter
        firfilt_crcf_push(decim_rx, y[i]);
        firfilt_crcf_execute(decim_rx, &z[i]);

        z[i] *= 2.0f * bw;

        // decimate output
        if ( (i%k)==0 ) {
            float phi_hat = cargf(conjf(z_prime) * z[i]);
            z_prime = z[i];

            printf("%3u : %12.8f + j%12.8f (%1u)", n, crealf(z[i]), cimagf(z[i]), phi_hat > 0 ? 1 : 0);
            if (n >= m+tx_delay) printf(" (%1u)\n", sym_in[n-m-tx_delay]);
            else                 printf("\n");
            n++;
        }
    }

    // compute number of errors
    unsigned int num_errors = 0;
    for (i=0; i<num_data_symbols; i++)
        num_errors += sym_in[i] == sym_out[i] ? 0 : 1;

    printf("errors : %3u / %3u\n", num_errors, num_data_symbols);

    // destroy objects
    interp_rrrf_destroy(interp_tx);
    firfilt_crcf_destroy(decim_rx);

    // compute power spectral density
    unsigned int nfft = 1024;<commit_msg>sandbox/mskmodem_test: fixing GMSK pulse timing, counting actual errors<commit_after>// 
// mskmodem_test.c
//

#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <math.h>
#include "liquid.h"

#define OUTPUT_FILENAME "mskmodem_test.m"

// print usage/help message
void usage()
{
    printf("mskmodem_test -- minimum-shift keying modem example\n");
    printf("options (default values in <>):\n");
    printf("  h     : print help\n");
    printf("  k     : samples/symbol,         default:  8\n");
    printf("  n     : number of data symbols, default: 10\n");
    printf("  s     : SNR [dB] <30>\n");
}

int main(int argc, char*argv[]) {
    // options
    unsigned int k=8;                   // filter samples/symbol
    unsigned int num_data_symbols = 20; // number of data symbols
    float SNRdB = 30.0f;                // signal-to-noise ratio [dB]

    int dopt;
    while ((dopt = getopt(argc,argv,"hk:n:s:")) != EOF) {
        switch (dopt) {
        case 'h': usage();                         return 0;
        case 'k': k = atoi(optarg);                break;
        case 'n': num_data_symbols = atoi(optarg); break;
        case 's': SNRdB = atof(optarg);            break;
        default:
            exit(1);
        }
    }

    unsigned int i;

    // derived values
    unsigned int num_symbols = num_data_symbols;
    unsigned int num_samples = k*num_symbols;
    float nstd = powf(10.0f, -SNRdB/20.0f);

    // arrays
    unsigned char sym_in[num_symbols];      // input symbols
    float phi[num_samples];                 // transmitted phase
    float complex x[num_samples];           // transmitted signal
    float complex y[num_samples];           // received signal
    float complex z[num_samples];           // output...
    //unsigned char sym_out[num_symbols];     // output symbols

    // create transmit/receive interpolator/decimator
    enum {
        TXFILT_SQUARE=0,
        TXFILT_RCOS_FULL,
        TXFILT_RCOS_HALF,
        TXFILT_GMSK,
    } tx_filter_type = TXFILT_GMSK;

    float theta = 0.0f;
    unsigned int ht_len = 0;
    unsigned int tx_delay = 0;
    float * ht = NULL;
    switch (tx_filter_type) {
    case TXFILT_SQUARE:
        // regular MSK
        theta = M_PI / 4.0f;
        ht_len = k;
        tx_delay = 1;
        ht = (float*) malloc(ht_len *sizeof(float));
        for (i=0; i<ht_len; i++)
            ht[i] = M_PI / (2.0f * k);
        break;
    case TXFILT_RCOS_FULL:
        // full-response raised-cosine pulse
        theta = M_PI / 4.0f;
        ht_len = k;
        tx_delay = 1;
        ht = (float*) malloc(ht_len *sizeof(float));
        for (i=0; i<ht_len; i++)
            ht[i] = M_PI / (2.0f*k) * (1.0f - cosf(2.0f*M_PI*i/(float)ht_len));
        break;
    case TXFILT_RCOS_HALF:
        // partial-response raised-cosine pulse
        theta = M_PI / 4.0f;
        ht_len = 3*k;
        tx_delay = 2;
        ht = (float*) malloc(ht_len *sizeof(float));
        for (i=0; i<ht_len; i++)
            ht[i] = 0.0f;
        for (i=0; i<2*k; i++)
            ht[i+k/2] = 0.5f * M_PI / (2.0f*k) * (1.0f - cosf(2.0f*M_PI*i/(float)(2*k)));
        break;
    case TXFILT_GMSK:
        theta = M_PI / 4.0f;
        ht_len = 2*k*3+1+k;
        tx_delay = 4;
        ht = (float*) malloc(ht_len *sizeof(float));
        for (i=0; i<ht_len; i++)
            ht[i] = 0.0f;
        liquid_firdes_gmsktx(k,3,0.9f,0.0f,&ht[k/2]);
        for (i=0; i<ht_len; i++)
            ht[i] *= 1.0f / (float)k;
        break;
    default:
        fprintf(stderr,"error: %s, invalid tx filter type\n", argv[0]);
        exit(1);
    }
    for (i=0; i<ht_len; i++)
        printf("ht(%3u) = %12.8f;\n", i+1, ht[i]);
    interp_rrrf interp_tx = interp_rrrf_create(k, ht, ht_len);
    free(ht);

    // generate symbols and interpolate
    // phase-accumulating filter (trapezoidal integrator)
    float b[2] = {0.5f,  0.5f};
    if (tx_filter_type == TXFILT_SQUARE) {
        b[0] = 0.0f;
        b[1] = 1.0f;
    }
    float a[2] = {1.0f, -1.0f};
    iirfilt_rrrf integrator = iirfilt_rrrf_create(b,2,a,2);
    iirfilt_rrrf_execute(integrator, theta, &theta);
    for (i=0; i<num_symbols; i++) {
        sym_in[i] = rand() % 2;
        interp_rrrf_execute(interp_tx, sym_in[i] ? 1.0f : -1.0f, &phi[k*i]);

        // accumulate phase
        unsigned int j;
        for (j=0; j<k; j++) {
            iirfilt_rrrf_execute(integrator, phi[i*k+j], &theta);
            x[i*k+j] = cexpf(_Complex_I*theta);
        }
    }
    iirfilt_rrrf_destroy(integrator);

    // push through channel
    for (i=0; i<num_samples; i++)
        y[i] = x[i] + nstd*(randnf() + _Complex_I*randnf())*M_SQRT1_2;
    
    // create decimator
    unsigned int m = 3;
    float bw = 0.9f / (float)k;
    firfilt_crcf decim_rx = firfilt_crcf_create_kaiser(2*k*m+1, bw, 60.0f, 0.0f);
    printf("bw = %f\n", bw);

    // run receiver
    unsigned int n=0;
    unsigned int num_errors = 0;
    unsigned int num_symbols_checked = 0;
    float complex z_prime = 0.0f;
    for (i=0; i<num_samples; i++) {
        // push through filter
        firfilt_crcf_push(decim_rx, y[i]);
        firfilt_crcf_execute(decim_rx, &z[i]);

        z[i] *= 2.0f * bw;

        // decimate output
        if ( (i%k)==0 ) {
            float phi_hat = cargf(conjf(z_prime) * z[i]);
            unsigned int sym_out = phi_hat > 0 ? 1 : 0; // estimated transmitted symbol
            z_prime = z[i];

            printf("%3u : %12.8f + j%12.8f (%1u)", n, crealf(z[i]), cimagf(z[i]), sym_out);
            if (n >= m+tx_delay) {
                num_errors += (sym_out == sym_in[n-m-tx_delay]) ? 0 : 1;
                num_symbols_checked++;
                printf(" (%1u)\n", sym_in[n-m-tx_delay]);
            } else {
                printf("\n");
            }
            n++;
        }
    }

    // print number of errors
    printf("errors : %3u / %3u\n", num_errors, num_symbols_checked);

    // destroy objects
    interp_rrrf_destroy(interp_tx);
    firfilt_crcf_destroy(decim_rx);

    // compute power spectral density
    unsigned int nfft = 1024;
{'fraction_non_alphanumeric': 0.15096251266464034, 'fraction_numerical': 0.02524484971293482, 'mean_word_length': 1.972396486825596, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6055178', 'n_tokens_mistral': 4784, 'n_tokens_neox': 4536, 'n_words': 975}