From bf8c6911ca6f7e77b49e8e6cb67ee1da2961b15b Mon Sep 17 00:00:00 2001 From: drowe67 Date: Fri, 7 Apr 2017 04:47:30 +0000 Subject: [PATCH] progressing rate Fs bpsk ofdm simulation git-svn-id: https://svn.code.sf.net/p/freetel/code@3093 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/bpsk_hf_fs.m | 293 +++++++++++++++++++++++++++++++++ 1 file changed, 293 insertions(+) create mode 100644 codec2-dev/octave/bpsk_hf_fs.m diff --git a/codec2-dev/octave/bpsk_hf_fs.m b/codec2-dev/octave/bpsk_hf_fs.m new file mode 100644 index 00000000..942c016b --- /dev/null +++ b/codec2-dev/octave/bpsk_hf_fs.m @@ -0,0 +1,293 @@ +% bpsk_hf_fs.m +% David Rowe Mar 2017 +% +% Rate Fs BPSK simulation, development of bpsk_hf_rs.m + +#{ + TODO: + [X] strip back experimental stuff to just features we need + [X] ZOH/integrator + [X] OFDM up and down conversion + [ ] rate Fs HF model and HF results + [ ] handle boarder carriers + [ ] start with phantom carriers + + but unhappy with 1800Hz bandwidth + [ ] also try interpolation or just single row + [ ] compute SNR and PAPR + [ ] acquisition & freq offset estimation + [ ] SSB bandpass filtering +#} + +1; + +function sim_out = run_sim(sim_in) + Rs = 100; + Fs = 8000; + M = Fs/Rs; + + EbNodB = sim_in.EbNodB; + verbose = sim_in.verbose; + hf_en = sim_in.hf_en; + + Ns = sim_in.Ns; % step size for pilots + Nc = sim_in.Nc; % Number of cols, aka number of carriers + + Nbitsperframe = (Ns-1)*Nc; + printf("Nbitsperframe: %d\n", Nbitsperframe); + Nrowsperframe = Nbitsperframe/Nc; + printf("Nrowsperframe: %d\n", Nrowsperframe); + + % Important to define run time in seconds so HF model will evolve the same way + % for different pilot insertion rates. So lets work backwards from approx + % seconds in run to get Nbits, the total number of payload data bits + + % frame has Ns-1 data symbols between pilots, e.g. for Ns=3: + % + % PPP + % DDD + % DDD + % PPP + + Nrows = sim_in.Nsec*Rs; + Nframes = floor((Nrows-1)/Ns); + Nbits = Nframes * Nbitsperframe; % number of payload data bits + + Nr = Nbits/Nc; % Number of data rows to get Nbits total + + if verbose + printf("Nc.....: %d\n", Nc); + printf("Ns.....: %d (step size for pilots, Ns-1 data symbols between pilots)\n", Ns); + printf("Nr.....: %d\n", Nr); + printf("Nbits..: %d\n", Nbits); + end + + % double check if Nbits fit neatly into carriers + + assert(Nbits/Nc == floor(Nbits/Nc), "Nbits/Nc 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); + + % simulate for each Eb/No point ------------------------------------ + + for nn=1:length(EbNodB) + rand('seed',1); + randn('seed',1); + + EbNo = 10 .^ (EbNodB(nn)/10); + variance = 1/(M*EbNo/2); + + % 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); + + % map to BPSK symbols + + tx_sym = 2*tx_bits - 1; + +#{ + % upsample to rate Fs using Zero Order Hold (ZOH) + + tx_sym_os = zeros(Nrp*M, Nc+2); + for r=1:Nrp + for rr=(r-1)*M+1:r*M + tx_sym_os(rr,:) = tx_sym(r,:); + end + end +#} + + % OFDM up conversion and upsampling to rate Fs + + w = (0:Nc+1)*2*pi*Rs/Fs; + W = zeros(M,Nc+2); + for c=1:Nc+2 + W(:,c) = exp(j*w(c)*(0:M-1)); + end + + Nsam = Nrp*M; + tx = zeros(Nrp*M,1); + for r=1:Nrp + for c=1:Nc+2 + acarrier = tx_sym(r,c) * W(:,c); + tx((r-1)*M+1:r*M) += acarrier/M; + end + end + + rx = tx; + + noise = sqrt(variance)*(0.5*randn(Nrp*M,1) + j*0.5*randn(Nrp*M,1)); + rx += noise; + + % downconvert, downsample and integrate using OFDM + + rx_sym = zeros(Nrp, Nc+2); + for r=1:Nrp + for c=1:Nc+2 + acarrier = rx((r-1)*M+1:r*M) .* conj(W(:,c)); + rx_sym(r,c) = sum(acarrier); + end + end + + % pilot based phase est, we use known tx symbols as pilots ---------- + + phase_est_pilot_log = 10*ones(Nrp,Nc+2); + phase_est_stripped_log = 10*ones(Nrp,Nc+2); + 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_pilot_rect = sum(rx_sym(r,cr)*tx_sym(r,cr)') + sum(rx_sym(r+Ns,cr)*tx_sym(r+Ns,cr)'); + + % use next step of pilots in past and future + + if r > Ns+1 + aphase_est_pilot_rect += sum(rx_sym(r-Ns,cr)*tx_sym(r-Ns,cr)'); + end + if r < Nrp - 2*Ns + aphase_est_pilot_rect += sum(rx_sym(r+2*Ns,cr)*tx_sym(r+2*Ns,cr)'); + end + + % correct phase offset using phase estimate + + for rr=r+1:r+Ns-1 + aphase_est_pilot = angle(aphase_est_pilot_rect); + phase_est_pilot_log(rr,c) = aphase_est_pilot; + rx_corr(rr,c) = rx_sym(rr,c) * exp(-j*aphase_est_pilot); + end + + end % r=1:Ns:Nrp-Ns + end % c=2:Nc+1 + + % remove pilots to give us just data symbols + + rx_np = []; + for r=1:Nrp + if mod(r-1,Ns) != 0 + rx_np = [rx_np; rx_corr(r,2:Nc+1)]; + end + end + + % calculate BER stats as a block, after pilots extracted + + 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, Nbits, Nerrs); + + if verbose + figure(1) + plot(real(tx)) + figure(2) + Tx = abs(fft(tx.*hanning(Nsam))); + Tx_dB = 20*log10(Tx); + dF = Fs/Nsam; + plot((1:Nsam)*dF, Tx_dB); + mx = max(Tx_dB); + axis([0 Fs/2 mx-60 mx]) + + figure(3); clf; + plot(rx_np,'+'); + axis([-2 2 -2 2]); + + if hf_en + figure(4); clf; + plot(abs(hf_model(:,2:Nc+1))); + end + + figure(5); clf; + plot(phase_est_log(:,2:Nc+1),'+', 'markersize', 10); + hold on; + plot(phase_est_pilot_log(:,2:Nc+1),'g+', 'markersize', 5); + if sim_in.hf_en + plot(angle(hf_model(:,2:Nc+1))); + end + axis([1 Nrp -pi pi]); + end + + sim_out.ber(nn) = sum(Nerrs)/Nbits; + sim_out.pilot_overhead = 10*log10(Ns/(Ns-1)); + end +endfunction + + +% Plot BER against Eb/No curves at various pilot insertion rates Ns +% using the HF multipath channel. Second set of curves includes Eb/No +% loss for pilot insertion, so small Ns means better tracking of phase +% but large pilot insertion loss + +% Target operating point Eb/No is 6dB, as this is where our rate 1/2 +% LDPC code gives good results (10% PER, 1% BER). However this means +% the Eb/No at the input is 10*log(1/2) or 3dB less, so we need to +% make sure phase est works at Eb/No = 6 - 3 = 3dB + +function run_curves_hf + sim_in.Nc = 8; + sim_in.Ns = 8; + sim_in.Nsec = 240; + sim_in.EbNodB = 1:8; + sim_in.verbose = 0; + sim_in.hf_en = 0; + + EbNoLin = 10.^(sim_in.EbNodB/10); + hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); + + hf_sim = run_sim(sim_in); + + figure(4); clf; + semilogy(sim_in.EbNodB, hf_theory,'b+-;HF theory;'); + hold on; + semilogy(sim_in.EbNodB, hf.ber,'c+-;Ns=8 HFsim;'); + hold off; + axis([1 8 4E-2 2E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); +end + + +function run_single + sim_in.Nc = 8; + sim_in.Ns = 8; + sim_in.Nsec = 10; + sim_in.EbNodB = 4; + sim_in.verbose = 1; + sim_in.hf_en = 0; + + run_sim(sim_in); +end + + +format; +more off; + + +run_single +%run_curves_hf + + + + -- 2.25.1