src/pulse_detect.c

   1 /**
   2  * Pulse detection functions
   3  *
   4  * Copyright (C) 2015 Tommy Vestermark
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 2 of the License, or
   8  * (at your option) any later version.
   9  */
  10
  11 #include "pulse_detect.h"
  12 #include "pulse_demod.h"
  13 #include "util.h"
  14 #include "rtl_433.h"
  15 #include <limits.h>
  16 #include <stdio.h>
  17 #include <stdlib.h>
  18
  19 void pulse_data_clear(pulse_data_t *data) {
  20         *data = (const pulse_data_t) {
  21                  0,
  22         };
  23 }
  24
  25
  26 void pulse_data_print(const pulse_data_t *data) {
  27     fprintf(stderr, "Pulse data: %u pulses\n", data->num_pulses);
  28         for(unsigned n = 0; n < data->num_pulses; ++n) {
  29                 fprintf(stderr, "[%3u] Pulse: %4u, Gap: %4u, Period: %4u\n", n, data->pulse[n], data->gap[n], data->pulse[n] + data->gap[n]);
  30         }
  31 }
  32
  33
  34 // OOK adaptive level estimator constants
  35 #define OOK_HIGH_LOW_RATIO      8                       // Default ratio between high and low (noise) level
  36 #define OOK_MIN_HIGH_LEVEL      1000            // Minimum estimate of high level
  37 #define OOK_MAX_HIGH_LEVEL      (128*128)       // Maximum estimate for high level (A unit phasor is 128, anything above is overdrive)
  38 #define OOK_MAX_LOW_LEVEL       (OOK_MAX_HIGH_LEVEL/2)  // Maximum estimate for low level
  39 #define OOK_EST_HIGH_RATIO      64                      // Constant for slowness of OOK high level estimator
  40 #define OOK_EST_LOW_RATIO       1024            // Constant for slowness of OOK low level (noise) estimator (very slow)
  41
  42 // FSK adaptive frequency estimator constants
  43 #define FSK_DEFAULT_FM_DELTA    6000    // Default estimate for frequency delta
  44 #define FSK_EST_RATIO           32                      // Constant for slowness of FSK estimators
  45
  46
  47 /// Internal state data for pulse_FSK_detect()
  48 typedef struct {
  49         unsigned int fsk_pulse_length;          // Counter for internal FSK pulse detection
  50         enum {
  51                 PD_FSK_STATE_INIT       = 0,    // Initial frequency estimation
  52                 PD_FSK_STATE_F1         = 1,    // High frequency (pulse)
  53                 PD_FSK_STATE_F2         = 2,    // Low frequency (gap)
  54                 PD_FSK_STATE_ERROR      = 3             // Error - stay here until cleared
  55         } fsk_state;
  56
  57         int fm_f1_est;                  // Estimate for the F1 frequency for FSK
  58         int fm_f2_est;                  // Estimate for the F2 frequency for FSK
  59
  60 } pulse_FSK_state_t;
  61
  62
  63 /// Demodulate Frequency Shift Keying (FSK) sample by sample
  64 ///
  65 /// Function is stateful between calls
  66 /// Builds estimate for initial frequency. When frequency deviates more than a
  67 /// threshold value it will determine whether the deviation is positive or negative
  68 /// to classify it as a pulse or gap. It will then transition to other state (F1 or F2)
  69 /// and build an estimate of the other frequency. It will then transition back and forth when current
  70 /// frequency is closer to other frequency estimate.
  71 /// Includes spurious suppression by coalescing pulses when pulse/gap widths are too short.
  72 /// Pulses equal higher frequency (F1) and Gaps equal lower frequency (F2)
  73 /// @param fm_n: One single sample of FM data
  74 /// @param *fsk_pulses: Will return a pulse_data_t structure for FSK demodulated data
  75 /// @param *s: Internal state
  76 void pulse_FSK_detect(int16_t fm_n, pulse_data_t *fsk_pulses, pulse_FSK_state_t *s) {
  77         const int fm_f1_delta = abs(fm_n - s->fm_f1_est);       // Get delta from F1 frequency estimate
  78         const int fm_f2_delta = abs(fm_n - s->fm_f2_est);       // Get delta from F2 frequency estimate
  79         s->fsk_pulse_length++;
  80
  81         switch(s->fsk_state) {
  82                 case PD_FSK_STATE_INIT:         // Initial frequency - High or low?
  83                         // Initial samples?
  84                         if (s->fsk_pulse_length < PD_MIN_PULSE_SAMPLES) {
  85                                 s->fm_f1_est = s->fm_f1_est/2 + fm_n/2;         // Quick initial estimator
  86                         // Above default frequency delta?
  87                         } else if (fm_f1_delta > (FSK_DEFAULT_FM_DELTA/2)) {
  88                                 // Positive frequency delta - Initial frequency was low (gap)
  89                                 if (fm_n > s->fm_f1_est) {
  90                                         s->fsk_state = PD_FSK_STATE_F1;
  91                                         s->fm_f2_est = s->fm_f1_est;    // Switch estimates
  92                                         s->fm_f1_est = fm_n;                    // Prime F1 estimate
  93                                         fsk_pulses->pulse[0] = 0;               // Initial frequency was a gap...
  94                                         fsk_pulses->gap[0] = s->fsk_pulse_length;               // Store gap width
  95                                         fsk_pulses->num_pulses++;
  96                                         s->fsk_pulse_length = 0;
  97                                 // Negative Frequency delta - Initial frequency was high (pulse)
  98                                 } else {
  99                                         s->fsk_state = PD_FSK_STATE_F2;
 100                                         s->fm_f2_est = fm_n;    // Prime F2 estimate
 101                                         fsk_pulses->pulse[0] = s->fsk_pulse_length;     // Store pulse width
 102                                         s->fsk_pulse_length = 0;
 103                                 }
 104                         // Still below threshold
 105                         } else {
 106                                 s->fm_f1_est += fm_n/FSK_EST_RATIO - s->fm_f1_est/FSK_EST_RATIO;        // Slow estimator
 107                         }
 108                         break;
 109                 case PD_FSK_STATE_F1:           // Pulse high at F1 frequency
 110                         // Closer to F2 than F1?
 111                         if (fm_f1_delta > fm_f2_delta) {
 112                                 s->fsk_state = PD_FSK_STATE_F2;
 113                                 // Store if pulse is not too short (suppress spurious)
 114                                 if (s->fsk_pulse_length >= PD_MIN_PULSE_SAMPLES) {
 115                                         fsk_pulses->pulse[fsk_pulses->num_pulses] = s->fsk_pulse_length;        // Store pulse width
 116                                         s->fsk_pulse_length = 0;
 117                                 // Else rewind to last gap
 118                                 } else {
 119                                         s->fsk_pulse_length += fsk_pulses->gap[fsk_pulses->num_pulses-1];       // Restore counter
 120                                         fsk_pulses->num_pulses--;               // Rewind one pulse
 121                                         // Are we back to initial frequency? (Was initial frequency a gap?)
 122                                         if ((fsk_pulses->num_pulses == 0) && (fsk_pulses->pulse[0] == 0)) {
 123                                                 s->fm_f1_est = s->fm_f2_est;    // Switch back estimates
 124                                                 s->fsk_state = PD_FSK_STATE_INIT;
 125                                         }
 126                                 }
 127                         // Still below threshold
 128                         } else {
 129                                 s->fm_f1_est += fm_n/FSK_EST_RATIO - s->fm_f1_est/FSK_EST_RATIO;        // Slow estimator
 130                         }
 131                         break;
 132                 case PD_FSK_STATE_F2:           // Pulse gap at F2 frequency
 133                         // Freq closer to F1 than F2 ?
 134                         if (fm_f2_delta > fm_f1_delta) {
 135                                 s->fsk_state = PD_FSK_STATE_F1;
 136                                 // Store if pulse is not too short (suppress spurious)
 137                                 if (s->fsk_pulse_length >= PD_MIN_PULSE_SAMPLES) {
 138                                         fsk_pulses->gap[fsk_pulses->num_pulses] = s->fsk_pulse_length;  // Store gap width
 139                                         fsk_pulses->num_pulses++;       // Go to next pulse
 140                                         s->fsk_pulse_length = 0;
 141                                         // When pulse buffer is full go to error state
 142                                         if (fsk_pulses->num_pulses >= PD_MAX_PULSES) {
 143                                                 fprintf(stderr, "pulse_FSK_detect(): Maximum number of pulses reached!\n");
 144                                                 s->fsk_state = PD_FSK_STATE_ERROR;
 145                                         }
 146                                 // Else rewind to last pulse
 147                                 } else {
 148                                         s->fsk_pulse_length += fsk_pulses->pulse[fsk_pulses->num_pulses];       // Restore counter
 149                                         // Are we back to initial frequency?
 150                                         if (fsk_pulses->num_pulses == 0) {
 151                                                 s->fsk_state = PD_FSK_STATE_INIT;
 152                                         }
 153                                 }
 154                         // Still below threshold
 155                         } else {
 156                                 s->fm_f2_est += fm_n/FSK_EST_RATIO - s->fm_f2_est/FSK_EST_RATIO;        // Slow estimator
 157                         }
 158                         break;
 159                 case PD_FSK_STATE_ERROR:                // Stay here until cleared
 160                         break;
 161                 default:
 162                         fprintf(stderr, "pulse_FSK_detect(): Unknown FSK state!!\n");
 163                         s->fsk_state = PD_FSK_STATE_ERROR;
 164         } // switch(s->fsk_state)
 165 }
 166
 167
 168 /// Wrap up FSK modulation and store last data at End Of Package
 169 ///
 170 /// @param fm_n: One single sample of FM data
 171 /// @param *fsk_pulses: Pulse_data_t structure for FSK demodulated data
 172 /// @param *s: Internal state
 173 void pulse_FSK_wrap_up(pulse_data_t *fsk_pulses, pulse_FSK_state_t *s) {
 174         if (fsk_pulses->num_pulses < PD_MAX_PULSES) {   // Avoid overflow
 175                 s->fsk_pulse_length++;
 176                 if(s->fsk_state == PD_FSK_STATE_F1) {
 177                         fsk_pulses->pulse[fsk_pulses->num_pulses] = s->fsk_pulse_length;        // Store last pulse
 178                         fsk_pulses->gap[fsk_pulses->num_pulses] = 0;    // Zero gap at end
 179                 } else {
 180                         fsk_pulses->gap[fsk_pulses->num_pulses] = s->fsk_pulse_length;  // Store last gap
 181                 }
 182                 fsk_pulses->num_pulses++;
 183         }
 184 }
 185
 186
 187 /// Internal state data for pulse_pulse_package()
 188 typedef struct {
 189         enum {
 190                 PD_OOK_STATE_IDLE               = 0,
 191                 PD_OOK_STATE_PULSE              = 1,
 192                 PD_OOK_STATE_GAP_START  = 2,
 193                 PD_OOK_STATE_GAP                = 3
 194         } ook_state;
 195         int pulse_length;               // Counter for internal pulse detection
 196         int max_pulse;                  // Size of biggest pulse detected
 197
 198         int data_counter;               // Counter for how much of data chunck is processed
 199         int lead_in_counter;    // Counter for allowing initial noise estimate to settle
 200
 201         int ook_low_estimate;           // Estimate for the OOK low level (base noise level) in the envelope data
 202         int ook_high_estimate;          // Estimate for the OOK high level
 203
 204         pulse_FSK_state_t       FSK_state;
 205
 206 } pulse_state_t;
 207 static pulse_state_t pulse_state;
 208
 209
 210 /// Demodulate On/Off Keying (OOK) and Frequency Shift Keying (FSK) from an envelope signal
 211 int pulse_detect_package(const int16_t *envelope_data, const int16_t *fm_data, int len, int16_t level_limit, uint32_t samp_rate, pulse_data_t *pulses, pulse_data_t *fsk_pulses) {
 212         const int samples_per_ms = samp_rate / 1000;
 213         pulse_state_t *s = &pulse_state;
 214         s->ook_high_estimate = max(s->ook_high_estimate, OOK_MIN_HIGH_LEVEL);   // Be sure to set initial minimum level
 215
 216         // Process all new samples
 217         while(s->data_counter < len) {
 218                 // Calculate OOK detection threshold and hysteresis
 219                 const int16_t am_n = envelope_data[s->data_counter];
 220                 int16_t ook_threshold = s->ook_low_estimate + (s->ook_high_estimate - s->ook_low_estimate) / 2;
 221                 if (level_limit != 0) ook_threshold = level_limit;      // Manual override
 222                 const int16_t ook_hysteresis = ook_threshold / 8;       // ±12%
 223
 224                 // OOK State machine
 225                 switch (s->ook_state) {
 226                         case PD_OOK_STATE_IDLE:
 227                                 if (am_n > (ook_threshold + ook_hysteresis)     // Above threshold?
 228                                  && s->lead_in_counter > OOK_EST_LOW_RATIO      // Lead in counter to stabilize noise estimate
 229                                  ) {
 230                                         // Initialize all data
 231                                         pulse_data_clear(pulses);
 232                                         pulse_data_clear(fsk_pulses);
 233                                         s->pulse_length = 0;
 234                                         s->max_pulse = 0;
 235                                         s->FSK_state = (pulse_FSK_state_t){0};
 236                                         s->ook_state = PD_OOK_STATE_PULSE;
 237                                 } else {        // We are still idle..
 238                                         // Estimate low (noise) level
 239                                         const int ook_low_delta = am_n - s->ook_low_estimate;
 240                                         s->ook_low_estimate += ook_low_delta / OOK_EST_LOW_RATIO;
 241                                         s->ook_low_estimate += ((ook_low_delta > 0) ? 1 : -1);  // Hack to compensate for lack of fixed-point scaling
 242                                         // Calculate default OOK high level estimate
 243                                         s->ook_high_estimate = OOK_HIGH_LOW_RATIO * s->ook_low_estimate;        // Default is a ratio of low level
 244                                         s->ook_high_estimate = max(s->ook_high_estimate, OOK_MIN_HIGH_LEVEL);
 245                                         s->ook_high_estimate = min(s->ook_high_estimate, OOK_MAX_HIGH_LEVEL);
 246                                         if (s->lead_in_counter <= OOK_EST_LOW_RATIO) s->lead_in_counter++;              // Allow inital estimate to settle
 247                                 }
 248                                 break;
 249                         case PD_OOK_STATE_PULSE:
 250                                 s->pulse_length++;
 251                                 // End of pulse detected?
 252                                 if (am_n  < (ook_threshold - ook_hysteresis)) { // Gap?
 253                                         // Check for spurious short pulses
 254                                         if (s->pulse_length < PD_MIN_PULSE_SAMPLES) {
 255                                                 s->ook_state = PD_OOK_STATE_IDLE;
 256                                         } else {
 257                                                 // Continue with OOK decoding
 258                                                 pulses->pulse[pulses->num_pulses] = s->pulse_length;    // Store pulse width
 259                                                 s->max_pulse = max(s->pulse_length, s->max_pulse);      // Find largest pulse
 260                                                 s->pulse_length = 0;
 261                                                 s->ook_state = PD_OOK_STATE_GAP_START;
 262                                         }
 263                                 // Still pulse
 264                                 } else {
 265                                         // Calculate OOK high level estimate
 266                                         s->ook_high_estimate += am_n / OOK_EST_HIGH_RATIO - s->ook_high_estimate / OOK_EST_HIGH_RATIO;
 267                                         s->ook_high_estimate = max(s->ook_high_estimate, OOK_MIN_HIGH_LEVEL);
 268                                         s->ook_high_estimate = min(s->ook_high_estimate, OOK_MAX_HIGH_LEVEL);
 269                                         // Estimate pulse carrier frequency
 270                                         pulses->fsk_f1_est += fm_data[s->data_counter] / OOK_EST_HIGH_RATIO - pulses->fsk_f1_est / OOK_EST_HIGH_RATIO;
 271                                 }
 272                                 // FSK Demodulation
 273                                 if(pulses->num_pulses == 0) {   // Only during first pulse
 274                                         pulse_FSK_detect(fm_data[s->data_counter], fsk_pulses, &s->FSK_state);
 275                                 }
 276                                 break;
 277                         case PD_OOK_STATE_GAP_START:    // Beginning of gap - it might be a spurious gap
 278                                 s->pulse_length++;
 279                                 // Pulse detected again already? (This is a spurious short gap)
 280                                 if (am_n  > (ook_threshold + ook_hysteresis)) { // New pulse?
 281                                         s->pulse_length += pulses->pulse[pulses->num_pulses];   // Restore counter
 282                                         s->ook_state = PD_OOK_STATE_PULSE;
 283                                 // Or this gap is for real?
 284                                 } else if (s->pulse_length >= PD_MIN_PULSE_SAMPLES) {
 285                                         s->ook_state = PD_OOK_STATE_GAP;
 286                                         // Determine if FSK modulation is detected
 287                                         if(fsk_pulses->num_pulses > PD_MIN_PULSES) {
 288                                                 // Store last pulse/gap
 289                                                 pulse_FSK_wrap_up(fsk_pulses, &s->FSK_state);
 290                                                 // Store estimates
 291                                                 fsk_pulses->fsk_f1_est = s->FSK_state.fm_f1_est;
 292                                                 fsk_pulses->fsk_f2_est = s->FSK_state.fm_f2_est;
 293                                                 fsk_pulses->ook_low_estimate = s->ook_low_estimate;
 294                                                 fsk_pulses->ook_high_estimate = s->ook_high_estimate;
 295                                                 s->ook_state = PD_OOK_STATE_IDLE;       // Ensure everything is reset
 296                                                 return 2;       // FSK package detected!!!
 297                                         }
 298                                 } // if
 299                                 // FSK Demodulation (continue during short gap - we might return...)
 300                                 if(pulses->num_pulses == 0) {   // Only during first pulse
 301                                         pulse_FSK_detect(fm_data[s->data_counter], fsk_pulses, &s->FSK_state);
 302                                 }
 303                                 break;
 304                         case PD_OOK_STATE_GAP:
 305                                 s->pulse_length++;
 306                                 // New pulse detected?
 307                                 if (am_n  > (ook_threshold + ook_hysteresis)) { // New pulse?
 308                                         pulses->gap[pulses->num_pulses] = s->pulse_length;      // Store gap width
 309                                         pulses->num_pulses++;   // Next pulse
 310
 311                                         // EOP if too many pulses
 312                                         if (pulses->num_pulses >= PD_MAX_PULSES) {
 313                                                 s->ook_state = PD_OOK_STATE_IDLE;
 314                                                 // Store estimates
 315                                                 pulses->ook_low_estimate = s->ook_low_estimate;
 316                                                 pulses->ook_high_estimate = s->ook_high_estimate;
 317                                                 return 1;       // End Of Package!!
 318                                         }
 319
 320                                         s->pulse_length = 0;
 321                                         s->ook_state = PD_OOK_STATE_PULSE;
 322                                 }
 323
 324                                 // EOP if gap is too long
 325                                 if (((s->pulse_length > (PD_MAX_GAP_RATIO * s->max_pulse))      // gap/pulse ratio exceeded
 326                                  && (s->pulse_length > (PD_MIN_GAP_MS * samples_per_ms)))       // Minimum gap exceeded
 327                                  || (s->pulse_length > (PD_MAX_GAP_MS * samples_per_ms))        // maximum gap exceeded
 328                                  ) {
 329                                         pulses->gap[pulses->num_pulses] = s->pulse_length;      // Store gap width
 330                                         pulses->num_pulses++;   // Store last pulse
 331                                         s->ook_state = PD_OOK_STATE_IDLE;
 332                                         // Store estimates
 333                                         pulses->ook_low_estimate = s->ook_low_estimate;
 334                                         pulses->ook_high_estimate = s->ook_high_estimate;
 335                                         return 1;       // End Of Package!!
 336                                 }
 337                                 break;
 338                         default:
 339                                 fprintf(stderr, "demod_OOK(): Unknown state!!\n");
 340                                 s->ook_state = PD_OOK_STATE_IDLE;
 341                 } // switch
 342                 s->data_counter++;
 343         } // while
 344
 345         s->data_counter = 0;
 346         return 0;       // Out of data
 347 }
 348
 349
 350 #define MAX_HIST_BINS 16
 351
 352 /// Histogram data for single bin
 353 typedef struct {
 354         unsigned count;
 355         int sum;
 356         int mean;
 357         int min;
 358         int max;
 359 } hist_bin_t;
 360
 361 /// Histogram data for all bins
 362 typedef struct {
 363         unsigned bins_count;
 364         hist_bin_t bins[MAX_HIST_BINS];
 365 } histogram_t;
 366
 367
 368 /// Generate a histogram (unsorted)
 369 void histogram_sum(histogram_t *hist, const int *data, unsigned len, float tolerance) {
 370         unsigned bin;   // Iterator will be used outside for!
 371
 372         for(unsigned n = 0; n < len; ++n) {
 373                 // Search for match in existing bins
 374                 for(bin = 0; bin < hist->bins_count; ++bin) {
 375                         int bn = data[n];
 376                         int bm = hist->bins[bin].mean;
 377                         if (abs(bn - bm) < (tolerance * max(bn, bm))) {
 378                                 hist->bins[bin].count++;
 379                                 hist->bins[bin].sum += data[n];
 380                                 hist->bins[bin].mean = hist->bins[bin].sum / hist->bins[bin].count;
 381                                 hist->bins[bin].min     = min(data[n], hist->bins[bin].min);
 382                                 hist->bins[bin].max     = max(data[n], hist->bins[bin].max);
 383                                 break;  // Match found! Data added to existing bin
 384                         }
 385                 }
 386                 // No match found? Add new bin
 387                 if(bin == hist->bins_count && bin < MAX_HIST_BINS) {
 388                         hist->bins[bin].count   = 1;
 389                         hist->bins[bin].sum             = data[n];
 390                         hist->bins[bin].mean    = data[n];
 391                         hist->bins[bin].min             = data[n];
 392                         hist->bins[bin].max             = data[n];
 393                         hist->bins_count++;
 394                 } // for bin
 395         } // for data
 396 }
 397
 398
 399 /// Delete bin from histogram
 400 void histogram_delete_bin(histogram_t *hist, unsigned index) {
 401         const hist_bin_t        zerobin = {0};
 402         if (hist->bins_count < 1) return;       // Avoid out of bounds
 403         // Move all bins afterwards one forward
 404         for(unsigned n = index; n < hist->bins_count-1; ++n) {
 405                 hist->bins[n] = hist->bins[n+1];
 406         }
 407         hist->bins_count--;
 408         hist->bins[hist->bins_count] = zerobin; // Clear previously last bin
 409 }
 410
 411
 412 /// Swap two bins in histogram
 413 void histogram_swap_bins(histogram_t *hist, unsigned index1, unsigned index2) {
 414         hist_bin_t      tempbin;
 415         if ((index1 < hist->bins_count) && (index2 < hist->bins_count)) {               // Avoid out of bounds
 416                 tempbin = hist->bins[index1];
 417                 hist->bins[index1] = hist->bins[index2];
 418                 hist->bins[index2] = tempbin;
 419         }
 420 }
 421
 422
 423 /// Sort histogram with mean value (order lowest to highest)
 424 void histogram_sort_mean(histogram_t *hist) {
 425         if (hist->bins_count < 2) return;               // Avoid underflow
 426         // Compare all bins (bubble sort)
 427         for(unsigned n = 0; n < hist->bins_count-1; ++n) {
 428                 for(unsigned m = n+1; m < hist->bins_count; ++m) {
 429                         if (hist->bins[m].mean < hist->bins[n].mean) {
 430                                 histogram_swap_bins(hist, m, n);
 431                         } // if
 432                 } // for m
 433         } // for n
 434 }
 435
 436
 437 /// Sort histogram with count value (order lowest to highest)
 438 void histogram_sort_count(histogram_t *hist) {
 439         if (hist->bins_count < 2) return;               // Avoid underflow
 440         // Compare all bins (bubble sort)
 441         for(unsigned n = 0; n < hist->bins_count-1; ++n) {
 442                 for(unsigned m = n+1; m < hist->bins_count; ++m) {
 443                         if (hist->bins[m].count < hist->bins[n].count) {
 444                                 histogram_swap_bins(hist, m, n);
 445                         } // if
 446                 } // for m
 447         } // for n
 448 }
 449
 450
 451 /// Fuse histogram bins with means within tolerance
 452 void histogram_fuse_bins(histogram_t *hist, float tolerance) {
 453         if (hist->bins_count < 2) return;               // Avoid underflow
 454         // Compare all bins
 455         for(unsigned n = 0; n < hist->bins_count-1; ++n) {
 456                 for(unsigned m = n+1; m < hist->bins_count; ++m) {
 457                         int bn = hist->bins[n].mean;
 458                         int bm = hist->bins[m].mean;
 459                         if (abs(bn - bm) < (tolerance * max(bn, bm))) {
 460                                 // Fuse data for bin[n] and bin[m]
 461                                 hist->bins[n].count += hist->bins[m].count;
 462                                 hist->bins[n].sum       += hist->bins[m].sum;
 463                                 hist->bins[n].mean      = hist->bins[n].sum / hist->bins[n].count;
 464                                 hist->bins[n].min       = min(hist->bins[n].min, hist->bins[m].min);
 465                                 hist->bins[n].max       = max(hist->bins[n].max, hist->bins[m].max);
 466                                 // Delete bin[m]
 467                                 histogram_delete_bin(hist, m);
 468                                 m--;    // Compare new bin in same place!
 469                         } // if within tolerance
 470                 } // for m
 471         } // for n
 472 }
 473
 474
 475 /// Print a histogram
 476 void histogram_print(const histogram_t *hist, uint32_t samp_rate) {
 477         for(unsigned n = 0; n < hist->bins_count; ++n) {
 478                 fprintf(stderr, " [%2u] count: %4u,  width: %5u [%2u;%2u]\t(%4.0f us)\n", n,
 479                         hist->bins[n].count,
 480                         hist->bins[n].mean,
 481                         hist->bins[n].min,
 482                         hist->bins[n].max,
 483                         1E6f * hist->bins[n].mean / samp_rate);
 484         }
 485 }
 486
 487
 488 #define TOLERANCE (0.2)         // 20% tolerance should still discern between the pulse widths: 0.33, 0.66, 1.0
 489
 490 /// Analyze the statistics of a pulse data structure and print result
 491 void pulse_analyzer(pulse_data_t *data, uint32_t samp_rate)
 492 {
 493         // Generate pulse period data
 494         int pulse_total_period = 0;
 495         pulse_data_t pulse_periods = {0};
 496         pulse_periods.num_pulses = data->num_pulses;
 497         for(unsigned n = 0; n < pulse_periods.num_pulses; ++n) {
 498                 pulse_periods.pulse[n] = data->pulse[n] + data->gap[n];
 499                 pulse_total_period += data->pulse[n] + data->gap[n];
 500         }
 501         pulse_total_period -= data->gap[pulse_periods.num_pulses-1];
 502
 503         histogram_t hist_pulses = {0};
 504         histogram_t hist_gaps = {0};
 505         histogram_t hist_periods = {0};
 506
 507         // Generate statistics
 508         histogram_sum(&hist_pulses, data->pulse, data->num_pulses, TOLERANCE);
 509         histogram_sum(&hist_gaps, data->gap, data->num_pulses-1, TOLERANCE);                                            // Leave out last gap (end)
 510         histogram_sum(&hist_periods, pulse_periods.pulse, pulse_periods.num_pulses-1, TOLERANCE);       // Leave out last gap (end)
 511
 512         // Fuse overlapping bins
 513         histogram_fuse_bins(&hist_pulses, TOLERANCE);
 514         histogram_fuse_bins(&hist_gaps, TOLERANCE);
 515         histogram_fuse_bins(&hist_periods, TOLERANCE);
 516
 517         fprintf(stderr, "Analyzing pulses...\n");
 518         fprintf(stderr, "Total count: %4u,  width: %5i\t\t(%4.1f ms)\n",
 519                 data->num_pulses, pulse_total_period, 1000.0f*pulse_total_period/samp_rate);
 520         fprintf(stderr, "Pulse width distribution:\n");
 521         histogram_print(&hist_pulses, samp_rate);
 522         fprintf(stderr, "Gap width distribution:\n");
 523         histogram_print(&hist_gaps, samp_rate);
 524         fprintf(stderr, "Pulse period distribution:\n");
 525         histogram_print(&hist_periods, samp_rate);
 526         fprintf(stderr, "Level estimates [high, low]: %6i, %6i\n",
 527                 data->ook_high_estimate, data->ook_low_estimate);
 528         fprintf(stderr, "Frequency offsets [F1, F2]:  %6i, %6i\t(%+.1f kHz, %+.1f kHz)\n",
 529                 data->fsk_f1_est, data->fsk_f2_est,
 530                 (float)data->fsk_f1_est/INT16_MAX*samp_rate/2.0/1000.0,
 531                 (float)data->fsk_f2_est/INT16_MAX*samp_rate/2.0/1000.0);
 532
 533         fprintf(stderr, "Guessing modulation: ");
 534         struct protocol_state device = { .name = "Analyzer Device", 0};
 535         histogram_sort_mean(&hist_pulses);      // Easier to work with sorted data
 536         histogram_sort_mean(&hist_gaps);
 537         if(hist_pulses.bins[0].mean == 0) { histogram_delete_bin(&hist_pulses, 0); }    // Remove FSK initial zero-bin
 538
 539         // Attempt to find a matching modulation
 540         if(data->num_pulses == 1) {
 541                 fprintf(stderr, "Single pulse detected. Probably Frequency Shift Keying or just noise...\n");
 542         } else if(hist_pulses.bins_count == 1 && hist_gaps.bins_count == 1) {
 543                 fprintf(stderr, "Un-modulated signal. Maybe a preamble...\n");
 544         } else if(hist_pulses.bins_count == 1 && hist_gaps.bins_count > 1) {
 545                 fprintf(stderr, "Pulse Position Modulation with fixed pulse width\n");
 546                 device.modulation       = OOK_PULSE_PPM_RAW;
 547                 device.short_limit      = (hist_gaps.bins[0].mean + hist_gaps.bins[1].mean) / 2;        // Set limit between two lowest gaps
 548                 device.long_limit       = hist_gaps.bins[1].max + 1;                                                            // Set limit above next lower gap
 549                 device.reset_limit      = hist_gaps.bins[hist_gaps.bins_count-1].max + 1;                       // Set limit above biggest gap
 550         } else if(hist_pulses.bins_count == 2 && hist_gaps.bins_count == 1) {
 551                 fprintf(stderr, "Pulse Width Modulation with fixed gap\n");
 552                 device.modulation       = OOK_PULSE_PWM_RAW;
 553                 device.short_limit      = (hist_pulses.bins[0].mean + hist_pulses.bins[1].mean) / 2;    // Set limit between two pulse widths
 554                 device.long_limit       = hist_gaps.bins[hist_gaps.bins_count-1].max + 1;                               // Set limit above biggest gap
 555                 device.reset_limit      = device.long_limit;
 556         } else if(hist_pulses.bins_count == 2 && hist_gaps.bins_count == 2 && hist_periods.bins_count == 1) {
 557                 fprintf(stderr, "Pulse Width Modulation with fixed period\n");
 558                 device.modulation       = OOK_PULSE_PWM_RAW;
 559                 device.short_limit      = (hist_pulses.bins[0].mean + hist_pulses.bins[1].mean) / 2;    // Set limit between two pulse widths
 560                 device.long_limit       = hist_gaps.bins[hist_gaps.bins_count-1].max + 1;                               // Set limit above biggest gap
 561                 device.reset_limit      = device.long_limit;
 562         } else if(hist_pulses.bins_count == 2 && hist_gaps.bins_count == 2 && hist_periods.bins_count == 3) {
 563                 fprintf(stderr, "Manchester coding\n");
 564                 device.modulation       = OOK_PULSE_MANCHESTER_ZEROBIT;
 565                 device.short_limit      = min(hist_pulses.bins[0].mean, hist_pulses.bins[1].mean);              // Assume shortest pulse is half period
 566                 device.long_limit       = 0; // Not used
 567                 device.reset_limit      = hist_gaps.bins[hist_gaps.bins_count-1].max + 1;                               // Set limit above biggest gap
 568         } else if((hist_pulses.bins_count >= 3 && hist_gaps.bins_count >= 3)
 569                 && (abs(hist_pulses.bins[1].mean - 2*hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)   // Pulses are multiples of shortest pulse
 570                 && (abs(hist_pulses.bins[2].mean - 3*hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)
 571                 && (abs(hist_gaps.bins[0].mean   -   hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)   // Gaps are multiples of shortest pulse
 572                 && (abs(hist_gaps.bins[1].mean   - 2*hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)
 573                 && (abs(hist_gaps.bins[2].mean   - 3*hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)
 574         ) {
 575                 fprintf(stderr, "Pulse Code Modulation (Not Return to Zero)\n");
 576                 device.modulation       = FSK_PULSE_PCM;
 577                 device.short_limit      = hist_pulses.bins[0].mean;                     // Shortest pulse is bit width
 578                 device.long_limit       = hist_pulses.bins[0].mean;                     // Bit period equal to pulse length (NRZ)
 579                 device.reset_limit      = hist_pulses.bins[0].mean*1024;        // No limit to run of zeros...
 580         } else if(hist_pulses.bins_count == 3) {
 581                 fprintf(stderr, "Pulse Width Modulation with startbit/delimiter\n");
 582                 device.modulation       = OOK_PULSE_PWM_TERNARY;
 583                 device.short_limit      = (hist_pulses.bins[0].mean + hist_pulses.bins[1].mean) / 2;    // Set limit between two lowest pulse widths
 584                 device.long_limit       = (hist_pulses.bins[1].mean + hist_pulses.bins[2].mean) / 2;    // Set limit between two next lowest pulse widths
 585                 device.reset_limit      = hist_gaps.bins[hist_gaps.bins_count-1].max  +1;                               // Set limit above biggest gap
 586                 // Re-sort to find lowest pulse count index (is probably delimiter)
 587                 histogram_sort_count(&hist_pulses);
 588                 if              (hist_pulses.bins[0].mean < device.short_limit) {       device.demod_arg = 0; } // Shortest pulse is delimiter
 589                 else if (hist_pulses.bins[0].mean < device.long_limit)  {       device.demod_arg = 1; } // Middle pulse is delimiter
 590                 else                                                                                                    {       device.demod_arg = 2; } // Longest pulse is delimiter
 591         } else {
 592                 fprintf(stderr, "No clue...\n");
 593         }
 594
 595         // Demodulate (if detected)
 596         if(device.modulation) {
 597                 fprintf(stderr, "Attempting demodulation... short_limit: %.0f, long_limit: %.0f, reset_limit: %.0f, demod_arg: %lu\n",
 598                         device.short_limit, device.long_limit, device.reset_limit, device.demod_arg);
 599                 switch(device.modulation) {
 600                         case FSK_PULSE_PCM:
 601                                 pulse_demod_pcm(data, &device);
 602                                 break;
 603                         case OOK_PULSE_PPM_RAW:
 604                                 data->gap[data->num_pulses-1] = device.reset_limit + 1; // Be sure to terminate package
 605                                 pulse_demod_ppm(data, &device);
 606                                 break;
 607                         case OOK_PULSE_PWM_RAW:
 608                                 data->gap[data->num_pulses-1] = device.reset_limit + 1; // Be sure to terminate package
 609                                 pulse_demod_pwm(data, &device);
 610                                 break;
 611                         case OOK_PULSE_PWM_TERNARY:
 612                                 data->gap[data->num_pulses-1] = device.reset_limit + 1; // Be sure to terminate package
 613                                 pulse_demod_pwm_ternary(data, &device);
 614                                 break;
 615                         case OOK_PULSE_MANCHESTER_ZEROBIT:
 616                                 data->gap[data->num_pulses-1] = device.reset_limit + 1; // Be sure to terminate package
 617                                 pulse_demod_manchester_zerobit(data, &device);
 618                                 break;
 619                         default:
 620                                 fprintf(stderr, "Unsupported\n");
 621                 }
 622         }
 623
 624         fprintf(stderr, "\n");
 625 }
 626