diff options
Diffstat (limited to 'octave/fsk_horus.m')
| -rw-r--r-- | octave/fsk_horus.m | 876 |
1 files changed, 876 insertions, 0 deletions
diff --git a/octave/fsk_horus.m b/octave/fsk_horus.m new file mode 100644 index 0000000..2b41cf1 --- /dev/null +++ b/octave/fsk_horus.m @@ -0,0 +1,876 @@ +% fsk_horus.m +% David Rowe 10 Oct 2015 +% +% Project Horus High Altitude Balloon (HAB) FSK demodulator +% See blog write up "All your modems are belong to us" +% http://www.rowetel.com/?p=4629 + + +fsk_lib; + + +% Basic modem set up for Horus + +function states = fsk_horus_init(Fs,Rs,M=2) + + states = fsk_init(Fs,Rs,M); + + % Freq. estimator limits - keep these narrow to stop errors with low SNR 4FSK + + states.fest_fmin = 300; + states.fest_fmax = 2800; +endfunction + + +% init rtty protocol specific states + +function rtty = fsk_horus_init_rtty + % Generate unique word that correlates against the ASCII "$$$$$" that + % is at the start of each frame. + % $ -> 36 decimal -> 0 1 0 0 1 0 0 binary + + dollar_bits = fliplr([0 1 0 0 1 0 0]); + mapped_db = 2*dollar_bits - 1; + sync_bits = [1 1 0]; + mapped_sb = 2*sync_bits - 1; + + mapped = [mapped_db mapped_sb]; + npad = rtty.npad = 3; % one start and two stop bits between 7 bit ascii chars + nfield = rtty.nfield = 7; % length of ascii character field + + rtty.uw = [mapped mapped mapped mapped mapped]; + rtty.uw_thresh = length(rtty.uw) - 2; % allow a few bit errors when looking for UW + rtty.max_packet_len = 1000; +endfunction + + +% I think this is the binary protocol work from Jan 2016 + +function binary = fsk_horus_init_binary + % Generate 16 bit "$$" unique word that is at the front of every horus binary + % packet + + dollar_bits = [0 0 1 0 0 1 0 0]; + mapped_db = 2*dollar_bits - 1; + + binary.uw = [mapped_db mapped_db]; + binary.uw_thresh = length(binary.uw)-2; % no bit errors when looking for UW + + binary.max_packet_len = 360; +endfunction + + +% Look for unique word and return index of first UW bit, or -1 if no +% UW found Sometimes there may be several matches, returns the +% position of the best match to UW. + +function [uw_start best_corr corr] = find_uw(states, start_bit, rx_bits) + uw = states.uw; + + mapped_rx_bits = 2*rx_bits - 1; + best_corr = 0; + uw_start = -1; + found_uw = 0; + + % first first UW in buffer that exceeds threshold + + for i=start_bit:length(rx_bits) - length(uw) + corr(i) = mapped_rx_bits(i:i+length(uw)-1) * uw'; + if (found_uw == 0) && (corr(i) >= states.uw_thresh) + uw_start = i; + best_corr = corr; + found_uw = 1; + end + end + +endfunction + + +% Extract ASCII string from a Horus frame of bits + +function [str crc_ok] = extract_ascii(states, rx_bits_buf, uw_loc) + nfield = states.nfield; + npad = states.npad; + + str = ""; str_dec = []; nstr = 0; ptx_crc = 1; rx_crc = ""; + endpacket = 0; + + st = uw_loc + length(states.uw); % first bit of first char + en = uw_loc + states.max_packet_len - nfield; + %printf("\nst: %d en: %d len: %d\n", st, en, length(rx_bits_buf)); + + for i=st:nfield+npad:en + field = rx_bits_buf(i:i+nfield-1); + ch_dec = field * (2.^(0:nfield-1))'; + + % filter out unlikely characters that bit errors may introduce, and ignore \n + + if (ch_dec > 31) && (ch_dec < 91) + str = [str char(ch_dec)]; + else + str = [str char(32)]; % space is "not sure" + end + nstr++; + + % build up array for CRC16 check + + if !endpacket && (ch_dec == 42) + endpacket = 1; + rx_crc = crc16(str_dec); % found a '*' so that's the end of the string for CRC calculations + ptx_crc = nstr+1; % this is where the transmit CRC starts + end + if !endpacket + str_dec = [str_dec ch_dec]; + end + end + + if (ptx_crc+3) <= length(str) + tx_crc = str(ptx_crc:ptx_crc+3); + crc_ok = strcmp(tx_crc, rx_crc); + else + crc_ok = 0; + end + + str = str(1:ptx_crc-2); + +endfunction + + +% Use soft decision information to find bits most likely in error. I think +% this is some form of maximum likelihood decoding. + +function [str crc_ok rx_bits_log_flipped] = sd_bit_flipping(states, rx_bits_log, rx_bits_sd_log, st, en); + + % force algorithm to ignore rs232 sync bits by marking them as "very likely", they have + % no input to crc algorithm + + nfield = states.nfield; + npad = states.npad; + for i=st:nfield+npad:en + rx_bits_sd_log(i+nfield:i+nfield+npad-1) = 1E6; + end + + % make a list of bits with smallest soft decn values + + [dodgy_bits_mag dodgy_bits_index] = sort(abs(rx_bits_sd_log(st+length(states.uw):en))); + dodgy_bits_index += length(states.uw) + st - 1; + nbits = 6; + ntries = 2^nbits; + str = ""; + crc_ok = 0; + + % try various combinations of these bits + + for i=1:ntries-1 + error_mask = zeros(1, length(rx_bits_log)); + for b=1:nbits + x = bitget(i,b); + bit_to_flip = dodgy_bits_index(b); + error_mask(bit_to_flip) = x; + %printf("st: %d i: %d b: %d x: %d index: %d\n", st, i,b,x,bit_to_flip); + end + rx_bits_log_flipped = xor(rx_bits_log, error_mask); + [str_flipped crc_ok_flipped] = extract_ascii(states, rx_bits_log_flipped, st); + if crc_ok_flipped + %printf("Yayy we fixed a packet by flipping with pattern %d\n", i); + str = str_flipped; + crc_ok = crc_ok_flipped; + end + end +endfunction + + +% Extract as many ASCII packets as we can from a great big buffer of bits + +function npackets = extract_and_print_rtty_packets(states, rtty, rx_bits_log, rx_bits_sd_log) + + % use UWs to delimit start and end of data packets + + bit = 1; + nbits = length(rx_bits_log); + nfield = rtty.nfield; + npad = rtty.npad; + npackets = 0; + + uw_loc = find_uw(rtty, bit, rx_bits_log, states.verbose); + + while (uw_loc != -1) + if bitand(states.verbose,0x8) + printf("nbits: %d max_packet_len: %d uw_loc: %d\n", nbits, rtty.max_packet_len, uw_loc); + end + + if (uw_loc + rtty.max_packet_len) < nbits + % Now start picking out 7 bit ascii chars from frame. It has some + % structure so we can guess where fields are. I hope we don't get + % RS232 idle bits stuck into it anywhere, ie "bit fields" don't + % change dynamically. + + % dump msg bits so we can use them as a test signal + %msg = rx_bits_log(st:uw_loc-1); + %save -ascii horus_msg.txt msg + + % simulate bit error for testing + %rx_bits_log(st+200) = xor(rx_bits_log(st+100),1); + %rx_bits_sd_log(st+100) = 0; + + [str crc_ok] = extract_ascii(rtty, rx_bits_log, uw_loc); + + if crc_ok == 0 + [str_flipped crc_flipped_ok rx_bits_log] = sd_bit_flipping(rtty, rx_bits_log, rx_bits_sd_log, uw_loc, uw_loc+rtty.max_packet_len); + end + + % update memory of previous packet, we use this to guess where errors may be + if crc_ok || crc_flipped_ok + states.prev_pkt = rx_bits_log(uw_loc+length(rtty.uw):uw_loc+rtty.max_packet_len); + end + + if crc_ok + str = sprintf("%s CRC OK", str); + npackets++; + else + if crc_flipped_ok + str = sprintf("%s fixed", str_flipped); + else + str = sprintf("%s CRC BAD", str); + end + end + printf("%s\n", str); + end + + % look for next packet + + bit = uw_loc + length(rtty.uw); + uw_loc = find_uw(rtty, bit, rx_bits_log, states.verbose); + + endwhile +endfunction + + +% Extract as many binary packets as we can from a great big buffer of bits, +% and send them to the C decoder for FEC decoding. +% horus_l2 can be compiled a bunch of different ways. You need to +% compile with: +% codec2-dev/src$ gcc horus_l2.c -o horus_l2 -Wall -DDEC_RX_BITS -DHORUS_L2_RX + +function corr_log = extract_and_decode_binary_packets(states, binary, rx_bits_log) + corr_log = []; + + % use UWs to delimit start and end of data packets + + bit = 1; + nbits = length(rx_bits_log); + + [uw_loc best_corr corr] = find_uw(binary, bit, rx_bits_log, states.verbose); + corr_log = [corr_log corr]; + + while (uw_loc != -1) + + if (uw_loc+binary.max_packet_len) < nbits + % printf("uw_loc: %d best_corr: %d\n", uw_loc, best_corr); + + % OK we have a packet delimited by two UWs. Lets convert the bit + % stream into bytes and save for decoding + + pin = uw_loc; + for i=1:45 + rx_bytes(i) = rx_bits_log(pin:pin+7) * (2.^(7:-1:0))'; + pin += 8; + %printf("%d 0x%02x\n", i, rx_bytes(i)); + end + + f=fopen("horus_rx_bits_binary.bin","wb"); + fwrite(f, rx_bytes, "uchar"); + fclose(f); + + % optionally write packet to disk to use as horus_tx_bits_binary.txt + f=fopen("horus_rx_bits_binary.txt","wt"); + for i=uw_loc:uw_loc+45*8-1 + fprintf(f, "%d ", rx_bits_log(i)); + end + fclose(f); + + system("../src/horus_l2"); % compile instructions above + end + + bit = uw_loc + length(binary.uw); + [uw_loc best_corr corr] = find_uw(binary, bit, rx_bits_log, states.verbose); + corr_log = [corr_log corr]; + + endwhile + +endfunction + + +% simulation of tx and rx side, add noise, channel impairments ---------------------- +% +% test_frame_mode Description +% 1 BER testing using known test frames +% 2 random bits +% 3 repeating sequence of all symbols +% 4 Horus RTTY +% 5 Horus Binary +% 6 Horus High Speed: A 8x oversampled modem, e.g. Fs=9600, Rs=1200 +% which is the same as Fs=921600 Rs=115200 +% Uses packet based BER counter + +function run_sim(test_frame_mode, M=2, frames = 10, EbNodB = 100, filename="fsk_horus.raw") + timing_offset = 0.0; % see resample() for clock offset below + fading = 0; % modulates tx power at 2Hz with 20dB fade depth, + % to simulate balloon rotating at end of mission + df = 0; % tx tone freq drift in Hz/s + dA = 1; % amplitude imbalance of tones (note this affects Eb so not a gd idea) + + more off + rand('state',1); + randn('state',1); + + % ---------------------------------------------------------------------- + + % sm2000 config ------------------------ + %states = fsk_horus_init(96000, 1200); + %states.f1_tx = 4000; + %states.f2_tx = 5200; + + if test_frame_mode < 4 + % horus rtty config --------------------- + states = fsk_horus_init(8000, 50, M); + end + + if test_frame_mode == 4 + % horus rtty config --------------------- + states = fsk_horus_init(8000, 100, 2); + states.tx_bits_file = "horus_payload_rtty.txt"; % Octave file of bits we FSK modulate + rtty = fsk_horus_init_rtty; + states.ntestframebits = rtty.max_packet_len; + end + + if test_frame_mode == 5 + % horus binary config --------------------- + states = fsk_horus_init(8000, 100, 4); + binary = fsk_horus_init_binary; + states.tx_bits_file = "horus_tx_bits_binary.txt"; % Octave file of bits we FSK modulate + states.ntestframebits = binary.max_packet_len; + end + + if test_frame_mode == 6 + % horus high speed --------------------- + states = fsk_horus_init(Fs=9600, Rs=1200, M=2, P=8, nsym=16); + states.tx_bits_file = "horus_high_speed.bin"; + end + + % Tones must be at least Rs apart for ideal non-coherent FSK + + states.ftx = 900 + 2*states.Rs*(1:states.M); + states.tx_tone_separation = states.ftx(2) - states.ftx(1); + + % ---------------------------------------------------------------------- + + states.verbose = 0x1; + M = states.M; + N = states.N; + P = states.P; + Rs = states.Rs; + nsym = states.nsym; + nbit = states.nbit; + Fs = states.Fs; + states.df(1:M) = df; + states.dA(1:M) = dA; + + % optional noise. Useful for testing performance of waveforms from real world modulators + + EbNo = 10^(EbNodB/10); + variance = states.Fs/(states.Rs*EbNo*states.bitspersymbol); + + % set up tx signal with payload bits based on test mode + + if (test_frame_mode == 1) + % test frame of bits, which we repeat for convenience when BER testing + states.ntestframebits = states.nbit; + test_frame = round(rand(1, states.ntestframebits)); + tx_bits = []; + for i=1:frames+1 + tx_bits = [tx_bits test_frame]; + end + end + + if test_frame_mode == 2 + % random bits, just to make sure sync algs work on random data + tx_bits = round(rand(1, states.nbit*(frames+1))); + end + + if test_frame_mode == 3 + % repeating sequence of all symbols + % great for initial test of demod if nothing else works, + % look for this pattern in rx_bits + if M == 2 + % ...10101... + tx_bits = zeros(1, states.nbit*(frames+1)); + tx_bits(1:2:length(tx_bits)) = 1; + else + % repeat each possible 4fsk symbol + pattern = [0 0 0 1 1 0 1 1]; + %pattern = [0 0 0 1 1 1 1 0]; + nrepeats = states.nbit*(frames+1)/length(pattern); + tx_bits = []; + for b=1:nrepeats + tx_bits = [tx_bits pattern]; + end + %tx_bits = zeros(1, states.nbit*(frames+1)); + end + end + + if (test_frame_mode == 4) || (test_frame_mode == 5) + + % load up a horus msg from disk and modulate that + + test_frame = load(states.tx_bits_file); + ltf = length(test_frame); + ntest_frames = ceil((frames+1)*nbit/ltf); + printf("Generating %d test packets\n", ntest_frames); + + % 1 second of random bits to let estimators lock on + + preamble = round(rand(1,states.Rs)); + + tx_bits = preamble; + for i=1:ntest_frames + tx_bits = [tx_bits test_frame]; + end + + % a packet len of random bits at end fill buffers to deocode final packet + + if test_frame_mode == 4 + postamble = round(rand(1,rtty.max_packet_len)); + else + postamble = round(rand(1,binary.max_packet_len)); + end + tx_bits = [tx_bits postamble]; + end + + if test_frame_mode == 6 + states.verbose += 0x4; + ftmp = fopen(states.tx_bits_file, "rb"); test_frame = fread(ftmp,Inf,"char")'; fclose(ftmp); + states.ntestframebits = length(test_frame); + printf("length test frame: %d\n", states.ntestframebits); + %test_frame = rand(1,states.ntestframebits) > 0.5; + + tx_bits = []; + for i=1:frames+1 + tx_bits = [tx_bits test_frame]; + end + end + + tx = fsk_mod(states, tx_bits); + + %tx = resample(tx, 1000, 1001); % simulated 1000ppm sample clock offset + + if fading + ltx = length(tx); + tx = tx .* (1.1 + cos(2*pi*2*(0:ltx-1)/Fs))'; % min amplitude 0.1, -20dB fade, max 3dB + end + + noise = sqrt(variance)*randn(length(tx),1); + rx = tx + noise; + printf("SNRdB meas: %4.1f\n", 10*log10(var(tx)/var(noise))); + + % dump simulated rx file + + ftx=fopen(filename,"wb"); rxg = rx*1000; fwrite(ftx, rxg, "short"); fclose(ftx); + + timing_offset_samples = round(timing_offset*states.Ts); + st = 1 + timing_offset_samples; + rx_bits_buf = zeros(1,nbit+states.ntestframebits); + x_log = []; + timing_nl_log = []; + norm_rx_timing_log = []; + f_int_resample_log = []; + f_log = []; + EbNodB_log = []; + rx_bits_log = []; + rx_bits_sd_log = []; + + % main loop --------------------------------------------------------------- + + run_frames = floor(length(rx)/N)-1; + for f=1:run_frames + + % extract nin samples from input stream + + nin = states.nin; + en = st + states.nin - 1; + + if en < length(rx) % due to nin variations its possible to overrun buffer + sf = rx(st:en); + st += nin; + + % demodulate to stream of bits + + states = est_freq(states, sf, states.M); + %states.f = 900 + 2*states.Rs*(1:states.M); + %states.f = [1200 1400 1600 1800]; + [rx_bits states] = fsk_demod(states, sf); + + rx_bits_buf(1:states.ntestframebits) = rx_bits_buf(nbit+1:states.ntestframebits+nbit); + rx_bits_buf(states.ntestframebits+1:states.ntestframebits+nbit) = rx_bits; + %rx_bits_buf(1:nbit) = rx_bits_buf(nbit+1:2*nbit); + %rx_bits_buf(nbit+1:2*nbit) = rx_bits; + + rx_bits_log = [rx_bits_log rx_bits]; + rx_bits_sd_log = [rx_bits_sd_log states.rx_bits_sd]; + + norm_rx_timing_log = [norm_rx_timing_log states.norm_rx_timing]; + x_log = [x_log states.x]; + timing_nl_log = [timing_nl_log states.timing_nl]; + f_int_resample_log = [f_int_resample_log abs(states.f_int_resample(:,:))]; + f_log = [f_log; states.f]; + EbNodB_log = [EbNodB_log states.EbNodB]; + + if test_frame_mode == 1 + states = ber_counter(states, test_frame, rx_bits_buf); + end + if test_frame_mode == 6 + states = ber_counter_packet(states, test_frame, rx_bits_buf); + end + end + end + + % print stats, count errors, decode packets ------------------------------------------ + + if (test_frame_mode == 1) || (test_frame_mode == 6) + printf("frames: %d EbNo: %3.2f Tbits: %d Terrs: %d BER %4.3f\n", frames, EbNodB, states.Tbits,states. Terrs, states.Terrs/states.Tbits); + end + + if test_frame_mode == 4 + npackets = extract_and_print_rtty_packets(states, rtty, rx_bits_log, rx_bits_sd_log); + printf("Received %d packets\n", npackets); + end + + if test_frame_mode == 5 + extract_and_decode_binary_packets(states, binary, rx_bits_log); + end + + figure(1); + plot(f_int_resample_log','+') + hold off; + + figure(2) + clf + m = max(abs(x_log)); + plot(x_log,'+') + axis([-m m -m m]) + title('fine timing metric') + + figure(3) + clf + subplot(211) + plot(norm_rx_timing_log); + axis([1 run_frames -1 1]) + title('norm fine timing') + subplot(212) + plot(states.nerr_log) + title('num bit errors each frame') + + figure(4) + clf + subplot(211) + one_sec_rx = rx(1:min(Fs,length(rx))); + plot(one_sec_rx) + title('rx signal at demod input') + subplot(212) + plot(abs(fft(one_sec_rx))) + + figure(5) + clf + plot(f_log,'+') + title('tone frequencies') + axis([1 run_frames 0 Fs/2]) + + figure(6) + clf + plot(EbNodB_log); + title('Eb/No estimate') + + figure(7) + clf + subplot(211) + X = abs(fft(timing_nl_log)); + plot(X(1:length(X)/2)) + subplot(212) + plot(abs(timing_nl_log(1:100))) + + endfunction + + +% --------------------------------------------------------------------- +% demodulate from a user-supplied file +% --------------------------------------------------------------------- + +function rx_bits_log = demod_file(filename, test_frame_mode=4, noplot=0, EbNodB=100, max_frames=1E32) + fin = fopen(filename,"rb"); + more off; + read_complex = 0; sample_size = 'int16'; shift_fs_on_4 = 0; + max_frames + + if test_frame_mode == 4 + % horus rtty config --------------------- + states = fsk_horus_init(8000, 100, 2); + rtty = fsk_horus_init_rtty; + states.ntestframebits = rtty.max_packet_len; + end + + if test_frame_mode == 5 + % horus binary config --------------------- + states = fsk_horus_init(8000, 100, 4); + binary = fsk_horus_init_binary; + states.ntestframebits = binary.max_packet_len; + end + + states.verbose = 0x1 + 0x8; + + if test_frame_mode == 6 + % Horus high speed config -------------- + states = fsk_horus_init(Fs=9600, Rs=1200, M=2, P=8, nsym=16); + states.tx_bits_file = "horus_high_speed.bin"; + states.verbose += 0x4; + ftmp = fopen(states.tx_bits_file, "rb"); test_frame = fread(ftmp,Inf,"char")'; fclose(ftmp); + states.ntestframebits = length(test_frame); + printf("length test frame: %d\n", states.ntestframebits); + end + + if test_frame_mode == 7 + % 800XA 4FSK modem -------------- + states = fsk_init(Fs=8000, Rs=400, M=4, P=10, nsym=256); + states.tx_bits_file = "horus_high_speed.bin"; + states.verbose += 0x4; + ftmp = fopen(states.tx_bits_file, "rb"); test_frame = fread(ftmp,Inf,"char")'; fclose(ftmp); + states.ntestframebits = length(test_frame); + printf("length test frame: %d\n", states.ntestframebits); + end + + if test_frame_mode == 8 + % test RS41 type balllon telemetry -------------- + states = fsk_init(96000, 4800, 2, 10, 16); + states.fest_fmin = 1000; + states.fest_fmax = 40000; + states.fest_min_spacing = 1000; + states.tx_bits_file = "../build_linux/src/tx_bit.bin"; + states.verbose += 0x4; + #ftmp = fopen(states.tx_bits_file, "rb"); test_frame = fread(ftmp,Inf,"char")'; fclose(ftmp); + #states.ntestframebits = length(test_frame); + #printf("length test frame: %d\n", states.ntestframebits); + states.ntestframebits = 1000; + read_complex = 1; + shift_fs_on_4 = 1; % get samples into range of current freq estimator + end + + if test_frame_mode == 9 + % Wenet high speed SSTV, we can just check raw demo here --------------------- + % despite the high sample rate the modem sees this as a 8:1 Fs/Rs configuration + states = fsk_init(8000, 1000, 2); + states.tx_tone_separation = 1000; + states.ntestframebits = (256+2+65)*8+40; % from src/drs232_lpc.c + states.freq_est_type = 'mask'; + read_complex=1; sample_size = 'uint8'; + printf("Wenet mode: ntestframebits: %d freq_est_type: %s\n", states.ntestframebits, states.freq_est_type); + %states.verbose = 0x8; + end + + N = states.N; + P = states.P; + Rs = states.Rs; + nsym = states.nsym; + nbit = states.nbit; + + frames = 0; + rx = []; + rx_bits_log = []; + rx_bits_sd_log = []; + norm_rx_timing_log = []; + f_int_resample_log = []; + EbNodB_log = []; + ppm_log = []; + f_log = []; Sf_log = []; + + rx_bits_buf = zeros(1,nbit + states.ntestframebits); + + % optional noise. Useful for testing performance of waveforms from real world modulators + % we need to pre-read the file to estimate the signal power + ftmp = fopen(filename,"rb"); s = fread(ftmp,Inf,sample_size); fclose(ftmp); + if sample_size == "uint8" s = (s - 127)/128; end + if read_complex s = s(1:2:end) + j*s(2:2:end); end + tx_pwr = var(s); + EbNo = 10^(EbNodB/10); + variance = (tx_pwr/2)*states.Fs/(states.Rs*EbNo*states.bitspersymbol); + + % First extract raw bits from samples ------------------------------------------------------ + + printf("demod of raw bits....\n"); + + finished = 0; ph = 1; + while (finished == 0) + + % extract nin samples from input stream + + nin = states.nin; + if read_complex + [sf count] = fread(fin, 2*nin, sample_size); + if sample_size == "uint8" sf = (sf - 127)/128; end + sf = sf(1:2:end) + j*sf(2:2:end); + count /= 2; + if shift_fs_on_4 + % optional shift up in freq by Fs/4 to get into freq est range + for i=1:count + ph = ph*exp(j*pi/4); + sf(i) *= ph; + end + end + else + [sf count] = fread(fin, nin, "short"); + end + rx = [rx; sf]; + + % add optional noise + + if count + noise = sqrt(variance)*randn(count,1); + sf += noise; + end + + if count == nin + frames++; + + % demodulate to stream of bits + + states = est_freq(states, sf, states.M); + if states.freq_est_type == 'mask' states.f = states.f2; end + [rx_bits states] = fsk_demod(states, sf); + + rx_bits_buf(1:states.ntestframebits) = rx_bits_buf(nbit+1:states.ntestframebits+nbit); + rx_bits_buf(states.ntestframebits+1:states.ntestframebits+nbit) = rx_bits; + + rx_bits_log = [rx_bits_log rx_bits]; + rx_bits_sd_log = [rx_bits_sd_log states.rx_bits_sd]; + norm_rx_timing_log = [norm_rx_timing_log states.norm_rx_timing]; + f_int_resample_log = [f_int_resample_log abs(states.f_int_resample)]; + EbNodB_log = [EbNodB_log states.EbNodB]; + ppm_log = [ppm_log states.ppm]; + f_log = [f_log; states.f]; + Sf_log = [Sf_log; states.Sf']; + + if (test_frame_mode == 1) + states = ber_counter(states, test_frame, rx_bits_buf); + if states.ber_state == 1 + states.verbose = 0; + end + end + if (test_frame_mode == 6) % || (test_frame_mode == 8) + states = ber_counter_packet(states, test_frame, rx_bits_buf); + end + else + finished = 1; + end + + if frames > max_frames finished=1; end + + end + printf("frames: %d\n", frames); + fclose(fin); + + if noplot == 0 + printf("plotting...\n"); + + figure(1); clf; + plot(f_log); + title('Tone Freq Estimates'); + + figure(2); + plot(f_int_resample_log','+') + title('Integrator outputs for each tone'); + + figure(3); clf + subplot(211) + plot(norm_rx_timing_log) + axis([1 frames -0.5 0.5]) + title('norm fine timing') + subplot(212) + plot(states.nerr_log) + title('num bit errors each frame') + + figure(4); clf + plot(EbNodB_log); + title('Eb/No estimate') + + figure(5); clf + rx_nowave = rx(1000:length(rx)); % skip past wav header if it's a wave file + subplot(211) + plot(real(rx_nowave)); + title('input signal to demod (1 sec)') + xlabel('Time (samples)'); + %axis([1 states.Fs -35000 35000]) + + % normalise spectrum to 0dB full scale with sine wave input + subplot(212); + if sample_size == "int16" max_value = 32767; end + if sample_size == "uint8" max_value = 127; end + RxdBFS = 20*log10(abs(fft(rx_nowave(1:states.Fs)))) - 20*log10((states.Fs/2)*max_value); + plot(RxdBFS) + axis([1 states.Fs/2 -80 0]) + xlabel('Frequency (Hz)'); + + figure(6); clf + plot(ppm_log) + title('Sample clock (baud rate) offset in PPM'); + + figure(7); clf; mesh(Sf_log(1:10,:)); + end + + if (test_frame_mode == 1) || (test_frame_mode == 6) + printf("frames: %d Tbits: %d Terrs: %d BER %4.3f EbNo: %3.2f\n", frames, states.Tbits,states. Terrs, states.Terrs/states.Tbits, mean(EbNodB_log)); + end + + % we can decode both protocols at the same time + + if (test_frame_mode == 4) + npackets = extract_and_print_rtty_packets(states, rtty, rx_bits_log, rx_bits_sd_log) + printf("Received %d packets\n", npackets); + end + + if (test_frame_mode == 5) + corr_log = extract_and_decode_binary_packets(states, binary, rx_bits_log); + + figure(8); + clf + plot(corr_log); + hold on; + plot([1 length(corr_log)],[binary.uw_thresh binary.uw_thresh],'g'); + hold off; + title('UW correlation'); + end + +endfunction + + +% Over the years this modem has been used for many different FSK signals ... + +if exist("fsk_horus_as_a_lib") == 0 + run_sim(test_frame_mode=4, M=2, frames=30, EbNodB = 20); + %run_sim(5, 4, 30, 100); + %rx_bits = demod_file("~/Desktop/115.wav",6,0,90); + %rx_bits = demod_file("~/Desktop/fsk_800xa_rx_hackrf.wav",7); + %rx_bits = demod_file("~/Desktop/2fsk_100_rx_rpi_rtlsdr_002_ledger.wav",4); + %rx_bits = demod_file("~/Desktop/phorus_binary_ascii.wav",4); + %rx_bits = demod_file("~/Desktop/binary/horus_160102_binary_rtty_2.wav",4); + %rx_bits = demod_file("~/Desktop/horus_160102_vk5ei_capture2.wav",4); + %rx_bits = demod_file("~/Desktop/horus_rtty_binary.wav",4); + %rx_bits = demod_file("~/Desktop/FSK_4FSK.wav",4); + %rx_bits = demod_file("t.raw",5); + %rx_bits = demod_file("~/Desktop/fsk_horus_10dB_1000ppm.wav",4); + %rx_bits = demod_file("~/Desktop/fsk_horus_6dB_0ppm.wav",4); + %rx_bits = demod_file("test.raw",1,1); + %rx_bits = .rawdemod_file("/dev/ttyACM0",1); + %rx_bits = demod_file("fsk_horus_rx_1200_96k.raw",1); + %rx_bits = demod_file("mp.raw",4); + %rx_bits = demod_file("~/Desktop/launchbox_v2_landing_8KHz_final.wav",4); + %rx_bits = demod_file("~/Desktop/fsk_800xa.wav",7); + %rx_bits = demod_file("~/Desktop/rs41_96k_10s.iq16",8); +end |