--- /dev/null
+% 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);
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);
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)));
projects / freetel-svn-tracking.git / commitdiff