From 3acc72bc0cac694ea6ac6e8dd9d54fcd269c86a6 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sat, 11 Oct 2014 04:20:38 +0000 Subject: [PATCH] combing pilot symbol assisted PSK with DSS, reasonable results git-svn-id: https://svn.code.sf.net/p/freetel/code@1882 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/test_dsss.m | 6 +- codec2-dev/octave/test_dsss_pilot.m | 437 ++++++++++++++++++++++++++++ codec2-dev/octave/test_ml.m | 54 ++-- 3 files changed, 475 insertions(+), 22 deletions(-) create mode 100644 codec2-dev/octave/test_dsss_pilot.m diff --git a/codec2-dev/octave/test_dsss.m b/codec2-dev/octave/test_dsss.m index 3ec08c43..78ed90bc 100644 --- a/codec2-dev/octave/test_dsss.m +++ b/codec2-dev/octave/test_dsss.m @@ -335,10 +335,11 @@ function test_curves sim_in.hf_sim = 1; sim_in.hf_mag_only = 1; sim_qpsk_hf = ber_test(sim_in, 'qpsk'); - sim_in.hf_mag_only = 0; sim_dqpsk_hf = ber_test(sim_in, 'dqpsk'); sim_in.Nchip = 4; sim_dqpsk_hf_dsss = ber_test(sim_in, 'dqpsk'); + sim_in.hf_mag_only = 1; + sim_qpsk_hf_dsss = ber_test(sim_in, 'qpsk'); figure(1); clf; @@ -348,6 +349,7 @@ function test_curves semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'b;QPSK HF;') semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'k;DQPSK HF;') semilogy(sim_dqpsk_hf_dsss.Ebvec, sim_dqpsk_hf_dsss.BERvec,'g;DQPSK DSSS HF;') + semilogy(sim_qpsk_hf_dsss.Ebvec, sim_qpsk_hf_dsss.BERvec,'r;QPSK DSSS HF;') hold off; xlabel('Eb/N0') @@ -404,5 +406,5 @@ endfunction more off; -test_1600_v_450(); +%test_1600_v_450(); test_curves(); diff --git a/codec2-dev/octave/test_dsss_pilot.m b/codec2-dev/octave/test_dsss_pilot.m new file mode 100644 index 00000000..1d3ec169 --- /dev/null +++ b/codec2-dev/octave/test_dsss_pilot.m @@ -0,0 +1,437 @@ +% test_dsss_pilot.m +% David Rowe Oct 2014 +% + +% Simulation to test FDM QPSK with pilot based coherent detection +% combined with DSSS. + +% reqd to make sure we can repeat tests exactly + +rand('state',1); +randn('state',1); + +1; + +% main test function + +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_mag_only = sim_in.hf_mag_only; + Nchip = sim_in.Nchip; + Np = sim_in.Np; + Ns = sim_in.Ns; + + bps = 2; + Nc = Nsymb = framesize/bps; + + prev_sym_tx = qpsk_mod([0 0])*ones(1,Nc*Nchip); + prev_sym_rx = qpsk_mod([0 0])*ones(1,Nc*Nchip); + + tx_bits_mem = zeros(Np*Ns+1, framesize); + tx_symb_mem = zeros(Np*Ns+1, Nc*Nchip); + s_ch_mem = zeros(Np*Ns+1, Nc*Nchip); + + % Init HF channel model from stored sample files of spreading signal ---------------------------------- + + % convert "spreading" samples from 1kHz carrier at Fs to complex + % baseband, generated by passing a 1kHz sine wave through PathSim + % with the ccir-poor model, enabling one path at a time. + + Fc = 1000; M = Fs/Rs; + fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb"); + spread1k = fread(fspread, "int16")/10000; + fclose(fspread); + fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb"); + spread1k_2ms = fread(fspread, "int16")/10000; + fclose(fspread); + + % down convert to complex baseband + spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))'); + spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))'); + + % remove -2000 Hz image + 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 + + spread = spread(1000:length(spread)); + spread_2ms = spread_2ms(1000:length(spread_2ms)); + + % decimate down to Rs + + spread = spread(1:M:length(spread)); + spread_2ms = spread_2ms(1:M:length(spread_2ms)); + + % Determine "gain" of HF channel model, so we can normalise + % carrier power during HF channel sim to calibrate SNR. I imagine + % different implementations of ccir-poor would do this in + % different ways, leading to different BER results. Oh Well! + + hf_gain = 1.0/sqrt(var(spread(1:Ntrials))+var(spread_2ms(1:Ntrials))); + + % Start Simulation ---------------------------------------------------------------- + + for ne = 1:length(Esvec) + EsNodB = Esvec(ne); + EsNo = 10^(EsNodB/10); + + variance = 1/EsNo; + if verbose > 1 + printf("EsNo (dB): %f EsNo: %f variance: %f\n", EsNodB, EsNo, variance); + end + + Terrs = 0; Tbits = 0; + + s_ch_tx_log = []; + rx_symb_log = []; + noise_log = []; + errors_log = []; + Nerrs_log = []; + + % init HF channel + + hf_n = 1; + phi_ = zeros(Ntrials+Np*Ns, Nc*Nchip); + + phase_offset = 0; + w_offset = 0; + + % simulation starts here----------------------------------- + + for nn = 1:Ntrials+Np*Ns + + tx_bits = round( rand( 1, framesize) ); + tx_bits_mem(1:Np*Ns,:) = tx_bits_mem(2:Np*Ns+1,:); + for b=1:framesize + tx_bits_mem(Np*Ns+1,b) = tx_bits(b); + end + + % modulate -------------------------------------------- + + tx_symb = zeros(1,Nc*Nchip); + + for i=1:Nc + tx_symb(i) = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); + end + + % Optionally copy to other carriers (spreading) + + for i=Nc+1:Nc:Nc*Nchip + tx_symb(i:i+Nc-1) = tx_symb(1:Nc); + end + + % Optionally DQPSK encode + + if strcmp(modulation,'dqpsk') + for c=1:Nc*Nchip + tx_symb(c) *= prev_sym_tx(c); + prev_sym_tx(c) = tx_symb(c); + end + end + + tx_symb_mem(1:Np*Ns,:) = tx_symb_mem(2:Np*Ns+1,:); + for c=1:Nc*Nchip + tx_symb_mem(Np*Ns+1,c) = tx_symb(c); + end + + s_ch = tx_symb/sqrt(Nchip); + + % HF channel simulation ------------------------------------ + + if hf_sim + + % 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 + + hf_model(hf_n, :) = zeros(1,Nc*Nchip); + + for i=1:Nchip + time_shift = floor(i*Rs/4); + for k=1:Nc + ahf_model = hf_gain*(spread(hf_n+time_shift) + exp(-j*k*wsep*nhfdelay)*spread_2ms(hf_n+time_shift)); + if hf_mag_only + s_ch((i-1)*Nc+k) *= abs(ahf_model); + else + s_ch((i-1)*Nc+k) *= ahf_model; + end + hf_model(hf_n, (i-1)*Nc+k) = ahf_model; + end + end + hf_n++; + end + + s_ch_tx_log = [s_ch_tx_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*Nchip) + j*randn(1,Nsymb*Nchip)); + noise_log = [noise_log noise]; + + s_ch = s_ch.*exp(j*phase_offset) + noise; + phase_offset += w_offset; + + s_ch_mem(1:Np*Ns,:) = s_ch_mem(2:Np*Ns+1,:); + for c=1:Nc*Nchip + s_ch_mem(Np*Ns+1,c) = s_ch(c); + end + + % pilot based phase estimation and correction + + if Np + for c=1:Nc*Nchip + pilots = tx_symb_mem(:,c); + pilots(floor(Np*Ns/2)+1) = 0; + phi_(nn,c) = angle(pilots(1:Ns:Np*Ns+1)'*s_ch_mem(1:Ns:Np*Ns+1,c)); + end + end + + % demodulate stage 1 + + for c=1:Nc*Nchip + rx_symb(c) = s_ch_mem(floor(Np*Ns/2)+1,c) *exp(-j*phi_(nn,c)); + if strcmp(modulation,'dqpsk') + tmp = rx_symb(c); + rx_symb(c) *= conj(prev_sym_rx(c)/abs(prev_sym_rx(c))); + prev_sym_rx(c) = tmp; + end + end + + % de-spread + + for i=Nc+1:Nc:Nchip*Nc + rx_symb(1:Nc) = rx_symb(1:Nc) + rx_symb(i:i+Nc-1); + end + + % demodulate stage 2 (start when we have Np*Ns symbols in memories) + + if nn > Np*Ns + rx_bits = zeros(1, framesize); + for c=1:Nc + rx_bits((2*(c-1)+1):(2*c)) = qpsk_demod(rx_symb(c)); + tx_bits((2*(c-1)+1):(2*c)) = tx_bits_mem(floor(Np*Ns/2)+1,(2*(c-1)+1):(2*c)); + end + rx_symb_log = [rx_symb_log rx_symb(1:Nc)]; + + % 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 + end + + TERvec(ne) = Terrs; + BERvec(ne) = Terrs/Tbits; + if verbose + av_tx_pwr = (s_ch_tx_log * s_ch_tx_log')/length(s_ch_tx_log); + + printf("EsNo (dB): %3.1f Terrs: %d BER %4.2f QPSK BER theory %4.2f 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); + clf; + scat = rx_symb_log .* exp(j*pi/4); + plot(real(scat), imag(scat),'+'); + title('Scatter plot'); + + if hf_sim + figure(3); + clf; + + y = 1:(hf_n-1); + x = 1:Nc*Nchip; + EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance); + EsNodBSurface(find(EsNodBSurface < -5)) = -5; + mesh(x,y,EsNodBSurface); + grid + axis([1 (Nc+1)*Nchip 1 Rs*5 -5 15]) + title('HF Channel Es/No'); + + if verbose + [m n] = size(hf_model); + av_hf_pwr = sum(sum(abs(hf_model(:,:)).^2))/(m*n); + printf("average HF power: %3.2f over %d symbols\n", av_hf_pwr, m*n); + end + + figure(5); + clf + subplot(211) + [m n] = size(hf_model); + plot(angle(hf_model(1:m,1)),'g;HF channel phase;') + hold on; + lphi_ = length(phi_); + plot(phi_(1+floor(Ns*Np/2):lphi_),'r+;Estimated HF channel phase;') + ylabel('Phase (rads)'); + subplot(212) + plot(abs(hf_model(1:m,1))) + ylabel('Amplitude'); + xlabel('Time (symbols)'); + end + + figure(4) + clf + stem(Nerrs_log) + end + +endfunction + +% Gray coded QPSK modulation function + +function symbol = qpsk_mod(two_bits) + two_bits_decimal = sum(two_bits .* [2 1]); + switch(two_bits_decimal) + case (0) symbol = 1; + case (1) symbol = j; + case (2) symbol = -j; + case (3) symbol = -1; + endswitch +endfunction + +% Gray coded QPSK demodulation function + +function two_bits = qpsk_demod(symbol) + if isscalar(symbol) == 0 + printf("only works with scalars\n"); + return; + end + bit0 = real(symbol*exp(j*pi/4)) < 0; + bit1 = imag(symbol*exp(j*pi/4)) < 0; + 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 = 2; + sim_in.Rs = 50; + + sim_in.phase_offset = 0; + sim_in.w_offset = 0; + sim_in.phase_noise_amp = 0; + + sim_in.hf_delay_ms = 2; + sim_in.hf_sim = 0; + sim_in.hf_mag_only = 0; + + sim_in.Nchip = 1; +endfunction + +function test_curves + + sim_in = standard_init(); + + sim_in.verbose = 1; + sim_in.plot_scatter = 1; + + sim_in.Esvec = 10; + sim_in.hf_sim = 1; + sim_in.Ntrials = 1000; + sim_in.Rs = 200; + sim_in.Np = 4; + sim_in.Ns = 8; + sim_in.Nchip = 1; + + sim_qpsk = ber_test(sim_in, 'qpsk'); + + sim_in.hf_sim = 0; + sim_in.plot_scatter = 0; + sim_in.Esvec = 10:20; + Ebvec = sim_in.Esvec - 10*log10(2); + BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); + + sim_in.Np = 0; + sim_in.Nchip = 1; + + sim_dqpsk = ber_test(sim_in, 'dqpsk'); + sim_in.hf_sim = 1; + sim_in.hf_mag_only = 1; + sim_qpsk_hf_ideal = ber_test(sim_in, 'qpsk'); + sim_in.hf_mag_only = 0; + sim_dqpsk_hf = ber_test(sim_in, 'dqpsk'); + sim_in.Np = 6; + sim_qpsk_hf_pilot = ber_test(sim_in, 'qpsk'); + sim_in.Nchip = 2; + sim_qpsk_hf_pilot_dsss = ber_test(sim_in, 'qpsk'); + + figure(1); + clf; + semilogy(Ebvec, BER_theory,'r;QPSK theory;') + hold on; + semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;') + semilogy(sim_qpsk_hf_ideal.Ebvec, sim_qpsk_hf_ideal.BERvec,'b;QPSK HF ideal;') + semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'k;DQPSK HF;') + semilogy(sim_qpsk_hf_pilot.Ebvec, sim_qpsk_hf_pilot.BERvec,'r;QPSK Np=6 HF;') + semilogy(sim_qpsk_hf_pilot_dsss.Ebvec, sim_qpsk_hf_pilot_dsss.BERvec,'g;QPSK Np=6 Nchip=2 HF;') + hold off; + + xlabel('Eb/N0') + ylabel('BER') + grid("minor") + axis([min(Ebvec) max(Ebvec) 1E-3 1]) +endfunction + +function test_single + + sim_in = standard_init(); + + sim_in.verbose = 1; + sim_in.plot_scatter = 1; + + sim_in.Ntrials = 1000; + sim_in.Esvec = 2; + sim_in.hf_sim = 0; + sim_in.hf_mag_only = 0; + sim_in.Nchip = 2; + sim_in.Np = 6; + sim_in.Ns = 8; + sim_in.Rs = 200; + + sim_qpsk_hf = ber_test(sim_in, 'qpsk'); +endfunction + + +% Start simulations --------------------------------------- + +more off; +test_curves(); +%test_single(); diff --git a/codec2-dev/octave/test_ml.m b/codec2-dev/octave/test_ml.m index 51dd8cb6..a735264b 100644 --- a/codec2-dev/octave/test_ml.m +++ b/codec2-dev/octave/test_ml.m @@ -2,10 +2,15 @@ % David Rowe Oct 2014 % -% Simulation to test FDM QPSK with maximum likelihood decoding on -% fading channels. +% Simulation to test FDM DQPSK with maximum likelihood detection on +% fading channels. Based on JPL publication 89-38 "Multiple Symbol +% Differential Detection of Uncoded and Trellis Coded MPSK" by +% Divsalar, Simon, Shahshahani. Thanks Johhn Gibbs NN7F for advice. -1; +% reqd to make sure we can repeat tests exactly + +rand('state',1); +randn('state',1); % main test function @@ -147,9 +152,9 @@ function sim_out = ber_test(sim_in, modulation) rx_bits = zeros(1, framesize); for c=1:Nc - r(c,1:3) = r(c, 2:4); - r(c,4) = s_ch(c); - + r(c,2:4) = r(c, 1:3); + r(c,1) = s_ch(c); + rx_symb(c) = s_ch(c); if strcmp(modulation,'dqpsk') tmp = rx_symb(c); @@ -157,8 +162,6 @@ function sim_out = ber_test(sim_in, modulation) prev_sym_rx(c) = tmp; end - % r(c,:) - if ml == 0 rx_bits((2*(c-1)+1):(2*c)) = qpsk_demod(rx_symb(c)); else @@ -167,8 +170,12 @@ function sim_out = ber_test(sim_in, modulation) 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; + if ml == 3 + eta = abs(r(c,3) + r(c,1)*tx(k)'*tx(k_1)' + r(c,2)*tx(k_1)')^2; + end + if ml == 4 + eta = abs(r(c,4) + r(c,1)*tx(k)'*tx(k_1)'*tx(k_2)' + r(c,2)*tx(k_1)'*tx(k_2)' + r(c,3)*tx(k_2)')^2; + end %printf(" %d %d %f \n", k_1, k, eta); if eta > max_eta max_eta = eta; @@ -197,7 +204,7 @@ function sim_out = ber_test(sim_in, modulation) BERvec(ne) = Terrs/Tbits; if verbose - av_tx_pwr = (tx_symb_log * tx_symb_log')/length(tx_symb_log) + 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); @@ -297,23 +304,30 @@ function test_curves Ebvec = sim_in.Esvec - 10*log10(2); BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); sim_in.ml = 0; - sim_dqpsk = ber_test(sim_in, 'dqpsk'); - sim_in.ml = 1; - sim_dqpsk_ml = ber_test(sim_in, 'dqpsk'); + %sim_dqpsk = ber_test(sim_in, 'dqpsk'); + sim_in.ml = 3; + %sim_dqpsk_ml3 = ber_test(sim_in, 'dqpsk'); + sim_in.ml = 4; + %sim_dqpsk_ml4 = ber_test(sim_in, 'dqpsk'); sim_in.hf_sim = 1; + sim_in.hf_mag_only = 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'); + sim_in.ml = 4; + sim_dqpsk_ml_hf3 = ber_test(sim_in, 'dqpsk'); + sim_in.ml = 4; + sim_dqpsk_ml_hf4 = ber_test(sim_in, 'dqpsk'); figure(1); clf; semilogy(Ebvec, BER_theory,'r;QPSK theory;') hold on; - 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.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;') + %semilogy(sim_dqpsk_ml3.Ebvec, sim_dqpsk_ml3.BERvec,'k;DQPSK ML N=3 AWGN;') + %semilogy(sim_dqpsk_ml4.Ebvec, sim_dqpsk_ml4.BERvec,'g;DQPSK ML N=4 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;') + semilogy(sim_dqpsk_ml_hf3.Ebvec, sim_dqpsk_ml_hf3.BERvec,'k;DQPSK ML N=3 HF;') + semilogy(sim_dqpsk_ml_hf4.Ebvec, sim_dqpsk_ml_hf4.BERvec,'g;DQPSK ML N=4 HF;') hold off; xlabel('Eb/N0') @@ -331,7 +345,7 @@ function test_single sim_in.Ntrials = 500; sim_in.hf_mag_only = 0; - sim_in.hf_sim = 1; + sim_in.hf_sim = 0; sim_in.ml = 0; sim_in.Esvec = 10; -- 2.25.1