rtl_433/src/pulse_detect.c

/**
 * Pulse detection functions
 *
 * Copyright (C) 2015 Tommy Vestermark
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 */

#include "pulse_detect.h"
#include "pulse_demod.h"
#include "util.h"
#include "rtl_433.h"
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>

void pulse_data_clear(pulse_data_t *data) {
	*data = (const pulse_data_t) {
		 0,
	};
}


void pulse_data_print(const pulse_data_t *data) {
    fprintf(stderr, "Pulse data: %u pulses\n", data->num_pulses);
	for(unsigned n = 0; n < data->num_pulses; ++n) {
		fprintf(stderr, "[%3u] Pulse: %4u, Gap: %4u, Period: %4u\n", n, data->pulse[n], data->gap[n], data->pulse[n] + data->gap[n]);
	}
}


// OOK adaptive level estimator constants
#define OOK_HIGH_LOW_RATIO	8			// Default ratio between high and low (noise) level
#define OOK_MIN_HIGH_LEVEL	1000		// Minimum estimate of high level
#define OOK_MAX_HIGH_LEVEL	(128*128)	// Maximum estimate for high level (A unit phasor is 128, anything above is overdrive)
#define OOK_MAX_LOW_LEVEL	(OOK_MAX_HIGH_LEVEL/2)	// Maximum estimate for low level
#define OOK_EST_HIGH_RATIO	64			// Constant for slowness of OOK high level estimator
#define OOK_EST_LOW_RATIO	1024		// Constant for slowness of OOK low level (noise) estimator (very slow)

// FSK adaptive frequency estimator constants
#define FSK_DEFAULT_FM_DELTA	6000	// Default estimate for frequency delta
#define FSK_EST_RATIO		32			// Constant for slowness of FSK estimators


/// Internal state data for pulse_FSK_detect()
typedef struct {
	unsigned int fsk_pulse_length;		// Counter for internal FSK pulse detection
	enum {
		PD_FSK_STATE_INIT	= 0,	// Initial frequency estimation
		PD_FSK_STATE_F1		= 1,	// High frequency (pulse)
		PD_FSK_STATE_F2		= 2,	// Low frequency (gap)
		PD_FSK_STATE_ERROR	= 3		// Error - stay here until cleared
	} fsk_state;

	int fm_f1_est;			// Estimate for the F1 frequency for FSK
	int fm_f2_est;			// Estimate for the F2 frequency for FSK

} pulse_FSK_state_t;


/// Demodulate Frequency Shift Keying (FSK) sample by sample
///
/// Function is stateful between calls
/// Builds estimate for initial frequency. When frequency deviates more than a
/// threshold value it will determine whether the deviation is positive or negative
/// to classify it as a pulse or gap. It will then transition to other state (F1 or F2)
/// and build an estimate of the other frequency. It will then transition back and forth when current
/// frequency is closer to other frequency estimate.
/// Includes spurious suppression by coalescing pulses when pulse/gap widths are too short.
/// Pulses equal higher frequency (F1) and Gaps equal lower frequency (F2)
/// @param fm_n: One single sample of FM data
/// @param *fsk_pulses: Will return a pulse_data_t structure for FSK demodulated data
/// @param *s: Internal state
void pulse_FSK_detect(int16_t fm_n, pulse_data_t *fsk_pulses, pulse_FSK_state_t *s) {
	const int fm_f1_delta = abs(fm_n - s->fm_f1_est);	// Get delta from F1 frequency estimate
	const int fm_f2_delta = abs(fm_n - s->fm_f2_est);	// Get delta from F2 frequency estimate
	s->fsk_pulse_length++;

	switch(s->fsk_state) {
		case PD_FSK_STATE_INIT:		// Initial frequency - High or low?
			// Initial samples?
			if (s->fsk_pulse_length < PD_MIN_PULSE_SAMPLES) {
				s->fm_f1_est = s->fm_f1_est/2 + fm_n/2;		// Quick initial estimator
			// Above default frequency delta?
			} else if (fm_f1_delta > (FSK_DEFAULT_FM_DELTA/2)) {
				// Positive frequency delta - Initial frequency was low (gap)
				if (fm_n > s->fm_f1_est) {
					s->fsk_state = PD_FSK_STATE_F1;
					s->fm_f2_est = s->fm_f1_est;	// Switch estimates
					s->fm_f1_est = fm_n;			// Prime F1 estimate
					fsk_pulses->pulse[0] = 0;		// Initial frequency was a gap...
					fsk_pulses->gap[0] = s->fsk_pulse_length;		// Store gap width
					fsk_pulses->num_pulses++;
					s->fsk_pulse_length = 0;
				// Negative Frequency delta - Initial frequency was high (pulse)
				} else {
					s->fsk_state = PD_FSK_STATE_F2;
					s->fm_f2_est = fm_n;	// Prime F2 estimate
					fsk_pulses->pulse[0] = s->fsk_pulse_length;	// Store pulse width
					s->fsk_pulse_length = 0;
				}
			// Still below threshold
			} else {
				s->fm_f1_est += fm_n/FSK_EST_RATIO - s->fm_f1_est/FSK_EST_RATIO;	// Slow estimator
			}
			break;
		case PD_FSK_STATE_F1:		// Pulse high at F1 frequency
			// Closer to F2 than F1?
			if (fm_f1_delta > fm_f2_delta) {
				s->fsk_state = PD_FSK_STATE_F2;
				// Store if pulse is not too short (suppress spurious)
				if (s->fsk_pulse_length >= PD_MIN_PULSE_SAMPLES) {
					fsk_pulses->pulse[fsk_pulses->num_pulses] = s->fsk_pulse_length;	// Store pulse width
					s->fsk_pulse_length = 0;
				// Else rewind to last gap
				} else {
					s->fsk_pulse_length += fsk_pulses->gap[fsk_pulses->num_pulses-1];	// Restore counter
					fsk_pulses->num_pulses--;		// Rewind one pulse
					// Are we back to initial frequency? (Was initial frequency a gap?)
					if ((fsk_pulses->num_pulses == 0) && (fsk_pulses->pulse[0] == 0)) {
						s->fm_f1_est = s->fm_f2_est;	// Switch back estimates
						s->fsk_state = PD_FSK_STATE_INIT;
					}
				}
			// Still below threshold
			} else {
				s->fm_f1_est += fm_n/FSK_EST_RATIO - s->fm_f1_est/FSK_EST_RATIO;	// Slow estimator
			}
			break;
		case PD_FSK_STATE_F2:		// Pulse gap at F2 frequency
			// Freq closer to F1 than F2 ?
			if (fm_f2_delta > fm_f1_delta) {
				s->fsk_state = PD_FSK_STATE_F1;
				// Store if pulse is not too short (suppress spurious)
				if (s->fsk_pulse_length >= PD_MIN_PULSE_SAMPLES) {
					fsk_pulses->gap[fsk_pulses->num_pulses] = s->fsk_pulse_length;	// Store gap width
					fsk_pulses->num_pulses++;	// Go to next pulse
					s->fsk_pulse_length = 0;
					// When pulse buffer is full go to error state
					if (fsk_pulses->num_pulses >= PD_MAX_PULSES) {
						fprintf(stderr, "pulse_FSK_detect(): Maximum number of pulses reached!\n");
						s->fsk_state = PD_FSK_STATE_ERROR;
					}
				// Else rewind to last pulse
				} else {
					s->fsk_pulse_length += fsk_pulses->pulse[fsk_pulses->num_pulses];	// Restore counter
					// Are we back to initial frequency?
					if (fsk_pulses->num_pulses == 0) {
						s->fsk_state = PD_FSK_STATE_INIT;
					}
				}
			// Still below threshold
			} else {
				s->fm_f2_est += fm_n/FSK_EST_RATIO - s->fm_f2_est/FSK_EST_RATIO;	// Slow estimator
			}
			break;
		case PD_FSK_STATE_ERROR:		// Stay here until cleared
			break;
		default:
			fprintf(stderr, "pulse_FSK_detect(): Unknown FSK state!!\n");
			s->fsk_state = PD_FSK_STATE_ERROR;
	} // switch(s->fsk_state)
}


/// Wrap up FSK modulation and store last data at End Of Package
///
/// @param fm_n: One single sample of FM data
/// @param *fsk_pulses: Pulse_data_t structure for FSK demodulated data
/// @param *s: Internal state
void pulse_FSK_wrap_up(pulse_data_t *fsk_pulses, pulse_FSK_state_t *s) {
	if (fsk_pulses->num_pulses < PD_MAX_PULSES) {	// Avoid overflow
		s->fsk_pulse_length++;
		if(s->fsk_state == PD_FSK_STATE_F1) {
			fsk_pulses->pulse[fsk_pulses->num_pulses] = s->fsk_pulse_length;	// Store last pulse
			fsk_pulses->gap[fsk_pulses->num_pulses] = 0;	// Zero gap at end
		} else {
			fsk_pulses->gap[fsk_pulses->num_pulses] = s->fsk_pulse_length;	// Store last gap
		}
		fsk_pulses->num_pulses++;
	}
}


/// Internal state data for pulse_pulse_package()
typedef struct {
	enum {
		PD_OOK_STATE_IDLE		= 0,
		PD_OOK_STATE_PULSE		= 1,
		PD_OOK_STATE_GAP_START	= 2,
		PD_OOK_STATE_GAP		= 3
	} ook_state;
	int pulse_length;		// Counter for internal pulse detection
	int max_pulse;			// Size of biggest pulse detected

	int data_counter;		// Counter for how much of data chunck is processed
	int lead_in_counter;	// Counter for allowing initial noise estimate to settle

	int ook_low_estimate;		// Estimate for the OOK low level (base noise level) in the envelope data
	int ook_high_estimate;		// Estimate for the OOK high level

	pulse_FSK_state_t	FSK_state;

} pulse_state_t;
static pulse_state_t pulse_state;


/// Demodulate On/Off Keying (OOK) and Frequency Shift Keying (FSK) from an envelope signal
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) {
	const int samples_per_ms = samp_rate / 1000;
	pulse_state_t *s = &pulse_state;
	s->ook_high_estimate = max(s->ook_high_estimate, OOK_MIN_HIGH_LEVEL);	// Be sure to set initial minimum level

	// Process all new samples
	while(s->data_counter < len) {
		// Calculate OOK detection threshold and hysteresis
		const int16_t am_n = envelope_data[s->data_counter];
		int16_t ook_threshold = s->ook_low_estimate + (s->ook_high_estimate - s->ook_low_estimate) / 2;
		if (level_limit != 0) ook_threshold = level_limit;	// Manual override
		const int16_t ook_hysteresis = ook_threshold / 8;	// ±12%

		// OOK State machine
		switch (s->ook_state) {
			case PD_OOK_STATE_IDLE:
				if (am_n > (ook_threshold + ook_hysteresis)	// Above threshold?
				 && s->lead_in_counter > OOK_EST_LOW_RATIO	// Lead in counter to stabilize noise estimate
				 ) {
					// Initialize all data
					pulse_data_clear(pulses);
					pulse_data_clear(fsk_pulses);
					s->pulse_length = 0;
					s->max_pulse = 0;
					s->FSK_state = (pulse_FSK_state_t){0};
					s->ook_state = PD_OOK_STATE_PULSE;
				} else {	// We are still idle..
					// Estimate low (noise) level
					const int ook_low_delta = am_n - s->ook_low_estimate;
					s->ook_low_estimate += ook_low_delta / OOK_EST_LOW_RATIO;
					s->ook_low_estimate += ((ook_low_delta > 0) ? 1 : -1);	// Hack to compensate for lack of fixed-point scaling
					// Calculate default OOK high level estimate
					s->ook_high_estimate = OOK_HIGH_LOW_RATIO * s->ook_low_estimate;	// Default is a ratio of low level
					s->ook_high_estimate = max(s->ook_high_estimate, OOK_MIN_HIGH_LEVEL);
					s->ook_high_estimate = min(s->ook_high_estimate, OOK_MAX_HIGH_LEVEL);
					if (s->lead_in_counter <= OOK_EST_LOW_RATIO) s->lead_in_counter++;		// Allow inital estimate to settle
				}
				break;
			case PD_OOK_STATE_PULSE:
				s->pulse_length++;
				// End of pulse detected?
				if (am_n  < (ook_threshold - ook_hysteresis)) {	// Gap?
					// Check for spurious short pulses
					if (s->pulse_length < PD_MIN_PULSE_SAMPLES) {
						s->ook_state = PD_OOK_STATE_IDLE;
					} else {
						// Continue with OOK decoding
						pulses->pulse[pulses->num_pulses] = s->pulse_length;	// Store pulse width
						s->max_pulse = max(s->pulse_length, s->max_pulse);	// Find largest pulse
						s->pulse_length = 0;
						s->ook_state = PD_OOK_STATE_GAP_START;
					}
				// Still pulse
				} else {
					// Calculate OOK high level estimate
					s->ook_high_estimate += am_n / OOK_EST_HIGH_RATIO - s->ook_high_estimate / OOK_EST_HIGH_RATIO;
					s->ook_high_estimate = max(s->ook_high_estimate, OOK_MIN_HIGH_LEVEL);
					s->ook_high_estimate = min(s->ook_high_estimate, OOK_MAX_HIGH_LEVEL);
					// Estimate pulse carrier frequency
					pulses->fsk_f1_est += fm_data[s->data_counter] / OOK_EST_HIGH_RATIO - pulses->fsk_f1_est / OOK_EST_HIGH_RATIO;
				}
				// FSK Demodulation
				if(pulses->num_pulses == 0) {	// Only during first pulse
					pulse_FSK_detect(fm_data[s->data_counter], fsk_pulses, &s->FSK_state);
				}
				break;
			case PD_OOK_STATE_GAP_START:	// Beginning of gap - it might be a spurious gap
				s->pulse_length++;
				// Pulse detected again already? (This is a spurious short gap)
				if (am_n  > (ook_threshold + ook_hysteresis)) {	// New pulse?
					s->pulse_length += pulses->pulse[pulses->num_pulses];	// Restore counter
					s->ook_state = PD_OOK_STATE_PULSE;
				// Or this gap is for real?
				} else if (s->pulse_length >= PD_MIN_PULSE_SAMPLES) {
					s->ook_state = PD_OOK_STATE_GAP;
					// Determine if FSK modulation is detected
					if(fsk_pulses->num_pulses > PD_MIN_PULSES) {
						// Store last pulse/gap
						pulse_FSK_wrap_up(fsk_pulses, &s->FSK_state);
						// Store estimates
						fsk_pulses->fsk_f1_est = s->FSK_state.fm_f1_est;
						fsk_pulses->fsk_f2_est = s->FSK_state.fm_f2_est;
						fsk_pulses->ook_low_estimate = s->ook_low_estimate;
						fsk_pulses->ook_high_estimate = s->ook_high_estimate;
						s->ook_state = PD_OOK_STATE_IDLE;	// Ensure everything is reset
						return 2;	// FSK package detected!!!
					}
				} // if
				// FSK Demodulation (continue during short gap - we might return...)
				if(pulses->num_pulses == 0) {	// Only during first pulse
					pulse_FSK_detect(fm_data[s->data_counter], fsk_pulses, &s->FSK_state);
				}
				break;
			case PD_OOK_STATE_GAP:
				s->pulse_length++;
				// New pulse detected?
				if (am_n  > (ook_threshold + ook_hysteresis)) {	// New pulse?
					pulses->gap[pulses->num_pulses] = s->pulse_length;	// Store gap width
					pulses->num_pulses++;	// Next pulse

					// EOP if too many pulses
					if (pulses->num_pulses >= PD_MAX_PULSES) {
						s->ook_state = PD_OOK_STATE_IDLE;
						// Store estimates
						pulses->ook_low_estimate = s->ook_low_estimate;
						pulses->ook_high_estimate = s->ook_high_estimate;
						return 1;	// End Of Package!!
					}

					s->pulse_length = 0;
					s->ook_state = PD_OOK_STATE_PULSE;
				}

				// EOP if gap is too long
				if (((s->pulse_length > (PD_MAX_GAP_RATIO * s->max_pulse))	// gap/pulse ratio exceeded
				 && (s->pulse_length > (PD_MIN_GAP_MS * samples_per_ms)))	// Minimum gap exceeded
				 || (s->pulse_length > (PD_MAX_GAP_MS * samples_per_ms))	// maximum gap exceeded
				 ) {
					pulses->gap[pulses->num_pulses] = s->pulse_length;	// Store gap width
					pulses->num_pulses++;	// Store last pulse
					s->ook_state = PD_OOK_STATE_IDLE;
					// Store estimates
					pulses->ook_low_estimate = s->ook_low_estimate;
					pulses->ook_high_estimate = s->ook_high_estimate;
					return 1;	// End Of Package!!
				}
				break;
			default:
				fprintf(stderr, "demod_OOK(): Unknown state!!\n");
				s->ook_state = PD_OOK_STATE_IDLE;
		} // switch
		s->data_counter++;
	} // while

	s->data_counter = 0;
	return 0;	// Out of data
}


#define MAX_HIST_BINS 16

/// Histogram data for single bin
typedef struct {
	unsigned count;
	int sum;
	int mean;
	int min;
	int max;
} hist_bin_t;

/// Histogram data for all bins
typedef struct {
	unsigned bins_count;
	hist_bin_t bins[MAX_HIST_BINS];
} histogram_t;


/// Generate a histogram (unsorted)
void histogram_sum(histogram_t *hist, const int *data, unsigned len, float tolerance) {
	unsigned bin;	// Iterator will be used outside for!

	for(unsigned n = 0; n < len; ++n) {
		// Search for match in existing bins
		for(bin = 0; bin < hist->bins_count; ++bin) {
			int bn = data[n];
			int bm = hist->bins[bin].mean;
			if (abs(bn - bm) < (tolerance * max(bn, bm))) {
				hist->bins[bin].count++;
				hist->bins[bin].sum += data[n];
				hist->bins[bin].mean = hist->bins[bin].sum / hist->bins[bin].count;
				hist->bins[bin].min	= min(data[n], hist->bins[bin].min);
				hist->bins[bin].max	= max(data[n], hist->bins[bin].max);
				break;	// Match found! Data added to existing bin
			}
		}
		// No match found? Add new bin
		if(bin == hist->bins_count && bin < MAX_HIST_BINS) {
			hist->bins[bin].count	= 1;
			hist->bins[bin].sum		= data[n];
			hist->bins[bin].mean	= data[n];
			hist->bins[bin].min		= data[n];
			hist->bins[bin].max		= data[n];
			hist->bins_count++;
		} // for bin
	} // for data
}


/// Delete bin from histogram
void histogram_delete_bin(histogram_t *hist, unsigned index) {
	const hist_bin_t	zerobin = {0};
	if (hist->bins_count < 1) return;	// Avoid out of bounds
	// Move all bins afterwards one forward
	for(unsigned n = index; n < hist->bins_count-1; ++n) {
		hist->bins[n] = hist->bins[n+1];
	}
	hist->bins_count--;
	hist->bins[hist->bins_count] = zerobin;	// Clear previously last bin
}


/// Swap two bins in histogram
void histogram_swap_bins(histogram_t *hist, unsigned index1, unsigned index2) {
	hist_bin_t	tempbin;
	if ((index1 < hist->bins_count) && (index2 < hist->bins_count)) {		// Avoid out of bounds
		tempbin = hist->bins[index1];
		hist->bins[index1] = hist->bins[index2];
		hist->bins[index2] = tempbin;
	}
}


/// Sort histogram with mean value (order lowest to highest)
void histogram_sort_mean(histogram_t *hist) {
	if (hist->bins_count < 2) return;		// Avoid underflow
	// Compare all bins (bubble sort)
	for(unsigned n = 0; n < hist->bins_count-1; ++n) {
		for(unsigned m = n+1; m < hist->bins_count; ++m) {
			if (hist->bins[m].mean < hist->bins[n].mean) {
				histogram_swap_bins(hist, m, n);
			} // if
		} // for m
	} // for n
}


/// Sort histogram with count value (order lowest to highest)
void histogram_sort_count(histogram_t *hist) {
	if (hist->bins_count < 2) return;		// Avoid underflow
	// Compare all bins (bubble sort)
	for(unsigned n = 0; n < hist->bins_count-1; ++n) {
		for(unsigned m = n+1; m < hist->bins_count; ++m) {
			if (hist->bins[m].count < hist->bins[n].count) {
				histogram_swap_bins(hist, m, n);
			} // if
		} // for m
	} // for n
}


/// Fuse histogram bins with means within tolerance
void histogram_fuse_bins(histogram_t *hist, float tolerance) {
	if (hist->bins_count < 2) return;		// Avoid underflow
	// Compare all bins
	for(unsigned n = 0; n < hist->bins_count-1; ++n) {
		for(unsigned m = n+1; m < hist->bins_count; ++m) {
			int bn = hist->bins[n].mean;
			int bm = hist->bins[m].mean;
			if (abs(bn - bm) < (tolerance * max(bn, bm))) {
				// Fuse data for bin[n] and bin[m]
				hist->bins[n].count += hist->bins[m].count;
				hist->bins[n].sum	+= hist->bins[m].sum;
				hist->bins[n].mean 	= hist->bins[n].sum / hist->bins[n].count;
				hist->bins[n].min 	= min(hist->bins[n].min, hist->bins[m].min);
				hist->bins[n].max 	= max(hist->bins[n].max, hist->bins[m].max);
				// Delete bin[m]
				histogram_delete_bin(hist, m);
				m--;	// Compare new bin in same place!
			} // if within tolerance
		} // for m
	} // for n
}


/// Print a histogram
void histogram_print(const histogram_t *hist, uint32_t samp_rate) {
	for(unsigned n = 0; n < hist->bins_count; ++n) {
		fprintf(stderr, " [%2u] count: %4u,  width: %5u [%2u;%2u]\t(%4.0f us)\n", n,
			hist->bins[n].count,
			hist->bins[n].mean,
			hist->bins[n].min,
			hist->bins[n].max,
			1E6f * hist->bins[n].mean / samp_rate);
	}
}


#define TOLERANCE (0.2)		// 20% tolerance should still discern between the pulse widths: 0.33, 0.66, 1.0

/// Analyze the statistics of a pulse data structure and print result
void pulse_analyzer(pulse_data_t *data, uint32_t samp_rate)
{
	// Generate pulse period data
	int pulse_total_period = 0;
	pulse_data_t pulse_periods = {0};
	pulse_periods.num_pulses = data->num_pulses;
	for(unsigned n = 0; n < pulse_periods.num_pulses; ++n) {
		pulse_periods.pulse[n] = data->pulse[n] + data->gap[n];
		pulse_total_period += data->pulse[n] + data->gap[n];
	}
	pulse_total_period -= data->gap[pulse_periods.num_pulses-1];

	histogram_t hist_pulses = {0};
	histogram_t hist_gaps = {0};
	histogram_t hist_periods = {0};

	// Generate statistics
	histogram_sum(&hist_pulses, data->pulse, data->num_pulses, TOLERANCE);
	histogram_sum(&hist_gaps, data->gap, data->num_pulses-1, TOLERANCE);						// Leave out last gap (end)
	histogram_sum(&hist_periods, pulse_periods.pulse, pulse_periods.num_pulses-1, TOLERANCE);	// Leave out last gap (end)

	// Fuse overlapping bins
	histogram_fuse_bins(&hist_pulses, TOLERANCE);
	histogram_fuse_bins(&hist_gaps, TOLERANCE);
	histogram_fuse_bins(&hist_periods, TOLERANCE);

	fprintf(stderr, "Analyzing pulses...\n");
	fprintf(stderr, "Total count: %4u,  width: %5i\t\t(%4.1f ms)\n",
		data->num_pulses, pulse_total_period, 1000.0f*pulse_total_period/samp_rate);
	fprintf(stderr, "Pulse width distribution:\n");
	histogram_print(&hist_pulses, samp_rate);
	fprintf(stderr, "Gap width distribution:\n");
	histogram_print(&hist_gaps, samp_rate);
	fprintf(stderr, "Pulse period distribution:\n");
	histogram_print(&hist_periods, samp_rate);
	fprintf(stderr, "Level estimates [high, low]: %6i, %6i\n",
		data->ook_high_estimate, data->ook_low_estimate);
	fprintf(stderr, "Frequency offsets [F1, F2]:  %6i, %6i\t(%+.1f kHz, %+.1f kHz)\n",
		data->fsk_f1_est, data->fsk_f2_est,
		(float)data->fsk_f1_est/INT16_MAX*samp_rate/2.0/1000.0,
		(float)data->fsk_f2_est/INT16_MAX*samp_rate/2.0/1000.0);

	fprintf(stderr, "Guessing modulation: ");
	struct protocol_state device = { .name = "Analyzer Device", 0};
	histogram_sort_mean(&hist_pulses);	// Easier to work with sorted data
	histogram_sort_mean(&hist_gaps);
	if(hist_pulses.bins[0].mean == 0) { histogram_delete_bin(&hist_pulses, 0); }	// Remove FSK initial zero-bin

	// Attempt to find a matching modulation
	if(data->num_pulses == 1) {
		fprintf(stderr, "Single pulse detected. Probably Frequency Shift Keying or just noise...\n");
	} else if(hist_pulses.bins_count == 1 && hist_gaps.bins_count == 1) {
		fprintf(stderr, "Un-modulated signal. Maybe a preamble...\n");
	} else if(hist_pulses.bins_count == 1 && hist_gaps.bins_count > 1) {
		fprintf(stderr, "Pulse Position Modulation with fixed pulse width\n");
		device.modulation	= OOK_PULSE_PPM_RAW;
		device.short_limit	= (hist_gaps.bins[0].mean + hist_gaps.bins[1].mean) / 2;	// Set limit between two lowest gaps
		device.long_limit	= hist_gaps.bins[1].max + 1;								// Set limit above next lower gap
		device.reset_limit	= hist_gaps.bins[hist_gaps.bins_count-1].max + 1;			// Set limit above biggest gap
	} else if(hist_pulses.bins_count == 2 && hist_gaps.bins_count == 1) {
		fprintf(stderr, "Pulse Width Modulation with fixed gap\n");
		device.modulation	= OOK_PULSE_PWM_RAW;
		device.short_limit	= (hist_pulses.bins[0].mean + hist_pulses.bins[1].mean) / 2;	// Set limit between two pulse widths
		device.long_limit	= hist_gaps.bins[hist_gaps.bins_count-1].max + 1;				// Set limit above biggest gap
		device.reset_limit	= device.long_limit;
	} else if(hist_pulses.bins_count == 2 && hist_gaps.bins_count == 2 && hist_periods.bins_count == 1) {
		fprintf(stderr, "Pulse Width Modulation with fixed period\n");
		device.modulation	= OOK_PULSE_PWM_RAW;
		device.short_limit	= (hist_pulses.bins[0].mean + hist_pulses.bins[1].mean) / 2;	// Set limit between two pulse widths
		device.long_limit	= hist_gaps.bins[hist_gaps.bins_count-1].max + 1;				// Set limit above biggest gap
		device.reset_limit	= device.long_limit;
	} else if(hist_pulses.bins_count == 2 && hist_gaps.bins_count == 2 && hist_periods.bins_count == 3) {
		fprintf(stderr, "Manchester coding\n");
		device.modulation	= OOK_PULSE_MANCHESTER_ZEROBIT;
		device.short_limit	= min(hist_pulses.bins[0].mean, hist_pulses.bins[1].mean);		// Assume shortest pulse is half period
		device.long_limit	= 0; // Not used
		device.reset_limit	= hist_gaps.bins[hist_gaps.bins_count-1].max + 1;				// Set limit above biggest gap
	} else if((hist_pulses.bins_count >= 3 && hist_gaps.bins_count >= 3)
		&& (abs(hist_pulses.bins[1].mean - 2*hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)	// Pulses are multiples of shortest pulse
		&& (abs(hist_pulses.bins[2].mean - 3*hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)
		&& (abs(hist_gaps.bins[0].mean   -   hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)	// Gaps are multiples of shortest pulse
		&& (abs(hist_gaps.bins[1].mean   - 2*hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)
		&& (abs(hist_gaps.bins[2].mean   - 3*hist_pulses.bins[0].mean) <= hist_pulses.bins[0].mean/8)
	) {
		fprintf(stderr, "Pulse Code Modulation (Not Return to Zero)\n");
		device.modulation	= FSK_PULSE_PCM;
		device.short_limit	= hist_pulses.bins[0].mean;			// Shortest pulse is bit width
		device.long_limit	= hist_pulses.bins[0].mean;			// Bit period equal to pulse length (NRZ)
		device.reset_limit	= hist_pulses.bins[0].mean*1024;	// No limit to run of zeros...
	} else if(hist_pulses.bins_count == 3) {
		fprintf(stderr, "Pulse Width Modulation with startbit/delimiter\n");
		device.modulation	= OOK_PULSE_PWM_TERNARY;
		device.short_limit	= (hist_pulses.bins[0].mean + hist_pulses.bins[1].mean) / 2;	// Set limit between two lowest pulse widths
		device.long_limit	= (hist_pulses.bins[1].mean + hist_pulses.bins[2].mean) / 2;	// Set limit between two next lowest pulse widths
		device.reset_limit	= hist_gaps.bins[hist_gaps.bins_count-1].max  +1;				// Set limit above biggest gap
		// Re-sort to find lowest pulse count index (is probably delimiter)
		histogram_sort_count(&hist_pulses);
		if		(hist_pulses.bins[0].mean < device.short_limit)	{	device.demod_arg = 0; }	// Shortest pulse is delimiter
		else if	(hist_pulses.bins[0].mean < device.long_limit)	{	device.demod_arg = 1; }	// Middle pulse is delimiter
		else													{	device.demod_arg = 2; }	// Longest pulse is delimiter
	} else {
		fprintf(stderr, "No clue...\n");
	}

	// Demodulate (if detected)
	if(device.modulation) {
		fprintf(stderr, "Attempting demodulation... short_limit: %.0f, long_limit: %.0f, reset_limit: %.0f, demod_arg: %lu\n",
			device.short_limit, device.long_limit, device.reset_limit, device.demod_arg);
		switch(device.modulation) {
			case FSK_PULSE_PCM:
				pulse_demod_pcm(data, &device);
				break;
			case OOK_PULSE_PPM_RAW:
				data->gap[data->num_pulses-1] = device.reset_limit + 1;	// Be sure to terminate package
				pulse_demod_ppm(data, &device);
				break;
			case OOK_PULSE_PWM_RAW:
				data->gap[data->num_pulses-1] = device.reset_limit + 1;	// Be sure to terminate package
				pulse_demod_pwm(data, &device);
				break;
			case OOK_PULSE_PWM_TERNARY:
				data->gap[data->num_pulses-1] = device.reset_limit + 1;	// Be sure to terminate package
				pulse_demod_pwm_ternary(data, &device);
				break;
			case OOK_PULSE_MANCHESTER_ZEROBIT:
				data->gap[data->num_pulses-1] = device.reset_limit + 1;	// Be sure to terminate package
				pulse_demod_manchester_zerobit(data, &device);
				break;
			default:
				fprintf(stderr, "Unsupported\n");
		}
	}

	fprintf(stderr, "\n");
}