From 315306d66629f452e45e7171c42e19ea8a6db1ac Mon Sep 17 00:00:00 2001 From: drowe67 Date: Tue, 7 Oct 2014 09:50:04 +0000 Subject: [PATCH] maximum likelihood dqpsk, works well for AWGN, not so good for HF git-svn-id: https://svn.code.sf.net/p/freetel/code@1879 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fdmdv.m | 2 +- codec2-dev/octave/ldpcut.m | 2 +- .../octave/{test_ucinter.m => test_ml.m} | 290 +++++++++--------- codec2-dev/octave/test_qpsk2.m | 24 +- 4 files changed, 154 insertions(+), 164 deletions(-) rename codec2-dev/octave/{test_ucinter.m => test_ml.m} (50%) diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m index c6aaf866..415a6d7d 100644 --- a/codec2-dev/octave/fdmdv.m +++ b/codec2-dev/octave/fdmdv.m @@ -1147,7 +1147,7 @@ for c=1:Nc/2 freq_pol(c) = 2*pi*carrier_freq/Fs; freq(c) = exp(j*freq_pol(c)); end -for c=Nc/2+1:Nc +for c=floor(Nc/2)+1:Nc %carrier_freq = (-Nc/2 + c)*Fsep + Fcentre; carrier_freq = (-Nc/2 + c)*Fsep; freq_pol(c) = 2*pi*carrier_freq/Fs; diff --git a/codec2-dev/octave/ldpcut.m b/codec2-dev/octave/ldpcut.m index 80142ba6..9dfdc94e 100644 --- a/codec2-dev/octave/ldpcut.m +++ b/codec2-dev/octave/ldpcut.m @@ -53,7 +53,7 @@ end for nn = 1: Ntrials st = (nn-1)*code_param.symbols_per_frame + 1; en = (nn)*code_param.symbols_per_frame; - detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r(st:en), EsNo); + detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r(st:en), EsNo, ones(1,code_param.symbols_per_frame)); st = (nn-1)*code_param.data_bits_per_frame + 1; en = (nn)*code_param.data_bits_per_frame; error_positions = xor( detected_data(1:code_param.data_bits_per_frame), data(st:en) ); diff --git a/codec2-dev/octave/test_ucinter.m b/codec2-dev/octave/test_ml.m similarity index 50% rename from codec2-dev/octave/test_ucinter.m rename to codec2-dev/octave/test_ml.m index 365bca48..51dd8cb6 100644 --- a/codec2-dev/octave/test_ucinter.m +++ b/codec2-dev/octave/test_ml.m @@ -1,27 +1,10 @@ -% test_ucinter.m -% David Rowe August 2014 +% test_ml.m +% David Rowe Oct 2014 % -% FDM QPSK modem simulation to test uncododed interleaving ideas on HF -% channels without building a full blown modem. - -% [X] baseline QPSK simulation AWGN -% [ ] Spreading function -% [ ] visualise different carriers with and without spreading -% [ ] prove perf same as AWGN when one carrier is knocked out -% + AWGN and "faded channel" with same average SNR -% + use contrived example -% + then try less contrived but still well behaived maths channel - -% Ideas: -% + decode quickly but then slow down while playing, using full interleave -% + exp type window so we can decode using current symbols alone in gd SNR, or extend window over greater time -% + like root nyq filtering, combining multiple symbols, weighted -% + spreading could make it worse, e.g. short term average might be very low -% + SSB spreads out over a long way. We could so this too! Like spread spectrum. Don't have to be related to -% carrier width. Minimise amount of stuff that gets nailed. But didn't we try this with "1600 wide"? OK, so - key issue was BER. That, is what we need to impove, using "high areas". - +% Simulation to test FDM QPSK with maximum likelihood decoding on +% fading channels. + 1; % main test function @@ -30,27 +13,22 @@ function sim_out = ber_test(sim_in, modulation) Fs = 8000; verbose = sim_in.verbose; - framesize = sim_in.framesize; Ntrials = sim_in.Ntrials; Esvec = sim_in.Esvec; - phase_offset = sim_in.phase_offset; - w_offset = sim_in.w_offset; plot_scatter = sim_in.plot_scatter; Rs = sim_in.Rs; hf_sim = sim_in.hf_sim; nhfdelay = sim_in.hf_delay_ms*Rs/1000; - hf_phase_only = sim_in.hf_phase_only; hf_mag_only = sim_in.hf_mag_only; - Nc = sim_in.Nc; + framesize = sim_in.framesize; + ml = sim_in.ml; bps = 2; - Nsymb = framesize/bps; - prev_sym_tx = qpsk_mod([0 0]); - prev_sym_rx = qpsk_mod([0 0]); + Nc = Nsymb = framesize/bps; % total number of symbols - tx_bits_buf = zeros(1,2*framesize); - rx_bits_buf = zeros(1,2*framesize); - rx_symb_buf = zeros(1,2*Nsymb); + prev_sym_tx = qpsk_mod([0 0])*ones(1,Nc); + prev_sym_rx = qpsk_mod([0 0])*ones(1,Nc); + r = qpsk_mod([0 0])*ones(Nc,4); % Init HF channel model from stored sample files of spreading signal ---------------------------------- @@ -74,7 +52,7 @@ function sim_out = ber_test(sim_in, modulation) b = fir1(50, 5/Fs); spread = filter(b,1,spreadbb); spread_2ms = filter(b,1,spreadbb_2ms); - + % discard first 1000 samples as these were near 0, probably as % PathSim states were ramping up @@ -108,125 +86,130 @@ function sim_out = ber_test(sim_in, modulation) tx_symb_log = []; rx_symb_log = []; - noise_log = []; + errors_log = []; + Nerrs_log = []; % init HF channel hf_n = 1; - hf_angle_log = []; - hf_fading = ones(1,Nsymb); % default input for ldpc dec - hf_model = ones(Ntrials*Nsymb/Nc, Nc); % defaults for plotting surface + % simulation starts here----------------------------------- + for nn = 1: Ntrials - tx_bits = round( rand( 1, framesize) ); - + tx_bits = round(rand(1,framesize)); + % modulate -------------------------------------------- - s = zeros(1, Nsymb); - for i=1:Nsymb - tx_symb = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); - if strcmp(modulation,'dqpsk') - tx_symb *= prev_sym_tx; - prev_sym_tx = tx_symb; - end - s(i) = tx_symb; + for c=1:Nc + tx_symb(c) = qpsk_mod(tx_bits((2*(c-1)+1):(2*c))); + if strcmp(modulation,'dqpsk') + tx_symb(c) *= prev_sym_tx(c); + prev_sym_tx(c) = tx_symb(c); + end end - s_ch = s; - + s_ch = tx_symb; + % HF channel simulation ------------------------------------ if hf_sim - % separation between carriers. Note this is - % effectively under samples at Rs, I dont think this - % matters. Equivalent to doing freq shift at Fs, then + % separation between carriers. Note this effectively + % under samples at Rs, I dont think this matters. + % Equivalent to doing freq shift at Fs, then % decimating to Rs. wsep = 2*pi*(1+0.5); % e.g. 75Hz spacing at Rs=50Hz, alpha=0.5 filters - if Nsymb/Nc != floor(Nsymb/Nc) - printf("Error: Nsymb/Nc must be an integrer\n") - return; - end - - % arrange symbols in Nsymb/Nc by Nc matrix - - for i=1:Nc:Nsymb - - % Determine HF channel at each carrier for this symbol - - for k=1:Nc - hf_model(hf_n, k) = hf_gain*(spread(hf_n) + exp(-j*k*wsep*nhfdelay)*spread_2ms(hf_n)); - hf_fading(i+k-1) = abs(hf_model(hf_n, k)); - if hf_mag_only - s_ch(i+k-1) *= abs(hf_model(hf_n, k)); - else - s_ch(i+k-1) *= hf_model(hf_n, k); - end + for c=1:Nc + ahf_model = hf_gain*(spread(hf_n) + exp(-j*c*wsep*nhfdelay)*spread_2ms(hf_n)); + if hf_mag_only + s_ch(c) *= abs(ahf_model); + else + s_ch(c) *= ahf_model; end - hf_n++; + hf_model(hf_n, c) = ahf_model; end + hf_n++; end - + tx_symb_log = [tx_symb_log s_ch]; % AWGN noise and phase/freq offset channel simulation % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im - noise = sqrt(variance*0.5)*(randn(1,Nsymb) + j*randn(1,Nsymb)); - noise_log = [noise_log noise]; - - % organise into carriers to apply frequency and phase offset + noise = sqrt(variance*0.5)*(randn(1,Nc) + j*randn(1,Nc)); + + s_ch = s_ch + noise; - for i=1:Nc:Nsymb - for k=1:Nc - s_ch(i+k-1) = s_ch(i+k-1)*exp(j*phase_offset) + noise(i+k-1); - end - phase_offset += w_offset; - end - % de-modulate rx_bits = zeros(1, framesize); - for i=1:Nsymb - rx_symb = s_ch(i); + for c=1:Nc + + r(c,1:3) = r(c, 2:4); + r(c,4) = s_ch(c); + + rx_symb(c) = s_ch(c); if strcmp(modulation,'dqpsk') - tmp = rx_symb; - rx_symb *= conj(prev_sym_rx/abs(prev_sym_rx)); - prev_sym_rx = tmp; + tmp = rx_symb(c); + rx_symb(c) *= conj(prev_sym_rx(c)/abs(prev_sym_rx(c))); + prev_sym_rx(c) = tmp; end - rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb); - rx_symb_log = [rx_symb_log rx_symb]; - end + % r(c,:) + + if ml == 0 + rx_bits((2*(c-1)+1):(2*c)) = qpsk_demod(rx_symb(c)); + else + tx = [1 j -j -1]; + max_eta = 0; max_symb = tx(1); + for k=1:4 + for k_1=1:4 + for k_2=1:4 + %eta = abs(r(c,1) + r(c,3)*tx(k)'*tx(k_1)' + r(c,2)*tx(k_1)')^2; + eta = abs(r(c,1) + r(c,4)*tx(k)'*tx(k_1)'*tx(k_2)' + r(c,3)*tx(k_1)'*tx(k_2)' + r(c,2)*tx(k_2)')^2; + %printf(" %d %d %f \n", k_1, k, eta); + if eta > max_eta + max_eta = eta; + max_symb = tx(k); + end + end + end + end + rx_bits((2*(c-1)+1):(2*c)) = qpsk_demod(max_symb); + end + + rx_symb_log = [rx_symb_log rx_symb(c)]; + end + % Measure BER error_positions = xor(rx_bits, tx_bits); Nerrs = sum(error_positions); Terrs += Nerrs; Tbits += length(tx_bits); + errors_log = [errors_log error_positions]; + Nerrs_log = [Nerrs_log Nerrs]; end TERvec(ne) = Terrs; BERvec(ne) = Terrs/Tbits; if verbose - printf("EsNo (dB): %f Terrs: %d BER %f BER theory %f", EsNodB, Terrs, - Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2))); + av_tx_pwr = (tx_symb_log * tx_symb_log')/length(tx_symb_log) + + printf("EsNo (dB): %3.1f Terrs: %d BER %4.3f QPSK BER theory %4.3f av_tx_pwr: %3.2f", EsNodB, Terrs, + Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), av_tx_pwr); printf("\n"); end - if verbose > 1 - printf("Terrs: %d BER %f BER theory %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs, - Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), var(tx_symb_log), var(noise_log), - var(tx_symb_log), var(noise_log), var(tx_symb_log)/var(noise_log)); - end end Ebvec = Esvec - 10*log10(bps); sim_out.BERvec = BERvec; sim_out.Ebvec = Ebvec; sim_out.TERvec = TERvec; + sim_out.errors_log = errors_log; if plot_scatter figure(2); @@ -235,34 +218,29 @@ function sim_out = ber_test(sim_in, modulation) plot(real(scat), imag(scat),'+'); title('Scatter plot'); - figure(3); - clf; + if hf_sim + figure(3); + clf; - y = 1:Rs*2; - x = 1:Nc; - EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance); - mesh(x,y,EsNodBSurface); - grid - axis([1 Nc 1 Rs*2 -5 15]) - title('HF Channel Es/No'); - - figure(4); - clf; - %mesh(x,y,unwrap(angle(hf_model(y,:)))); - subplot(211) - plot(y,abs(hf_model(y,1))) - title('HF Channel Carrier 1 Mag'); - subplot(212) - plot(y,angle(hf_model(y,1))) - title('HF Channel Carrier 1 Phase'); - - figure(6) - tmp = []; - for i = 1:hf_n-1 - tmp = [tmp abs(hf_model(i,:))]; + y = 1:(hf_n-1); + x = 1:Nc; + EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance); + EsNodBSurface(find(EsNodBSurface < -5)) = -5; + mesh(x,y,EsNodBSurface); + grid + axis([1 Nc 1 Rs*5 -5 15]) + title('HF Channel Es/No'); + + if verbose + av_hf_pwr = sum(abs(hf_model(y)).^2)/(hf_n-1); + printf("average HF power: %3.2f over %d symbols\n", av_hf_pwr, hf_n-1); + end end - hist(tmp); - end + + figure(4) + clf + stem(Nerrs_log) + end endfunction @@ -290,53 +268,52 @@ function two_bits = qpsk_demod(symbol) two_bits = [bit1 bit0]; endfunction + function sim_in = standard_init sim_in.verbose = 1; sim_in.plot_scatter = 0; sim_in.Esvec = 5; sim_in.Ntrials = 30; - sim_in.framesize = 576; - sim_in.Rs = 100; - sim_in.Nc = 8; - - sim_in.phase_offset = 0; - sim_in.w_offset = 0; - sim_in.phase_noise_amp = 0; + sim_in.Rs = 50; + sim_in.framesize = 8; + sim_in.ml = 0; sim_in.hf_delay_ms = 2; sim_in.hf_sim = 0; - sim_in.hf_phase_only = 0; - sim_in.hf_mag_only = 1; + sim_in.hf_mag_only = 0; endfunction -function test_ideal - sim_in = standard_init(); +function test_curves - sim_in.verbose = 1; - sim_in.plot_scatter = 1; - - sim_in.Esvec = 5; - sim_in.hf_sim = 1; - sim_in.Ntrials = 30; + sim_in = standard_init(); - sim_qpsk_hf = ber_test(sim_in, 'qpsk'); + sim_in.Ntrials = 1000; sim_in.hf_sim = 0; sim_in.plot_scatter = 0; - sim_in.Esvec = 2:5; + sim_in.Esvec = 5:15; Ebvec = sim_in.Esvec - 10*log10(2); BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); - sim_qpsk = ber_test(sim_in, 'qpsk'); + sim_in.ml = 0; sim_dqpsk = ber_test(sim_in, 'dqpsk'); + sim_in.ml = 1; + sim_dqpsk_ml = ber_test(sim_in, 'dqpsk'); + sim_in.hf_sim = 1; + sim_in.ml = 0; + sim_dqpsk_hf = ber_test(sim_in, 'dqpsk'); + sim_in.ml = 1; + sim_dqpsk_ml_hf = ber_test(sim_in, 'dqpsk'); figure(1); clf; semilogy(Ebvec, BER_theory,'r;QPSK theory;') hold on; - semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec,'g;QPSK AWGN;') semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;') + semilogy(sim_dqpsk_ml.Ebvec, sim_dqpsk_ml.BERvec,'k;DQPSK ML AWGN;') + semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'c;DQPSK HF;') + semilogy(sim_dqpsk_ml_hf.Ebvec, sim_dqpsk_ml_hf.BERvec,'k;DQPSK ML HF;') hold off; xlabel('Eb/N0') @@ -346,8 +323,25 @@ function test_ideal endfunction +function test_single + sim_in = standard_init(); + + sim_in.verbose = 1; + sim_in.plot_scatter = 1; + sim_in.Ntrials = 500; + + sim_in.hf_mag_only = 0; + sim_in.hf_sim = 1; + sim_in.ml = 0; + sim_in.Esvec = 10; + + sim_qpsk_hf = ber_test(sim_in, 'dqpsk'); +endfunction + + % Start simulations --------------------------------------- more off; -test_ideal(); +test_curves(); +%test_single(); diff --git a/codec2-dev/octave/test_qpsk2.m b/codec2-dev/octave/test_qpsk2.m index a17e1ee4..bf7bd8e9 100644 --- a/codec2-dev/octave/test_qpsk2.m +++ b/codec2-dev/octave/test_qpsk2.m @@ -515,7 +515,7 @@ function ideal sim_in.hf_sim = 0; sim_in.plot_scatter = 0; - sim_in.Esvec = 2:15; + sim_in.Esvec = 2:10; sim_in.ldpc_code = 0; Ebvec = sim_in.Esvec - 10*log10(2); BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); @@ -523,16 +523,13 @@ function ideal sim_dqpsk = ber_test(sim_in, 'dqpsk'); sim_in.hf_sim = 1; - sim_in.Esvec = 2:15; sim_qpsk_hf = ber_test(sim_in, 'qpsk'); sim_dqpsk_hf = ber_test(sim_in, 'dqpsk'); - sim_in.ldpc_code = 1; - sim_qpsk_hf_ldpc1 = ber_test(sim_in, 'qpsk'); sim_in.ldpc_code_rate = 1/2; - sim_qpsk_hf_ldpc2 = ber_test(sim_in, 'qpsk'); - sim_in.ldpc_code_rate = 3/4; - sim_in.hf_sim = 0; - sim_qpsk_awgn_ldpc = ber_test(sim_in, 'qpsk'); + sim_in.ldpc_code = 1; + sim_qpsk_hf_ldpc = ber_test(sim_in, 'qpsk'); + sim_in.hf_mag_only = 0; + sim_dqpsk_hf_ldpc = ber_test(sim_in, 'dqpsk'); figure(1); clf; @@ -540,11 +537,10 @@ function ideal hold on; semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec,'g;QPSK AWGN;') semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'r;QPSK HF;') - semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;') - semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'m;DQPSK HF;') - semilogy(sim_qpsk_hf_ldpc1.Ebvec, sim_qpsk_hf_ldpc1.BERldpcvec,'k;QPSK HF LDPC 3/4;') - semilogy(sim_qpsk_hf_ldpc2.Ebvec, sim_qpsk_hf_ldpc2.BERldpcvec,'b;QPSK HF LDPC 1/2;') - semilogy(sim_qpsk_awgn_ldpc.Ebvec, sim_qpsk_awgn_ldpc.BERldpcvec,'k;QPSK AWGN LDPC 3/4;') + semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'g;DQPSK AWGN;') + semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'r;DQPSK HF;') + semilogy(sim_qpsk_hf_ldpc.Ebvec, sim_qpsk_hf_ldpc.BERldpcvec,'b;QPSK HF LDPC 1/2;') + semilogy(sim_dqpsk_hf_ldpc.Ebvec, sim_dqpsk_hf_ldpc.BERldpcvec,'b;DQPSK HF LDPC 1/2;') hold off; xlabel('Eb/N0') @@ -636,4 +632,4 @@ endfunction more off; -test_phase_est(); +ideal(); -- 2.25.1