it's all broken but it's almost not broken

This commit is contained in:
Fluora Eigenwire 2020-03-12 09:25:07 -05:00
parent 5fffcbce43
commit 777d1ada81
6 changed files with 429 additions and 144 deletions

View File

@ -35,7 +35,7 @@ pub fn generate
new_chirp.push(
(
32000.0 *
32767.0 *
env *
(
PI * t *

View File

@ -2,28 +2,34 @@
use crate::posi;
use crate::chirp;
use crate::gauss;
pub struct Demodulator
{
symbol_rise: Vec<i16>, // chirp from lower_freq to upper_freq, encoding a 1
symbol_fall: Vec<i16>, // chirp from upper_freq to lower_freq, encoding a 0
signal_bufsize: usize,
symbol_rise_norm: i64,
symbol_fall_norm: i64,
signal_buffer: posi::Buffer<i16>, // buffer holding the incoming signal to be convolved with the rising and falling chirps
signal_abs_sum: i64,
signal_buffer: posi::Buffer<i16>,
correl_bufsize: usize,
correl_rise_lopass: gauss::Filter,
correl_fall_lopass: gauss::Filter,
correl_rise_buffer: posi::Buffer<i64>, // buffer holding the values of the convolution with the rising chirp
correl_rise_abs_sum: i64,
correl_rise_abs_max: Vec<i64>, // structure for tracking the absolute maximum value currently contained in this buffer
correl_rise_hipass: gauss::Filter,
correl_fall_hipass: gauss::Filter,
correl_fall_buffer: posi::Buffer<i64>, // buffer holding the values of the convolution with the falling chirp
correl_fall_abs_sum: i64,
correl_fall_abs_max: Vec<i64>,
correl_rise_buffer: posi::Buffer<i64>,
correl_fall_buffer: posi::Buffer<i64>,
correl_rise_maxima: Vec<i64>,
correl_fall_maxima: Vec<i64>,
holdoff: usize,
accumulator_byte: u8,
accumulator_bit: u8,
samples_since_bit: usize,
abandon_samples: usize,
}
impl Demodulator
@ -31,7 +37,7 @@ impl Demodulator
pub fn new
(
center_freq: f32, // center frequency of the signal, referenced to the Nyquist limit
deviation: f32, // half-width of the signal chirps, referenced to the center frequency. (note: actual bandwidth will be higher due to sidebands.)
deviation: f32, // half-width of the signal chirps, referenced to the Nyquist limit. (note: actual bandwidth will be higher due to sidebands.)
symbol_len: f32, // length of a symbol (a chirp) in samples
) ->
Demodulator
@ -39,58 +45,72 @@ impl Demodulator
let lower_freq = (center_freq - deviation).max(0.0).min(1.0);
let upper_freq = (center_freq + deviation).max(0.0).min(1.0);
// the signals are, in fact, just mirror images of each other, but we're coding this
// as though they are not in case this changes later.
let mut symbol_rise:Vec<i16> = chirp::generate(lower_freq, upper_freq, symbol_len);
let mut symbol_fall:Vec<i16> = chirp::generate(upper_freq, lower_freq, symbol_len);
// important!!
// since the posi buffers shift from low to high index, their signal contents is reversed.
// we need to reverse our templates, too, or our convolutions will be backwards and produce inverted data.
symbol_rise.reverse();
symbol_fall.reverse();
// the signals are, in fact, just mirror images of each other, but we're coding this
// as though they are not in case this changes later.
let symbol_rise:Vec<i16> = chirp::generate(upper_freq, lower_freq, symbol_len);
let symbol_fall:Vec<i16> = chirp::generate(lower_freq, upper_freq, symbol_len);
let symbol_rise_norm:i64 =
(
symbol_rise
.iter()
.map(
|&x| (x as i64).abs()
)
.sum()
);
let symbol_fall_norm:i64 =
(
symbol_fall
.iter()
.map(
|&x| (x as i64).abs()
)
.sum()
);
let signal_bufsize = symbol_len.round() as usize;
let correl_bufsize = (symbol_len * 0.9).round() as usize;
let correl_bufsize = (symbol_len * 1.9).round() as usize;
Demodulator
{
symbol_rise: symbol_rise,
symbol_fall: symbol_fall,
signal_bufsize: signal_bufsize,
symbol_rise_norm: symbol_rise_norm,
symbol_fall_norm: symbol_fall_norm,
signal_buffer: posi::Buffer::new(signal_bufsize),
signal_abs_sum: 0,
correl_bufsize: correl_bufsize,
correl_rise_lopass: gauss::Filter::new(deviation / 2.0),
correl_fall_lopass: gauss::Filter::new(deviation / 2.0),
correl_rise_hipass: gauss::Filter::new(deviation / 16.0),
correl_fall_hipass: gauss::Filter::new(deviation / 16.0),
correl_rise_buffer: posi::Buffer::new(correl_bufsize),
correl_rise_abs_sum: 0,
correl_rise_abs_max: Vec::new(),
correl_fall_buffer: posi::Buffer::new(correl_bufsize),
correl_fall_abs_sum: 0,
correl_fall_abs_max: Vec::new(),
holdoff: correl_bufsize,
correl_rise_maxima: Vec::new(),
correl_fall_maxima: Vec::new(),
accumulator_byte: 0x00,
accumulator_bit: 0,
samples_since_bit: 0,
abandon_samples: (symbol_len * 17.0 / 16.0).round() as usize,
}
}
fn ingest
fn ingest_signal
(
&mut self,
new_sample: i16,
) {
let old_sample = self.signal_buffer.last();
self.signal_abs_sum += new_sample.abs() as i64;
self.signal_abs_sum -= old_sample.abs() as i64;
self.signal_buffer.shift_in(new_sample);
if self.signal_abs_sum == 0
{
self.holdoff = self.signal_bufsize;
}
}
fn correlate
@ -130,83 +150,95 @@ impl Demodulator
fall_value as i64
);
}
let (new_value_rise, new_value_fall) = (correl_rise, correl_fall);
let (old_value_rise, old_value_fall) = (self.correl_rise_buffer.last(), self.correl_fall_buffer.last());
self.correl_rise_abs_sum += new_value_rise.abs();
self.correl_rise_abs_sum -= old_value_rise.abs();
self.correl_fall_abs_sum += new_value_fall.abs();
self.correl_fall_abs_sum -= old_value_fall.abs();
correl_rise *= 32767;
correl_rise /= self.symbol_rise_norm;
if let Some(index) = self.correl_rise_abs_max.iter().position(|x| x == &old_value_rise.abs())
correl_fall *= 32767;
correl_fall /= self.symbol_fall_norm;
correl_rise = correl_rise.abs();
correl_fall = correl_fall.abs();
correl_rise = self.correl_rise_lopass.lo_process(correl_rise);
correl_fall = self.correl_fall_lopass.lo_process(correl_fall);
correl_rise = self.correl_rise_hipass.hi_process(correl_rise);
correl_fall = self.correl_fall_hipass.hi_process(correl_fall);
let old_rise = self.correl_rise_buffer.last();
let old_fall = self.correl_fall_buffer.last();
if let Ok(index) = self.correl_rise_maxima.binary_search(&old_rise)
{
self.correl_rise_abs_max.remove(index);
self.correl_rise_maxima.remove(index);
}
if let Some(index) = self.correl_fall_abs_max.iter().position(|x| x == &old_value_fall.abs())
if let Ok(index) = self.correl_fall_maxima.binary_search(&old_fall)
{
self.correl_fall_abs_max.remove(index);
self.correl_fall_maxima.remove(index);
}
if (
new_value_rise >= *self.correl_rise_abs_max.last().unwrap_or(&0)
&&
new_value_rise != 0
self.correl_rise_maxima.len() == 0
||
Some(&correl_rise) >= self.correl_rise_maxima.last()
) {
self.correl_rise_abs_max.push(new_value_rise);
self.correl_rise_maxima.push(correl_rise);
}
if (
new_value_fall >= *self.correl_fall_abs_max.last().unwrap_or(&0)
&&
new_value_fall != 0
self.correl_fall_maxima.len() == 0
||
Some(&correl_fall) >= self.correl_fall_maxima.last()
) {
self.correl_fall_abs_max.push(new_value_fall);
self.correl_fall_maxima.push(correl_fall);
}
self.correl_rise_buffer.shift_in(new_value_rise);
self.correl_fall_buffer.shift_in(new_value_fall);
self.correl_rise_buffer.shift_in(correl_rise);
self.correl_fall_buffer.shift_in(correl_fall);
}
fn determine
(
&self,
&mut self,
) ->
Option<bool>
{
let rise_avg = (self.correl_rise_abs_sum / self.correl_bufsize as i64);
let fall_avg = (self.correl_fall_abs_sum / self.correl_bufsize as i64);
let rise_value:i64 = self.correl_rise_buffer.mid();
let fall_value:i64 = self.correl_fall_buffer.mid();
// these values will stay at zero if no interval-maximum is found for their corresponding symbol template.
// otherwise, they will be set to that symbol's convolution magnitude.
let mut rise_weight:i64 = 0;
let mut fall_weight:i64 = 0;
let rise_max = *self.correl_rise_maxima.last().unwrap_or(&0);
let fall_max = *self.correl_fall_maxima.last().unwrap_or(&0);
let current_rise = self.correl_rise_buffer.mid().abs();
if (
current_rise > rise_avg * 4
&&
self.correl_rise_abs_max.last() == Some(&current_rise)
) {
rise_weight = current_rise;
}
let current_fall = self.correl_fall_buffer.mid().abs();
if (
current_fall > fall_avg * 4
&&
self.correl_fall_abs_max.last() == Some(&current_fall)
) {
fall_weight = current_fall;
}
if rise_weight != 0 || fall_weight != 0
{
return Some(fall_weight < rise_weight);
} else {
if rise_max == 0 && fall_max == 0
{
return None;
}
let rise_apex:bool = rise_value * 9 / 8 >= rise_max;
let fall_apex:bool = fall_value * 9 / 8 >= fall_max;
let rise_more:bool = rise_value >= fall_value * 9 / 8;
let fall_more:bool = fall_value >= rise_value * 9 / 8;
return match (rise_apex, rise_more, fall_apex, fall_more)
{
(false, _ , false, _ ) => None,
( true, true , _, false) => Some(true ),
( _, false, true, true ) => Some(false),
( _, _ , _, _ ) => None,
};
}
pub fn dump_correl
(
&self,
) ->
(i64, i64)
{
(
self.correl_rise_buffer.mid(),
self.correl_fall_buffer.mid(),
)
}
pub fn process
@ -214,28 +246,40 @@ impl Demodulator
&mut self,
new_sample: i16,
) ->
Option<bool>
Option<u8>
{
// general principle of operation:
// 1. convolve input signal with 1-symbol and 0-symbol template signals.
// 2. divide the convolution values by the total magnitude of the input signal over the convolution interval.
// 3. continuously check the normalized convolution curve for peaks which are, for an interval of one
// symbol-length centered on the peak in question, both:
// a. the tallest peak in the interval.
// b. at least three times the average value in the interval.
// 4. upon finding such a peak for either symbol template, emit a bit of demodulated data.
// if peaks coincide for both symbols, emit whichever bit has a higher convolution value.
self.ingest_signal(new_sample);
self.ingest(new_sample);
self.correlate();
self.correlate();
if self.holdoff > 0
if let Some(bit) = self.determine()
{
self.holdoff -= 1;
return None;
} else {
return self.determine();
self.accumulator_byte |= (bit as u8).reverse_bits() >> self.accumulator_bit;
self.accumulator_bit += 1;
self.samples_since_bit = 0;
}
else
{
self.samples_since_bit += 1;
if self.samples_since_bit > self.abandon_samples
{
self.accumulator_byte = 0x00;
self.accumulator_bit = 0;
}
}
let mut result:Option<u8> = None;
if self.accumulator_bit >= 8
{
result = Some(self.accumulator_byte);
self.accumulator_byte = 0x00;
self.accumulator_bit = 0;
}
return result;
}
}

View File

@ -2,7 +2,7 @@
use crate::chirp;
pub struct Modulator
pub struct Enmodulator
{
symbol_len: usize,
guard_pad: usize,
@ -10,7 +10,7 @@ pub struct Modulator
symbol_fall: Vec<i16>,
}
impl Modulator
impl Enmodulator
{
pub fn new
(
@ -18,7 +18,7 @@ impl Modulator
deviation: f32,
symbol_len: f32,
) ->
Modulator
Enmodulator
{
let lower_freq = (center_freq - deviation).max(0.0).min(1.0);
let upper_freq = (center_freq + deviation).max(0.0).min(1.0);
@ -26,7 +26,7 @@ impl Modulator
let symbol_rise:Vec<i16> = chirp::generate(lower_freq, upper_freq, symbol_len);
let symbol_fall:Vec<i16> = chirp::generate(upper_freq, lower_freq, symbol_len);
Modulator
Enmodulator
{
symbol_len: symbol_len.round() as usize,
guard_pad: (symbol_len / 16.0).round() as usize,

137
src/gauss.rs Normal file
View File

@ -0,0 +1,137 @@
use std::f32::consts::PI;
use crate::posi;
fn gaussian
(
width: f32,
) ->
Vec<i64>
{
let radius = (width * 9.4).round().max(1.0) as isize;
let mut output:Vec<i64> = Vec::with_capacity(radius as usize * 2 + 1);
for t in
(
(
-radius
..=
radius
).map(
|x| x as f32
)
) {
output.push(
(
32768.0 *
(
-0.5 * (t / width).powi(2)
)
.exp()
/
(
width * (2.0 * PI).sqrt()
)
)
.round()
as i64
);
}
return output;
}
pub struct Filter
{
kernel: Vec<i64>,
buffer: posi::Buffer<i64>,
}
impl Filter
{
pub fn new
(
frequency: f32,
) ->
Filter
{
let width:f32 = 1.0 / (PI * frequency);
let kernel = gaussian(width);
let buffer = posi::Buffer::new(kernel.len());
Filter
{
kernel: kernel,
buffer: buffer,
}
}
pub fn hi_process
(
&mut self,
new_value: i64,
) ->
i64
{
self.buffer.shift_in(new_value);
let mut output_value:i64 = 0;
for (&buffer_value, &kernel_value) in
(
self.buffer.to_vec().iter()
).zip(
self.kernel.iter()
) {
output_value +=
(
buffer_value
*
kernel_value
);
}
output_value /= 32768;
output_value =
(
self.buffer.mid() as i64
-
output_value
);
return output_value;
}
pub fn lo_process
(
&mut self,
new_value: i64,
) ->
i64
{
self.buffer.shift_in(new_value);
let mut output_value:i64 = 0;
for (&buffer_value, &kernel_value) in
(
self.buffer.to_vec().iter()
).zip(
self.kernel.iter()
) {
output_value +=
(
buffer_value
*
kernel_value
);
}
output_value /= 32768;
return output_value;
}
}

View File

@ -1,14 +1,19 @@
#![allow(unused_parens)]
use std::io::prelude::*;
use std::io::ErrorKind;
use std::io::BufReader;
use std::fs::File;
use std::time::Duration;
use std::time::Instant;
use std::process::exit;
extern crate lapp;
mod chirp;
mod gauss;
mod posi;
mod modulation;
mod enmodulation;
mod demodulation;
fn main
@ -28,13 +33,15 @@ fn main
Data input/output is arbitrary binary data, which can be text, file contents, or whatever else.
Usage:
<mode> (string) - Either 'e' to encode/modulate, or 'd' to decode/demodulate.
<mode> (string) - Either 'e' to encode/enmodulate, or 'd' to decode/demodulate.
-i, --input (default stdin) - Filelike object from which to read input.
-o, --output (default stdout) - Filelike object into which to write output.
-s, --sample-rate (default 48000.0) - Sample rate in Hertz, used for calculating modulation frequencies
-b, --bit-rate (default 480.0) - Data rate, in chirps (bits) per second. Be aware that lower values will increase system load.
-f, --center-frequency (default 12000.0) - Center frequency of the signal in Hertz.
-d, --deviation (default 1000.0) - Maximum distance traversed above and below the center frequency by the signal.
-v, --debug (default 0) - if specified when demodulating, emits every nth correlation value in f32 form to threnodyne_rise.dat and threnodyne_fall.dat.
"
)
);
@ -48,11 +55,13 @@ fn main
let arg_bitrate= args.get_float("bit-rate");
let arg_centerfreq = args.get_float("center-frequency");
let arg_deviation = args.get_float("deviation");
let arg_debug = args.get_integer("debug");
if !(
&["e", "d"].contains(&mode.as_str())
) {
eprintln!("threnodyne: please specify either 'e' to encode/modulate or 'd' to decode/demodulate.");
eprintln!("threnodyne: please specify either 'e' to encode/enmodulate or 'd' to decode/demodulate.");
exit(1);
}
@ -99,9 +108,13 @@ fn main
if mode.as_str() == "e"
{
let modulator = modulation::Modulator::new(center_freq, deviation, symbol_len);
let enmodulator = enmodulation::Enmodulator::new(center_freq, deviation, symbol_len);
let mut output_buffer:Vec<u8> = Vec::with_capacity(symbol_size);
let wall_start = Instant::now();
let mut process_time = Duration::new(0,0);
let mut samples_processed:usize = 0;
loop
{
@ -116,8 +129,10 @@ fn main
},
Ok(_) => input_buffer[0],
};
let iter_time = Instant::now();
let signal = modulator.process(byte);
let signal = enmodulator.process(byte);
output_buffer.clear();
@ -126,6 +141,9 @@ fn main
output_buffer.extend_from_slice(&sample.to_le_bytes());
}
process_time += iter_time.elapsed();
samples_processed += signal.len();
match output_file.write_all(&output_buffer)
{
Err(why) =>
@ -136,50 +154,104 @@ fn main
Ok(_) => (),
};
}
eprintln!();
eprintln!("samples processed: {}", samples_processed);
eprintln!("wall time: {} ms", wall_start.elapsed().as_millis());
eprintln!("process time: {} ms", process_time.as_millis());
eprintln!("process rate: {} samples/s", samples_processed as f32 / process_time.as_secs_f32());
}
if mode.as_str() == "d"
{
let mut demodulator = demodulation::Demodulator::new(center_freq, deviation, symbol_len);
let mut accumulator_byte:u8 = 0x00;
let mut accumulator_bit:u8 = 0;
let mut samples_since_bit:usize = 0;
let wall_start = Instant::now();
let mut process_time = Duration::new(0,0);
let mut samples_processed:usize = 0;
let guard_pad = symbol_size / 16;
let mut debug_files:Option<(File,File)> = None;
if arg_debug > 0
{
debug_files =
Some(
(
match File::create("threnodyne_rise.dat")
{
Ok(file) => file,
Err(why) =>
{
eprintln!("threnodyne: could not create debug output file for rising correlation: {}", why);
exit(1);
},
},
match File::create("threnodyne_fall.dat")
{
Ok(file) => file,
Err(why) =>
{
eprintln!("threnodyne: could not create debug output file for falling correlation: {}", why);
exit(1);
},
},
)
);
}
let mut ending:bool = false;
let mut ending_counter:usize = symbol_len as usize;
loop
{
let mut sample_bytes:[u8;2] = [0x00, 0x00];
match input_reader.read_exact(&mut sample_bytes)
if !ending
{
Err(why) if why.kind() == ErrorKind::UnexpectedEof => break,
Err(why) =>
match input_reader.read_exact(&mut sample_bytes)
{
eprintln!("threnodyne: input read error: {}", why);
exit(1);
},
Ok(_) => (),
Err(why) if why.kind() == ErrorKind::UnexpectedEof =>
{
ending = true
},
Err(why) =>
{
eprintln!("threnodyne: input read error: {}", why);
exit(1);
},
Ok(_) => (),
}
}
else if ending_counter > 0
{
ending_counter -= 1;
}
else
{
break;
}
let iter_time = Instant::now();
let sample = i16::from_le_bytes(sample_bytes);
if let Some(bit) = demodulator.process(sample)
let demod_result = demodulator.process(sample);
if let Some((risefile, fallfile)) = &mut debug_files
{
accumulator_byte |= [0b0000_0000, 0b1000_0000] [bit as usize] >> accumulator_bit;
accumulator_bit += 1;
samples_since_bit = 0;
} else {
samples_since_bit += 1;
if samples_since_bit >= (symbol_size + guard_pad)
if samples_processed % (arg_debug as usize) == 0
{
accumulator_bit = 0;
accumulator_byte = 0x00;
let (riseval, fallval) = demodulator.dump_correl();
let riseval_f = (riseval as f32) / 32768.0;
let fallval_f = (fallval as f32) / 32768.0;
risefile.write_all(&riseval_f.to_be_bytes()).expect("write failed on rise debug file");
fallfile.write_all(&fallval_f.to_be_bytes()).expect("write failed on fall debug file");
}
}
if accumulator_bit >= 8
process_time += iter_time.elapsed();
samples_processed += 1;
if let Some(byte) = demod_result
{
match output_file.write_all(&[accumulator_byte])
match output_file.write_all(&[byte])
{
Err(why) =>
{
@ -188,12 +260,15 @@ fn main
},
Ok(_) => (),
};
accumulator_bit = 0;
accumulator_byte = 0x00;
let _ = output_file.flush();
}
}
eprintln!();
eprintln!("samples processed: {}", samples_processed);
eprintln!("wall time: {} ms", wall_start.elapsed().as_millis());
eprintln!("process time: {} ms", process_time.as_millis());
eprintln!("process rate: {} samples/s", samples_processed as f32 / process_time.as_secs_f32());
}
}

View File

@ -59,6 +59,15 @@ impl<T> Buffer<T>
}
}
pub fn first
(
&self,
) ->
T
{
self.get(0)
}
pub fn mid
(
&self,
@ -67,6 +76,26 @@ impl<T> Buffer<T>
{
self.get(self.size / 2)
}
pub fn above_mid
(
&self,
index: usize
) ->
T
{
self.get((self.size / 2).wrapping_add(index))
}
pub fn below_mid
(
&self,
index: usize
) ->
T
{
self.get((self.size / 2).wrapping_sub(index))
}
pub fn last
(