From: drowe67 Date: Fri, 2 Jun 2017 01:19:22 +0000 (+0000) Subject: made timing/freq est window indep of Rs X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=1dbe3f5b062a27a0e343848bfc3db2578453ef4b;p=freetel-svn-tracking.git made timing/freq est window indep of Rs git-svn-id: https://svn.code.sf.net/p/freetel/code@3150 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/fsk_eme.m b/codec2-dev/octave/fsk_eme.m new file mode 100644 index 00000000..47e57a23 --- /dev/null +++ b/codec2-dev/octave/fsk_eme.m @@ -0,0 +1,377 @@ +% fsk_eme.m +% David Rowe July 2017 +% +% FSK modem for digital voice over 1296 MHz Earth-Moon-Earth (EME) links + +fsk_lib; + +function run_sim(frames = 10, EbNodB = 100) + Fs = 8000; Rs = 8; M = 16; Fsep = 50; + + 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); + + states = fsk_init(Fs, 50, M); + + states.ftx = 1000 + Fsep*(-M/2:M/2-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); + + % 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 + + + 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 + + g = 1000; + ftx=fopen("fsk_horus_tx.raw","wb"); txg = tx*g; fwrite(ftx, txg, "short"); fclose(ftx); + frx=fopen("fsk_horus_rx.raw","wb"); rxg = rx*g; fwrite(frx, rxg, "short"); fclose(frx); + + 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 = states.ftx; + [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]; + + states = ber_counter(states, test_frame, rx_bits_buf); + end + end + + % print stats, count errors, decode packets ------------------------------------------ + + printf("frames: %d EbNo: %3.2f Tbits: %d Terrs: %d BER %4.3f\n", frames, EbNodB, states.Tbits,states. Terrs, states.Terrs/states.Tbits); + + 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 a file of 8kHz 16bit short samples -------------------------------- + +function rx_bits_log = demod_file(filename, test_frame_mode, noplot=0, EbNodB=100) + fin = fopen(filename,"rb"); + more off; + + %states = fsk_horus_init(96000, 1200); + + if test_frame_mode == 4 + % horus rtty config --------------------- + states = fsk_horus_init(8000, 100, 2); + uwstates = fsk_horus_init_rtty_uw(states); + states.ntestframebits = states.nbits; + end + + if test_frame_mode == 5 + % horus binary config --------------------- + states = fsk_horus_init(8000, 50, 4); + uwstates = fsk_horus_init_binary_uw; + states.ntestframebits = states.nbits; + end + + states.verbose = 0x1 + 0x8; + + if test_frame_mode == 6 + % Horus high speed config -------------- + states = fsk_horus_init_hbr(9600, 8, 1200, 2, 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_horus_init_hbr(8000, 10, 400, 4, 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 + + 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 = []; + rx_bits_buf = zeros(1,nbit + states.ntestframebits); + + % optional noise. Useful for testing performance of waveforms from real world modulators + + EbNo = 10^(EbNodB/10); + ftmp = fopen(filename,"rb"); s = fread(ftmp,Inf,"short"); fclose(ftmp); tx_pwr = var(s); + 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; + while (finished == 0) + + % extract nin samples from input stream + + nin = states.nin; + [sf count] = fread(fin, nin, "short"); + 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); + %states.f = [1450 1590 1710 1850]; + [rx_bits states] = fsk_horus_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]; + + 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 + states = ber_counter_packet(states, test_frame, rx_bits_buf); + end + else + finished = 1; + end + end + fclose(fin); + + if noplot == 0 + printf("plotting...\n"); + + figure(1); + plot(f_log); + hold off; + + figure(2); + plot(f_int_resample_log','+') + + figure(3) + clf + subplot(211) + plot(norm_rx_timing_log) + axis([1 frames -0.5 0.5]) + title('norm fine timing') + grid + 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)); + subplot(211) + plot(rx_nowave(1:states.Fs)); + title('input signal to demod (1 sec)') + xlabel('Time (samples)'); + axis([1 states.Fs -35000 35000]) + + % normalise spectrum to 0dB full scale with a 32767 sine wave input + + subplot(212) + RxdBFS = 20*log10(abs(fft(rx_nowave(1:states.Fs)))) - 20*log10((states.Fs/2)*32767); + 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'); + 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) || (test_frame_mode == 5) + extract_and_print_rtty_packets(states, rx_bits_log, rx_bits_sd_log) + corr_log = extract_and_decode_binary_packets(states, rx_bits_log); + + figure(8); + clf + plot(corr_log); + hold on; + plot([1 length(corr_log)],[states.binary.uw_thresh states.binary.uw_thresh],'g'); + hold off; + title('UW correlation'); + end + +endfunction + + +% Start simulations here ------------------------------------------------------------- + +more off; format; + +%run_sim(1, 2, 100, 9); +run_sim(10, 20); diff --git a/codec2-dev/octave/fsk_horus.m b/codec2-dev/octave/fsk_horus.m index 3f1ed936..eb44fda9 100644 --- a/codec2-dev/octave/fsk_horus.m +++ b/codec2-dev/octave/fsk_horus.m @@ -298,8 +298,6 @@ function corr_log = extract_and_decode_binary_packets(states, rx_bits_log) endfunction -% Alternative stateless BER counter that works on packets that may have gaps between them - % simulation of tx and rx side, add noise, channel impairments ---------------------- % % test_frame_mode Description @@ -332,7 +330,7 @@ function run_sim(test_frame_mode, M=2, frames = 10, EbNodB = 100) if test_frame_mode < 4 % horus rtty config --------------------- - states = fsk_horus_init(8000, 8, M); + states = fsk_horus_init(8000, 50, M); %states = fsk_horus_init_hbr(8000, 10, 400, 4); % EME end @@ -784,7 +782,7 @@ endfunction % run test functions from here during development if exist("fsk_horus_as_a_lib") == 0 - run_sim(1, 16, 100, 100); + run_sim(1, 2, 100, 9); %rx_bits = demod_file("~/Desktop/115.wav",6,0,90); %rx_bits = demod_file("fsk_horus.raw",5); %rx_bits = demod_file("~/Desktop/4FSK_Binary_NoLock.wav",4); diff --git a/codec2-dev/octave/fsk_lib.m b/codec2-dev/octave/fsk_lib.m index 46c68db2..87252428 100644 --- a/codec2-dev/octave/fsk_lib.m +++ b/codec2-dev/octave/fsk_lib.m @@ -19,25 +19,26 @@ function states = fsk_init(Fs, Rs, M=2) states.M = M; states.bitspersymbol = log2(M); states.Fs = Fs; - N = states.N = Fs; % processing buffer size, nice big window for timing est - %states.Ndft = 2.^ceil(log2(N)); % find nearest power of 2 for efficient FFT - states.Ndft = 1024; % find nearest power of 2 for efficient FFT states.Rs = Rs; - Ts = states.Ts = Fs/Rs; + + states.nsym = 50; % need enough symbols for good timing and freq offset est + Ts = states.Ts = Fs/Rs; % number of samples per symbol assert(Ts == floor(Ts), "Fs/Rs must be an integer"); - states.nsym = N/Ts; % number of symbols in one processing frame + + N = states.N = Ts*states.nsym; % processing buffer size, nice big window for timing est + states.Ndft = min(1024, 2.^ceil(log2(N))); % find nearest power of 2 for efficient FFT states.nbit = states.nsym*states.bitspersymbol; % number of bits per processing frame - Nmem = states.Nmem = N+2*Ts; % two symbol memory in down converted signals to allow for timing adj + Nmem = states.Nmem = N+2*Ts; % two symbol memory in down converted signals to allow for timing adj states.Sf = zeros(states.Ndft/2,1); % current memory of dft mag samples states.f_dc = zeros(M,Nmem); - states.P = 8; % oversample rate out of filter + states.P = 8; % oversample rate out of filter assert(Ts/states.P == floor(Ts/states.P), "Ts/P must be an integer"); - states.nin = N; % can be N +/- Ts/P samples to adjust for sample clock offsets + states.nin = N; % can be N +/- Ts/P samples to adjust for sample clock offsets states.verbose = 0; - states.phi = zeros(1, M); % keep down converter osc phase continuous + states.phi = zeros(1, M); % keep down converter osc phase continuous %printf("M: %d Fs: %d Rs: %d Ts: %d nsym: %d nbit: %d\n", states.M, states.Fs, states.Rs, states.Ts, states.nsym, states.nbit); @@ -387,7 +388,7 @@ function [rx_bits states] = fsk_demod(states, sf) states.f_int_resample = f_int_resample; states.rx_bits_sd = rx_bits_sd; - % Eb/No estimation (todo: this needs some work, like calibrartion, low Eb/No perf) + % Eb/No estimation (todo: this needs some work, like calibration, low Eb/No perf) tone_max = abs(tone_max); states.EbNodB = -6 + 20*log10(1E-6+mean(tone_max)/(1E-6+std(tone_max)));