From 9807e9663cf626b5451c2910c97b8c44e78489ec Mon Sep 17 00:00:00 2001 From: drowe67 Date: Mon, 13 Feb 2017 07:01:10 +0000 Subject: [PATCH] all curves up to HF working, about to change from 1 sample per symbol to M samples/symbol for diversity git-svn-id: https://svn.code.sf.net/p/freetel/code@3030 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/hf_modem_curves.m | 140 +++++++++++++++++++++++++--- 1 file changed, 125 insertions(+), 15 deletions(-) diff --git a/codec2-dev/octave/hf_modem_curves.m b/codec2-dev/octave/hf_modem_curves.m index aaf43dfe..5cc68216 100644 --- a/codec2-dev/octave/hf_modem_curves.m +++ b/codec2-dev/octave/hf_modem_curves.m @@ -5,11 +5,14 @@ % generate plots for a blog post. #{ - [ ] ideal AWGN/HF curves - [ ] exp AWGN QPSK curves - [ ] exp AWGN DQPSK curves - [ ] exp HF channel model - [ ] with COHPSK frames + [X] ideal AWGN/HF curves + [X] exp AWGN QPSK curves + [X] exp AWGN DQPSK curves + [X] exp HF channel model + [ ] diversity + [ ] COHPSK frames + + would require multiple carriers + + filtering or OFDM #} 1; @@ -36,70 +39,176 @@ function two_bits = qpsk_demod(symbol) endfunction +% Rate Rs modem simulation model ------------------------------------------------------- + function sim_out = ber_test(sim_in) bps = 2; % two bits/symbol for QPSK verbose = sim_in.verbose; - nbits = sim_in.nbits; - nsymb = nbits/bps; EbNovec = sim_in.EbNovec; + hf_en = sim_in.hf_en; + + % user can supply number of bits per point to get good results + % at high Eb/No + + if length(sim_in.nbits) > 1 + nbitsvec = sim_in.nbits; + nbitsvec += 100 - mod(nbitsvec,100); % round up to nearest 100 + else + nbitsvec(1:length(EbNovec)) = sim_in.nbits; + end + + % init HF model + + if hf_en + % some typical values + Rs = 50; dopplerSpreadHz = 1; path_delay = 0.001/Rs; + + nsymb = max(nbitsvec)/2; + spread1 = doppler_spread(dopplerSpreadHz, Rs, nsymb); + spread2 = doppler_spread(dopplerSpreadHz, Rs, nsymb); + hf_gain = 1.0/sqrt(var(spread1)+var(spread2)); + %printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1)); + end for ne = 1:length(EbNovec) + + % work out noise power ------------- + EbNodB = EbNovec(ne); EsNodB = EbNodB + 10*log10(bps); EsNo = 10^(EsNodB/10); - variance = 1/EsNo; + nbits = nbitsvec(ne); + nsymb = nbits/bps; + + % modulator ------------------------ tx_bits = rand(1,nbits) > 0.5; tx_symb = []; + prev_tx_symb = 1; for s=1:nsymb atx_symb = qpsk_mod(tx_bits(2*s-1:2*s)); + if sim_in.dqpsk + atx_symb *= prev_tx_symb; + prev_tx_symb = atx_symb; + end tx_symb = [tx_symb atx_symb]; end + % channel --------------------------- + + rx_symb = tx_symb; + + if hf_en + + % simplified rate Rs simulation model that doesn't include + % ISI, just freq filtering. We assume perfect phase estimation + % so it's just amplitude distortion. + + for s=1:nsymb + hf_model = hf_gain*(spread1(s) + exp(-j*path_delay)*spread2(s)); + rx_symb(s) = rx_symb(s)*abs(hf_model); + end + end + % variance is noise power, which is divided equally between real and % imag components of noise noise = sqrt(variance*0.5)*(randn(1,nsymb) + j*randn(1,nsymb)); - rx_symb = tx_symb + noise; + rx_symb += noise; + + % demodulator ------------------------------------------ rx_bits = []; + prev_rx_symb = 1; for s=1:nsymb - two_bits = qpsk_demod(rx_symb(s)); + arx_symb = rx_symb(s); + if sim_in.dqpsk + tmp = arx_symb; + arx_symb *= prev_rx_symb'; + prev_rx_symb = tmp; + end + two_bits = qpsk_demod(arx_symb); rx_bits = [rx_bits two_bits]; end + % count errors ----------------------------------------- + error_pattern = xor(tx_bits, rx_bits); nerrors = sum(error_pattern); bervec(ne) = nerrors/nbits; if verbose printf("EbNodB: %3.1f nerrors: %5d ber: %4.3f\n", EbNodB, nerrors, bervec(ne)); + if verbose == 2 + figure(2); clf; + plot(rx_symb,'+'); + end end end sim_out.bervec = bervec; endfunction + % ------------------------------------------------------------- + +function run_single + sim_in.verbose = 2; + sim_in.nbits = 10000; + sim_in.EbNovec = 8; + sim_in.dqpsk = 0; + sim_in.hf_en = 1; + + sim_qpsk = ber_test(sim_in); +endfunction + + function run_curves sim_in.verbose = 1; - sim_in.nbits = 100000; - sim_in.EbNovec = 0:7; + sim_in.EbNovec = 0:10; + sim_in.dqpsk = 0; + sim_in.hf_en = 0; + + % AWGN ----------------------------- + + ber_awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNovec/10))); + sim_in.nbits = min(1E5, floor(500 ./ ber_awgn_theory)); - ber_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNovec/10))); sim_qpsk = ber_test(sim_in); + sim_in.dqpsk = 1; + sim_dqpsk = ber_test(sim_in); + + % HF ----------------------------- + + hf_sim_in = sim_in; hf_sim_in.dqpsk = 0; hf_sim_in.hf_en = 1; + hf_sim_in.EbNovec = 0:16; + + EbNoLin = 10.^(hf_sim_in.EbNovec/10); + ber_hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); + + hf_sim_in.nbits = min(1E5, floor(500 ./ ber_hf_theory)); + sim_qpsk_hf = ber_test(hf_sim_in); + + hf_sim_in.dqpsk = 1; + sim_dqpsk_hf = ber_test(hf_sim_in); + + % Plot results -------------------- figure(1); clf; - semilogy(sim_in.EbNovec, ber_theory,'r+-;QPSK AWGN theory;','markersize', 10, 'linewidth', 2) + semilogy(sim_in.EbNovec, ber_awgn_theory,'r+-;QPSK AWGN theory;','markersize', 10, 'linewidth', 2) hold on; + semilogy(hf_sim_in.EbNovec, ber_hf_theory,'r+-;QPSK HF theory;','markersize', 10, 'linewidth', 2) semilogy(sim_in.EbNovec, sim_qpsk.bervec,'g+-;QPSK AWGN simulated;','markersize', 10, 'linewidth', 2) + semilogy(sim_in.EbNovec, sim_dqpsk.bervec,'b+-;DQPSK AWGN simulated;','markersize', 10, 'linewidth', 2) + semilogy(hf_sim_in.EbNovec, sim_qpsk_hf.bervec,'c+-;QPSK HF simulated;','markersize', 10, 'linewidth', 2) + semilogy(hf_sim_in.EbNovec, sim_dqpsk_hf.bervec,'k+-;DQPSK HF simulated;','markersize', 10, 'linewidth', 2) hold off; xlabel('Eb/N0') ylabel('BER') grid("minor") - axis([min(sim_in.EbNovec) max(sim_in.EbNovec) 1E-3 1]) + axis([min(hf_sim_in.EbNovec) max(hf_sim_in.EbNovec) 1E-3 1]) endfunction @@ -108,3 +217,4 @@ endfunction more off; rand('seed',1); randn('seed', 1); run_curves +%run_single -- 2.25.1