diff --git a/dsd/CMakeLists.txt b/dsd/CMakeLists.txt index 9fd50a33f..b2d9d08fb 100644 --- a/dsd/CMakeLists.txt +++ b/dsd/CMakeLists.txt @@ -16,6 +16,7 @@ set(dsd_SOURCES dsd_mbe.c dsd_nocarrier.c dsd_opts.c + dsd_state.c dsd_symbol.c dsd_upsample.c dstar_const.c @@ -28,6 +29,7 @@ set(dsd_SOURCES nxdn96_const.c p25_lcw.c p25p1_const.c + p25p1_heuristics.c p25p1_hdu.c p25p1_ldu1.c p25p1_ldu2.c diff --git a/dsd/dsd_livescanner.c b/dsd/dsd_livescanner.c index 0b0490b98..72092280f 100644 --- a/dsd/dsd_livescanner.c +++ b/dsd/dsd_livescanner.c @@ -21,15 +21,7 @@ void liveScanner(dsd_opts * opts, dsd_state * state) { - if (opts->audio_in_fd == -1) - { - if (pthread_mutex_lock(&state->input_mutex)) - { - printf("liveScanner -> Unable to lock mutex\n"); - } - } - - while (state->dsd_running) + while (1) // FIXME: loop while input available { noCarrier(opts, state); state->synctype = getFrameSync(opts, state); @@ -39,7 +31,7 @@ void liveScanner(dsd_opts * opts, dsd_state * state) state->umid = (((state->max) - state->center) * 5 / 8) + state->center; state->lmid = (((state->min) - state->center) * 5 / 8) + state->center; - while ((state->synctype != -1) && (state->dsd_running)) + while (state->synctype != -1) // TODO: make sure it ends { processFrame(opts, state); state->synctype = getFrameSync(opts, state); diff --git a/dsd/dsd_state.h b/dsd/dsd_state.h index 5b2477e7b..807fd71c1 100644 --- a/dsd/dsd_state.h +++ b/dsd/dsd_state.h @@ -18,7 +18,6 @@ #ifndef INCLUDE_DSD_STATE_H_ #define INCLUDE_DSD_STATE_H_ -#include #include #include "p25p1_heuristics.h" @@ -104,22 +103,16 @@ typedef struct int exitflag; // the former global that cannot be a global within SDRangel and is not of much use with it anyway // New from original DSD for in-memory processing support with SDRangel: - pthread_mutex_t input_mutex; //!< mutex to communicate with input thread - pthread_cond_t input_ready; //!< signals that input demodulator samples are available for processing const short *input_samples; //!< demodulator samples int input_length; //!< 0: data not ready, >0: data ready for this amount of demodulator samples int input_offset; //!< consumer pointer - pthread_mutex_t output_mutex; //!< mutex to communicate with output (audio) thread - pthread_cond_t output_ready; //!< signals that output audio samples are ready short *output_buffer; //!< Output of decoder single S16LE int output_offset; //!< producer pointer short *output_samples; //!< L+R channels S16LE ready for writing to audio FIFO int output_num_samples; //!< Number of L+R samples available in the above buffer int output_length; //!< L+R buffer size (fixed) int output_finished; //!< 0: not ready, 1: ready - - int dsd_running; //!< True while DSD thread is running } dsd_state; #ifdef __cplusplus diff --git a/dsd/dsd_symbol.c b/dsd/dsd_symbol.c index 15a144719..c1f82a3d9 100644 --- a/dsd/dsd_symbol.c +++ b/dsd/dsd_symbol.c @@ -86,16 +86,8 @@ int getSymbol(dsd_opts * opts, dsd_state * state, int have_sync) { if (opts->audio_in_fd == -1) { - while (state->input_length == 0) - { - // If the buffer is empty, wait for more samples to arrive. - if (pthread_cond_wait(&state->input_ready, &state->input_mutex)) - { - printf("getSymbol -> Error waiting for input condition\n"); - } - } // Get the next sample from the buffer - sample = state->input_samples[state->input_offset++]; + sample = state->input_samples[state->input_offset++]; // FIXME: get sample only if available if (state->input_offset == state->input_length) // all available samples have been read { @@ -104,12 +96,6 @@ int getSymbol(dsd_opts * opts, dsd_state * state, int have_sync) // We've reached the end of the buffer. Wait for more next time. state->input_length = 0; - // make output samples availabele - if (pthread_mutex_lock(&state->output_mutex)) - { - printf("Unable to lock output mutex\n"); - } - state->output_num_samples = state->output_offset; // GNUradio drivel @@ -148,17 +134,6 @@ int getSymbol(dsd_opts * opts, dsd_state * state, int have_sync) } state->output_finished = 1; - - // Wake up audio out thread - if (pthread_cond_signal(&state->output_ready)) - { - printf("Unable to signal output ready\n"); - } - - if (pthread_mutex_unlock(&state->output_mutex)) - { - printf("Unable to unlock output mutex\n"); - } } } else diff --git a/dsd/p25p1_heuristics.c b/dsd/p25p1_heuristics.c new file mode 100644 index 000000000..f8fa4aaf6 --- /dev/null +++ b/dsd/p25p1_heuristics.c @@ -0,0 +1,386 @@ + +#define _USE_MATH_DEFINES +#include + +#include "dsd.h" + +#include "p25p1_heuristics.h" + +/** + * This module is dedicated to improve the accuracy of the digitizer. The digitizer is the piece of code that + * translates an analog value to an actual symbol, in the case of P25, a dibit. + * It implements a simple Gaussian classifier. It's based in the assumption that the analog values from the + * signal follow normal distributions, one single distribution for each symbol. + * Analog values for the dibit "0" will fit into a Gaussian bell curve with a characteristic mean and std + * distribution and the same goes for dibits "1", "2" and "3." + * Hopefully those bell curves are well separated from each other so we can accurately discriminate dibits. + * If we could model the Gaussian of each dibit, then given an analog value, the dibit whose Gaussian fits + * better is the most likely interpretation for that value. By better fit we can calculate the PDF + * (probability density function) for the Gaussian, the one with the highest value is the best fit. + * + * The approach followed here to model the Gaussian for each dibit is to use the error corrected information + * as precise oracles. P25 uses strong error correction, some dibits are doubly protected by Hamming/Golay + * and powerful Reed-Solomon codes. If a sequence of dibits clears the last Reed-Solomon check, we can be + * quite confident that those values are correct. We can use the analog values for those cleared dibits to + * calculate mean and std deviation of our Gaussians. With this we are ready to calculate the PDF of a new + * unknown analog value when required. + * Values that don't clear the Reed-Solomon check are discarded. + * This implementation uses a circular buffer to keep track of the N latest cleared analog dibits so we can + * adapt to changes in the signal. + * A modification was made to improve results for C4FM signals. See next block comment. + */ + +/** + * In the C4FM P25 recorded files from the "samples" repository, it can be observed that there is a + * correlation between the correct dibit associated for a given analog value and the value of the previous + * dibit. For instance, in one P25 recording, the dibits "0" come with an average analog signal of + * 3829 when the previous dibit was also "0," but if the previous dibit was a "3" then the average + * analog signal is 6875. These are the mean and std deviations for the full 4x4 combinations of previous and + * current dibits: + * + * 00: count: 200 mean: 3829.12 sd: 540.43 <- + * 01: count: 200 mean: 13352.45 sd: 659.74 + * 02: count: 200 mean: -5238.56 sd: 1254.70 + * 03: count: 200 mean: -13776.50 sd: 307.41 + * 10: count: 200 mean: 3077.74 sd: 1059.00 + * 11: count: 200 mean: 11935.11 sd: 776.20 + * 12: count: 200 mean: -6079.46 sd: 1003.94 + * 13: count: 200 mean: -13845.43 sd: 264.42 + * 20: count: 200 mean: 5574.33 sd: 1414.71 + * 21: count: 200 mean: 13687.75 sd: 727.68 + * 22: count: 200 mean: -4753.38 sd: 765.95 + * 23: count: 200 mean: -12342.17 sd: 1372.77 + * 30: count: 200 mean: 6875.23 sd: 1837.38 <- + * 31: count: 200 mean: 14527.99 sd: 406.85 + * 32: count: 200 mean: -3317.61 sd: 1089.02 + * 33: count: 200 mean: -12576.08 sd: 1161.77 + * || | | | + * || | | \_std deviation + * || | | + * || | \_mean of the current dibit + * || | + * || \_number of dibits used to calculate mean and std deviation + * || + * |\_current dibit + * | + * \_previous dibit + * + * This effect is not observed on QPSK or GFSK signals, there the mean values are quite consistent regardless + * of the previous dibit. + * + * The following define enables taking the previous dibit into account for C4FM signals. Comment out + * to disable. + */ +#define USE_PREVIOUS_DIBIT + +/** + * There is a minimum of cleared analog values we need to produce a meaningful mean and std deviation. + */ +#define MIN_ELEMENTS_FOR_HEURISTICS 10 + + +//Uncomment to disable the behaviour of this module. +//#define DISABLE_HEURISTICS + +/** + * The value of the previous dibit is only taken into account on the C4FM modulation. QPSK and GFSK are + * not improved by this technique. + */ +static int use_previous_dibit(int rf_mod) +{ + // 0: C4FM modulation + // 1: QPSK modulation + // 2: GFSK modulation + + // Use previoud dibit information when on C4FM + return (rf_mod == 0)? 1 : 0; +} + +/** + * Update the model of a symbol's Gaussian with new information. + * \param heuristics Pointer to the P25Heuristics module with all the needed state information. + * \param previous_dibit The cleared previous dibit value. + * \param original_dibit The current dibit as it was interpreted initially. + * \param dibit The current dibit. Will be different from original_dibit if the FEC fixed it. + * \param analog_value The actual analog signal value from which the original_dibit was derived. + */ +static void update_p25_heuristics(P25Heuristics* heuristics, int previous_dibit, int original_dibit, int dibit, int analog_value) +{ + float mean; + int old_value; + float old_mean; + + SymbolHeuristics* sh; + int number_errors; + +#ifndef USE_PREVIOUS_DIBIT + previous_dibit = 0; +#endif + + // Locate the Gaussian (SymbolHeuristics structure) we are going to update + sh = &(heuristics->symbols[previous_dibit][dibit]); + + // Update the circular buffers of values + old_value = sh->values[sh->index]; + old_mean = sh->means[sh->index]; + + // Update the BER statistics + number_errors = 0; + if (original_dibit != dibit) { + if ((original_dibit == 0 && dibit == 3) || (original_dibit == 3 && dibit == 0) || + (original_dibit == 1 && dibit == 2) || (original_dibit == 2 && dibit == 1)) { + // Interpreting a "00" as "11", "11" as "00", "01" as "10" or "10" as "01" counts as 2 errors + number_errors = 2; + } else { + // The other 8 combinations count (where original_dibit != dibit) as 1 error. + number_errors = 1; + } + } + update_error_stats(heuristics, 2, number_errors); + + // Update the running mean and variance. This is to calculate the PDF faster when required + if (sh->count >= HEURISTICS_SIZE) { + sh->sum -= old_value; + sh->var_sum -= (((float)old_value) - old_mean) * (((float)old_value) - old_mean); + } + sh->sum += analog_value; + + sh->values[sh->index] = analog_value; + if (sh->count < HEURISTICS_SIZE) { + sh->count++; + } + mean = sh->sum / ((float)sh->count); + sh->means[sh->index] = mean; + if (sh->index >= (HEURISTICS_SIZE-1)) { + sh->index = 0; + } else { + sh->index++; + } + + sh->var_sum += (((float)analog_value) - mean) * (((float)analog_value) - mean); +} + +void contribute_to_heuristics(int rf_mod, P25Heuristics* heuristics, AnalogSignal* analog_signal_array, int count) +{ + int i; + int use_prev_dibit; + +#ifdef USE_PREVIOUS_DIBIT + use_prev_dibit = use_previous_dibit(rf_mod); +#else + use_prev_dibit = 0; +#endif + + for (i=0; icount = 0; + sh->index = 0; + sh->sum = 0; + sh->var_sum = 0; +} + +void initialize_p25_heuristics(P25Heuristics* heuristics) +{ + int i, j; + for (i=0; i<4; i++) { + for (j=0; j<4; j++) { + initialize_symbol_heuristics(&(heuristics->symbols[i][j])); + } + } + heuristics->bit_count = 0; + heuristics->bit_error_count = 0; +} + +/** + * Important method to calculate the PDF (probability density function) of the Gaussian. + * TODO: improve performance. Since we are calculating this value to compare it with other PDF we can + * simplify very much. We don't really need to know the actual PDF value, just which Gaussian's got the + * highest PDF, which is a simpler problem. + */ +static float evaluate_pdf(SymbolHeuristics* se, int value) +{ + float x = (se->count*((float)value) - se->sum); + float y = -0.5F*x*x/(se->count*se->var_sum); + float pdf = sqrtf(se->count / se->var_sum) * expf(y) / sqrtf(2.0F * ((float) M_PI)); + + return pdf; +} + + +/** + * Logging of the internal PDF values for a given analog value and previous dibit. + */ +static void debug_log_pdf(P25Heuristics* heuristics, int previous_dibit, int analog_value) +{ + int i; + float pdfs[4]; + + for (i=0; i<4; i++) { + pdfs[i] = evaluate_pdf(&(heuristics->symbols[previous_dibit][i]), analog_value); + } + + fprintf(stderr, "v: %i, (%e, %e, %e, %e)\n", analog_value, pdfs[0], pdfs[1], pdfs[2], pdfs[3]); +} + +int estimate_symbol(int rf_mod, P25Heuristics* heuristics, int previous_dibit, int analog_value, int* dibit) +{ + int valid; + int i; + float pdfs[4]; + + +#ifdef USE_PREVIOUS_DIBIT + int use_prev_dibit = use_previous_dibit(rf_mod); + + if (use_prev_dibit == 0) + { + // Ignore + previous_dibit = 0; + } +#else + // Use previous_dibit as it comes. +#endif + + valid = 1; + + // Check if we have enough values to model the Gaussians for each symbol involved. + for (i=0; i<4; i++) { + if (heuristics->symbols[previous_dibit][i].count >= MIN_ELEMENTS_FOR_HEURISTICS) { + pdfs[i] = evaluate_pdf(&(heuristics->symbols[previous_dibit][i]), analog_value); + } else { + // Not enough data, we don't trust this result + valid = 0; + break; + } + } + + if (valid) { + // Find the highest pdf + int max_index; + float max; + + max_index = 0; + max = pdfs[0]; + for (i=1; i<4; i++) { + if (pdfs[i] > max) { + max_index = i; + max = pdfs[i]; + } + } + + // The symbol is the one with the highest pdf + *dibit = max_index; + } + +#ifdef DISABLE_HEURISTICS + valid = 0; +#endif + + return valid; +} + +/** + * Logs the internal state of the heuristic's state. Good for debugging. + */ +static void debug_print_symbol_heuristics(int previous_dibit, int dibit, SymbolHeuristics* sh) +{ + float mean, sd; + int k; + int n; + + n = sh->count; + if (n == 0) + { + mean = 0; + sd = 0; + } + else + { + mean = sh->sum/n; + sd = sqrtf(sh->var_sum / ((float) n)); + } + fprintf(stderr, "%i%i: count: %2i mean: % 10.2f sd: % 10.2f", previous_dibit, dibit, sh->count, mean, sd); + /* + fprintf(stderr, "("); + for (k=0; kvalues[k]); + } + fprintf(stderr, ")"); + */ + fprintf(stderr, "\n"); + +} + +void debug_print_heuristics(P25Heuristics* heuristics) +{ + int i,j; + + fprintf(stderr, "\n"); + + for(i=0; i<4; i++) + { + for(j=0; j<4; j++) + { + debug_print_symbol_heuristics(i, j, &(heuristics->symbols[i][j])); + } + } +} + +void update_error_stats(P25Heuristics* heuristics, int bits, int errors) +{ + heuristics->bit_count += bits; + heuristics->bit_error_count += errors; + + // Normalize to avoid overflow in the counters + if ((heuristics->bit_count & 1) == 0 && (heuristics->bit_error_count & 1) == 0) { + // We can divide both values by 2 safely. We just care about their ratio, not the actual value + heuristics->bit_count >>= 1; + heuristics->bit_error_count >>= 1; + } +} + +float get_P25_BER_estimate(P25Heuristics* heuristics) +{ + float ber; + if (heuristics->bit_count == 0) { + ber = 0.0F; + } else { + ber = ((float)heuristics->bit_error_count) * 100.0F / ((float)heuristics->bit_count); + } + return ber; +} diff --git a/plugins/channel/demoddsd/dsddecoder.cpp b/plugins/channel/demoddsd/dsddecoder.cpp index 45386f757..4184d9961 100644 --- a/plugins/channel/demoddsd/dsddecoder.cpp +++ b/plugins/channel/demoddsd/dsddecoder.cpp @@ -49,25 +49,9 @@ DSDDecoder::DSDDecoder() m_dsdParams.opts.mod_gfsk = 1; m_dsdParams.state.rf_mod = 0; - // Initialize the conditions - if(pthread_cond_init(&m_dsdParams.state.input_ready, NULL)) - { - qCritical("DSDDecoder::DSDDecoder: Unable to initialize input condition"); - } - if(pthread_cond_init(&m_dsdParams.state.output_ready, NULL)) - { - qCritical("DSDDecoder::DSDDecoder: Unable to initialize output condition"); - } - m_dsdParams.state.input_length = 0; m_dsdParams.state.input_offset = 0; - // Lock output mutex - if (pthread_mutex_lock(&m_dsdParams.state.output_mutex)) - { - qCritical("DSDDecoder::DSDDecoder: Unable to lock output mutex"); - } - m_dsdParams.state.output_buffer = (short *) malloc(1<<18); // Raw output buffer with single S16LE samples @ 8k (max: 128 kS) m_dsdParams.state.output_offset = 0; @@ -84,20 +68,10 @@ DSDDecoder::DSDDecoder() qCritical("DSDDecoder::DSDDecoder: Unable to allocate audio L+R buffer."); } - m_dsdThread = pthread_self(); // dummy initialization - m_dsdParams.state.dsd_running = 0; // wait for start() } DSDDecoder::~DSDDecoder() { - stop(); - - // Unlock output mutex - if (pthread_mutex_unlock(&m_dsdParams.state.output_mutex)) - { - qCritical("DSDDecoder::~DSDDecoder: Unable to unlock output mutex"); - } - free(m_dsdParams.state.output_samples); free(m_dsdParams.state.output_buffer); } @@ -109,61 +83,6 @@ void DSDDecoder::setInBuffer(const short *inBuffer) void DSDDecoder::pushSamples(int nbSamples) { - if (pthread_mutex_lock(&m_dsdParams.state.input_mutex)) - { - qCritical("DSDDecoder::pushSamples: Unable to lock input mutex"); - } - m_dsdParams.state.input_length = nbSamples; m_dsdParams.state.input_offset = 0; - - if (pthread_cond_signal(&m_dsdParams.state.input_ready)) - { - qCritical("DSDDecoder::pushSamples: Unable to signal input ready"); - } - - if (pthread_mutex_unlock(&m_dsdParams.state.input_mutex)) - { - qCritical("DSDDecoder::pushSamples: Unable to unlock input mutex"); - } } - -void DSDDecoder::start() -{ - if (m_dsdParams.state.dsd_running == 1) - { - m_dsdParams.state.dsd_running = 0; - - if (pthread_join(m_dsdThread, NULL)) { - qCritical("DSDDecoder::start: error joining DSD thread. Not starting"); - return; - } - } - - m_dsdParams.state.dsd_running = 1; - - if (pthread_create(&m_dsdThread, NULL, &run_dsd, &m_dsdParams)) - { - qCritical("DSDDecoder::start: Unable to spawn DSD thread"); - } -} - -void DSDDecoder::stop() -{ - if (m_dsdParams.state.dsd_running == 1) - { - m_dsdParams.state.dsd_running = 0; - - if (pthread_join(m_dsdThread, NULL)) { - qCritical("DSDDecoder::stop: error joining DSD thread. Not starting"); - } - } -} - -void* DSDDecoder::run_dsd (void *arg) -{ - dsd_params *params = (dsd_params *) arg; - liveScanner(¶ms->opts, ¶ms->state); - return NULL; -} - diff --git a/plugins/channel/demoddsd/dsddecoder.h b/plugins/channel/demoddsd/dsddecoder.h index ea44ca82f..1c4fd931d 100644 --- a/plugins/channel/demoddsd/dsddecoder.h +++ b/plugins/channel/demoddsd/dsddecoder.h @@ -29,11 +29,6 @@ public: void setInBuffer(const short *inBuffer); void pushSamples(int nbSamples); // Push this amount of samples to the DSD decoder thread - void start(); - void stop(); - - static void* run_dsd (void *arg); - private: typedef struct { @@ -42,12 +37,6 @@ private: } dsd_params; dsd_params m_dsdParams; - -// dsd_opts m_dsdOpts; -// dsd_state m_dsdState; - pthread_t m_dsdThread; }; - - #endif /* PLUGINS_CHANNEL_DEMODDSD_DSDDECODER_H_ */ diff --git a/plugins/channel/demoddsd/dsddemod.cpp b/plugins/channel/demoddsd/dsddemod.cpp index fa9de27fd..9d30478aa 100644 --- a/plugins/channel/demoddsd/dsddemod.cpp +++ b/plugins/channel/demoddsd/dsddemod.cpp @@ -201,12 +201,10 @@ void DSDDemod::start() { m_audioFifo.clear(); m_phaseDiscri.reset(); - m_dsdDecoder.start(); } void DSDDemod::stop() { - m_dsdDecoder.stop(); } bool DSDDemod::handleMessage(const Message& cmd)