From 49f5842cf4cb686861fa8f2dbf47e3a768feab77 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Thu, 30 Mar 2017 04:40:35 +0000 Subject: [PATCH] reasonable results with HF channel model, about 1dB IL on HF at 6dB Eb/No. Next step try interpolating or pilots from > 2 frames git-svn-id: https://svn.code.sf.net/p/freetel/code@3085 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/bpsk_hf.m | 222 ++++++++++++++++++++++++------------ 1 file changed, 147 insertions(+), 75 deletions(-) diff --git a/codec2-dev/octave/bpsk_hf.m b/codec2-dev/octave/bpsk_hf.m index 7bf0f5ed..fa5ffcd6 100644 --- a/codec2-dev/octave/bpsk_hf.m +++ b/codec2-dev/octave/bpsk_hf.m @@ -6,8 +6,16 @@ #{ TODO: - [ ] sim pilot based phase est using known symbols - [ ] test AWGN BER with averaging pilots from adj carriers + [X] sim pilot based phase est using known symbols + [X] test AWGN BER with averaging pilots from adj carriers + [X] refactor to insert pilot rows + [X] add border cols, not used for data + [X] centre est on current carrier, extend to > 3 + [ ] test single points + [ ] first pass curves taking into account pilot losses + [ ] consider interpolation + [ ] consider combining mod stripping phase est techniques inside frame + [ ] remove border carriers, interpolate edge carrier #} 1; @@ -19,13 +27,44 @@ function sim_out = run_sim(sim_in) EbNodB = sim_in.EbNodB; verbose = sim_in.verbose; hf_en = sim_in.hf_en; + hf_phase = sim_in.hf_phase; Ns = sim_in.Ns; % step size for pilots Nc = sim_in.Nc; % Number of cols, aka number of carriers Nr = Nbits/Nc; % Number of rows to get Nbits total phase_offset = sim_in.phase_offset; + if verbose + printf("Nbits: %d\n", Nbits); + printf("Nc: %d\n", Nc); + printf("Nr: %d\n", Nr); + printf("Ns: %d (step size for pilots, Ns-1 data symbols between pilots)\n", Ns); + end + + % check if Nbits fit neatly into carriers + assert(Nbits/Nc == floor(Nbits/Nc), "Nbits/Nc must be an integer"); + + % check if bits fit neatly into frames with rows of pilots above and below + % PPP + % DDD + % DDD + % PPP + Nbitsperframe = (Ns-1)*Nc; + printf("Nbitsperframe: %d\n", Nbitsperframe); + Nframes = Nbits/Nbitsperframe; + Nrowsperframe = Nbitsperframe/Nc; + printf("Nrowsperframe: %d\n", Nrowsperframe); + + % check if Nbits fit neatly into frames delineated by pilots + + assert(Nframes == floor(Nframes), "Nbits/Nbits/frame must be an integer"); + printf("Nframes: %d\n", Nframes); + + Nrp = Nr + Nframes + 1; % number of rows once pilots inserted + % extra row of pilots at end + printf("Nrp: %d (number of rows including pilots)\n", Nrp); + % set up HF model if hf_en @@ -34,13 +73,31 @@ function sim_out = run_sim(sim_in) dopplerSpreadHz = 1.0; path_delay = 1E-3*Rs; - spread1 = doppler_spread(dopplerSpreadHz, Rs, Nr*2); - spread2 = doppler_spread(dopplerSpreadHz, Rs, Nr*2); + spread1 = doppler_spread(dopplerSpreadHz, Rs, Nrp+10); + spread2 = doppler_spread(dopplerSpreadHz, Rs, Nrp+10); + + % sometimes doppler_spread() doesn't return exactly the number of samples we need + + assert(length(spread1) >= Nrp, "not enough doppler spreading samples"); + assert(length(spread2) >= Nrp, "not enough doppler spreading samples"); + hf_gain = 1.0/sqrt(var(spread1)+var(spread2)); % printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1)); end - % simulate for each Eb/No point + % construct an artificial phase countour for testing, linear across freq and time + + if sim_in.phase_test + phase_test = ones(Nrp, Nc+2); + for r=1:Nrp + for c=1:Nc+2 + phase_test(r,c) = -pi/2 + c*pi/(Nc+2) + r*0.01*2*pi; + phase_test(r,c) = phase_test(r,c) - 2*pi*floor((phase_test(r,c)+pi)/(2*pi)); + end + end + end + + % simulate for each Eb/No point ------------------------------------ for nn=1:length(EbNodB) rand('seed',1); @@ -48,62 +105,56 @@ function sim_out = run_sim(sim_in) EbNo = 10 .^ (EbNodB(nn)/10); variance = 1/(EbNo/2); - noise = sqrt(variance)*(0.5*randn(Nr,Nc) + j*0.5*randn(Nr,Nc)); + noise = sqrt(variance)*(0.5*randn(Nrp,Nc+2) + j*0.5*randn(Nrp,Nc+2)); - tx_bits = rand(Nr,Nc) > 0.5; - tx = 1 - 2*tx_bits; + % generate tx bits and insert pilot rows and border cols + + tx_bits = []; tx_bits_np = []; + for f=1:Nframes + tx_bits = [tx_bits; ones(1,Nc+2)]; + for r=1:Nrowsperframe + arowofbits = rand(1,Nc) > 0.5; + tx_bits = [tx_bits; [1 arowofbits 1]]; + tx_bits_np = [tx_bits_np; arowofbits]; + end + end + tx_bits = [tx_bits; [1 ones(1,Nc) 1]]; % final row of pilots + [nr nc] = size(tx_bits); + assert(nr == Nrp); + %tx_bits + + tx = 2*tx_bits - 1; rx = tx * exp(j*phase_offset); + if sim_in.phase_test + rx = rx .* exp(j*phase_test); + end + 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. + % ISI, just freq filtering. % Note Rs carrier spacing, sample rate is Rs - hf_model = zeros(Nr,Nc); phase_est = zeros(Nr,Nc); - for i=1:Nr - for c=1:Nc + hf_model = zeros(Nr,Nc+2); phase_est = zeros(Nr,Nc); + for r=1:Nrp + for c=1:Nc+2 w = 2*pi*c*Rs/Rs; - hf_model(i,c) = hf_gain*(spread1(i) + exp(-j*w*path_delay)*spread2(i)); + hf_model(r,c) = hf_gain*(spread1(r) + exp(-j*w*path_delay)*spread2(r)); end - - if sim_in.av_phase - - % Simulate obtaining phase est by averaging hf channel model - % samples at carrier frequencies adjacent to each carrier. - % In case of edge carriers we sample HF channel model to one - % side of carrier for now, for a real world estimate we'd - % need to interpolate phase to edge carriers somehow - - for c=1:Nc - phase_sum = 0; - for cc = c-1:c+1 - w = 2*pi*cc*Rs/Rs; - ahf_model = hf_gain*(spread1(i) + exp(-j*w*path_delay)*spread2(i)); - phase_sum += ahf_model; - end - phase_est(i,c) = angle(phase_sum); - end - rx(i,:) = rx(i,:) .* hf_model(i,:) .* exp(-j*phase_est(i,:)); - + + if hf_phase + rx(r,:) = rx(r,:) .* hf_model(r,:); else - - % just apply amplitude fading, assume we have ideal phase est - - rx(i,:) = rx(i,:) .* abs(hf_model(i,:)); + rx(r,:) = rx(r,:) .* abs(hf_model(r,:)); end end end rx += noise; - Nbits_np = Nbits; - rx_np = rx; - tx_bits_np = tx_bits; - % pilot based phase est, we use known tx symbols as pilots ---------- if sim_in.pilot_phase_est @@ -112,14 +163,19 @@ function sim_out = run_sim(sim_in) % adjust phase of data symbol % demodulate and count errors of just data - phase_est_log = zeros(Nr,Nc); - for c=1:Nc - for r=1:Ns:Nr-Ns - - % estimate phase using average of two pilots at either end of frame - - %aphase_est_rect = rx(r,c)*tx(r,c)' + rx(r+Ns,c)*tx(r+Ns,c)'; - aphase_est_rect = sum(rx(r,:)*tx(r,:)') + sum(rx(r+Ns,:)*tx(r+Ns,:)'); + phase_est_log = 10*ones(Nrp,Nc+2); + for c=2:Nc+1 + for r=1:Ns:Nrp-Ns + + % estimate phase using average of 6 pilots in a rect 2D window centred + % on this carrier + % PPP + % DDD + % DDD + % PPP + + cr = c-1:c+1; + aphase_est_rect = sum(rx(r,cr)*tx(r,cr)') + sum(rx(r+Ns,cr)*tx(r+Ns,cr)'); aphase_est = angle(aphase_est_rect); % apply this phase estimate to symbols in frame, and check for error @@ -130,41 +186,57 @@ function sim_out = run_sim(sim_in) end end end + end - % build up tx_bits and rx symbols matrices without pilots + % strip pilots to give us just data symbols - tx_bits_np = []; rx_np = []; Nbits_np = 0; - for r=1:Nr - if mod(r-1,Ns) != 0 - tx_bits_np = [tx_bits_np; tx_bits(r,:)]; - rx_np = [rx_np; rx(r,:)]; - Nbits_np += Nc; - end + rx_np = []; + for r=1:Nrp + if mod(r-1,Ns) != 0 + rx_np = [rx_np; rx(r,2:Nc+1)]; end end + %phase_test + %phase_est_log + % calculate BER stats as a block, after pilots extracted - rx_bits_np = real(rx_np) < 0; + rx_bits_np = real(rx_np) > 0; errors = xor(tx_bits_np, rx_bits_np); Nerrs = sum(sum(errors)); - printf("EbNodB: %3.2f BER: %4.3f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/Nbits_np, Nbits_np, Nerrs); + printf("EbNodB: %3.2f BER: %4.3f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/Nbits, Nbits, Nerrs); if verbose figure(1); clf; plot(rx_np,'+'); axis([-2 2 -2 2]); - if sim_in.pilot_phase_est - figure(2); clf; - plot(phase_est_log,'+'); - end if hf_en figure(2); clf; - plot(abs(hf_model)); - figure(3); clf; - plot(angle(hf_model)); - hold on; plot(phase_est,'+'); hold off; + plot(abs(hf_model(:,2:Nc+1))); + if sim_in.hf_phase + figure(3); clf; + plot(angle(hf_model(:,2:Nc+1))); + if sim_in.pilot_phase_est + hold on; plot(phase_est_log(:,2:Nc+1),'+', 'markersize', 10); hold off; + end + axis([1 Nrp -pi pi]); + end + else + if sim_in.pilot_phase_est + figure(3); clf; + plot(phase_est_log,'+'); + end + end + + if sim_in.phase_test + figure(3); clf; + plot(phase_test(:,2:Nc+1)); + if sim_in.pilot_phase_est + hold on; plot(phase_est_log(:,2:Nc+1),'+', 'markersize', 10); hold off; + end + axis([1 Nrp -pi pi]); end end @@ -199,16 +271,16 @@ end function run_single - data_rows = 1000; sim_in.Nc = 3; - sim_in.Ns = 6; - sim_in.Nbits = (data_rows*sim_in.Ns+1)*sim_in.Nc; - sim_in.EbNodB = 2; + sim_in.Ns = 8; + sim_in.Nbits = 5000*sim_in.Nc*(sim_in.Ns-1); + sim_in.EbNodB = 6; sim_in.verbose = 1; sim_in.pilot_phase_est = 1; - sim_in.phase_offset = -pi; - sim_in.hf_en = 0; - sim_in.av_phase = 0; + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 1; + sim_in.hf_phase = 1; run_sim(sim_in); end -- 2.25.1