From: drowe67 Date: Sun, 23 Apr 2017 08:04:48 +0000 (+0000) Subject: starting to refactor rate Fs ofdm simulation, and changed name of rate fs and rate rs X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=d3af0c792bc9e665531551795eb6a4772be5918d;p=freetel-svn-tracking.git starting to refactor rate Fs ofdm simulation, and changed name of rate fs and rate rs git-svn-id: https://svn.code.sf.net/p/freetel/code@3101 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/bpsk_hf_fs.m b/codec2-dev/octave/bpsk_hf_fs.m deleted file mode 100644 index 4ace45fb..00000000 --- a/codec2-dev/octave/bpsk_hf_fs.m +++ /dev/null @@ -1,697 +0,0 @@ -% 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 - [X] rate Fs HF model and HF results - [X] add QPSK - [X] add CP - [ ] adjust waveform parameters for real world - [ ] Nsec run time take into account CP - [ ] plot phase ests - [ ] handle border carriers - [ ] start with phantom carriers - + but unhappy with 1800Hz bandwidth - [ ] also try interpolation or just single row - [ ] compute SNR and PAPR - [ ] timing estimator - [ ] acquisition & freq offset estimation - [ ] SSB bandpass filtering -#} - -1; - -% 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) - bit0 = real(symbol*exp(j*pi/4)) < 0; - bit1 = imag(symbol*exp(j*pi/4)) < 0; - two_bits = [bit1 bit0]; -endfunction - - -% Correlates the OFDM pilot symbol samples with a window of received -% samples to determine the most likely timing offset. Combines two -% frames pilots so we need at least Nsamperframe+M+Ncp samples in rx. -% Also determines frequency offset at maximimum correlation. Can be -% used for acquisition (coarse timing a freq offset), and fine timing - -function [t_est foff_est] = coarse_sync(states, rx, rate_fs_pilot_samples) - Nsamperframe = states.Nsamperframe; Fs = states.Fs; - Npsam = length(rate_fs_pilot_samples); - verbose = states.verbose; - - Ncorr = length(rx) - (Nsamperframe+Npsam) + 1; - assert(Ncorr > 0); - corr = zeros(1,Ncorr); - for i=1:Ncorr - corr(i) = abs(rx(i:i+Npsam-1) * rate_fs_pilot_samples'); - corr(i) += abs(rx(i+Nsamperframe:i+Nsamperframe+Npsam-1) * rate_fs_pilot_samples'); - end - - [mx t_est] = max(abs(corr)); - - C = abs(fft(rx(t_est:t_est+Npsam-1) .* conj(rate_fs_pilot_samples), Fs)); - C += abs(fft(rx(t_est+Nsamperframe:t_est+Nsamperframe+Npsam-1) .* conj(rate_fs_pilot_samples), Fs)); - - fmax = 30; - [mx_pos foff_est_pos] = max(C(1:fmax)); - [mx_neg foff_est_neg] = max(C(Fs-fmax+1:Fs)); - - if mx_pos > mx_neg - foff_est = foff_est_pos - 1; - else - foff_est = foff_est_neg - fmax - 1; - end - - if verbose - %printf("t_est: %d\n", t_est); - figure(7); clf; - plot(abs(corr)) - figure(8) - plot(C) - axis([0 200 0 0.4]) - end -endfunction - - -function [sim_out rate_fs_pilot_samples rx] = run_sim(sim_in) - Fs = 8000; - Rs = sim_in.Rs; - Tcp = sim_in.Tcp; - M = Fs/Rs; - Ncp = Tcp*Fs; - woffset = 2*pi*sim_in.foff_hz/Fs; - - bps = sim_in.bps; - EbNodB = sim_in.EbNodB; - verbose = sim_in.verbose; - hf_en = sim_in.hf_en; - timing_en = sim_in.timing_en; - - Ns = sim_in.Ns; % step size for pilots - Nc = sim_in.Nc; % Number of cols, aka number of carriers - - Nbitsperframe = (Ns-1)*Nc*bps; - Nrowsperframe = Nbitsperframe/(Nc*bps); - if verbose - printf("Rs:.........: %4.2f\n", Rs); - printf("M:..........: %d\n", M); - printf("bps:.........: %d\n", bps); - printf("Nbitsperframe: %d\n", Nbitsperframe); - printf("Nrowsperframe: %d\n", Nrowsperframe); - end - - % 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*bps); % Number of data rows to get Nbits total - - % double check if Nbits fit neatly into carriers - - assert(Nbits/(Nc*bps) == floor(Nbits/(Nc*bps)), "Nbits/(Nc*bps) must be an integer"); - - Nrp = Nr + Nframes + 1; % number of rows once pilots inserted - % extra row of pilots at end - Nsamperframe = (Nrowsperframe+1)*(M+Ncp); - - 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); - printf("Nframes: %d\n", Nframes); - printf("Nrp....: %d (number of rows including pilots)\n", Nrp); - end - - % generate same pilots each time - - rand('seed',1); - pilots = 1 - 2*(rand(1,Nc+2) > 0.5); - - % set up HF model --------------------------------------------------------------- - - if hf_en - - % some typical values, or replace with user supplied - - dopplerSpreadHz = 1.0; path_delay_ms = 1; - - if isfield(sim_in, "dopplerSpreadHz") - dopplerSpreadHz = sim_in.dopplerSpreadHz; - end - if isfield(sim_in, "path_delay_ms") - path_delay_ms = sim_in.path_delay_ms; - end - path_delay_samples = path_delay_ms*Fs/1000; - printf("Doppler Spread: %3.2f Hz Path Delay: %3.2f ms %d samples\n", dopplerSpreadHz, path_delay_ms, path_delay_samples); - - % generate same fading pattern for every run - - randn('seed',1); - - spread1 = doppler_spread(dopplerSpreadHz, Fs, (sim_in.Nsec*(M+Ncp)/M+0.2)*Fs); - spread2 = doppler_spread(dopplerSpreadHz, Fs, (sim_in.Nsec*(M+Ncp)/M+0.2)*Fs); - - % sometimes doppler_spread() doesn't return exactly the number of samples we need - - assert(length(spread1) >= Nrp*M, "not enough doppler spreading samples"); - assert(length(spread2) >= Nrp*M, "not enough doppler spreading samples"); - - hf_gain = 1.0/sqrt(var(spread1)+var(spread2)); - % printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1)); - end - - % init timing est states - - tstates.Nsamperframe = Nsamperframe; tstates.M = M; tstates.Ncp = Ncp; - tstates.verbose = 0; tstates.Fs = Fs; - delta_t = []; - window_width = 11; % search window_width/2 from current timing instant - - % simulate for each Eb/No point ------------------------------------ - - for nn=1:length(EbNodB) - rand('seed',1); - randn('seed',1); - - EbNo = bps * (10 .^ (EbNodB(nn)/10)); - variance = 1/(M*EbNo/2); - - % generate tx bits - - tx_bits = rand(1,Nbits) > 0.5; - - % map to symbols in linear array - - if bps == 1 - tx_sym_lin = 2*tx_bits - 1; - end - if bps == 2 - for s=1:Nbits/bps - tx_sym_lin(s) = qpsk_mod(tx_bits(2*(s-1)+1:2*s)); - end - end - - % place symbols in multi-carrier frame with pilots and boundary carriers - - tx_sym = []; s = 1; - for f=1:Nframes - aframe = zeros(Nrowsperframe,Nc+2); - aframe(1,:) = pilots; - for r=1:Nrowsperframe - arowofsymbols = tx_sym_lin(s:s+Nc-1); - s += Nc; - aframe(r+1,2:Nc+1) = arowofsymbols; - end - tx_sym = [tx_sym; aframe]; - end - tx_sym = [tx_sym; pilots]; % final row of pilots - [nr nc] = size(tx_sym); - assert(nr == Nrp); - - % OFDM up conversion and upsampling to rate Fs --------------------- - - w = (0:Nc+1)*2*pi*Rs/Fs; - W = zeros(Nc+2,M); - for c=1:Nc+2 - W(c,:) = exp(j*w(c)*(0:M-1)); - end - - Nsam = Nrp*(M+Ncp); - tx = []; - - % OFDM upconvert symbol by symbol so we can add CP - - for r=1:Nrp - asymbol = tx_sym(r,:) * W/M; - asymbol_cp = [asymbol(M-Ncp+1:M) asymbol]; - tx = [tx asymbol_cp]; - end - assert(length(tx) == Nsam); - - % OFDM symbol of pilots is used for coarse timing and freq during acquisition, and fine timing - - rate_fs_pilot_samples = tx(1:M+Ncp); - - % channel simulation --------------------------------------------------------------- - - tx = resample(tx, 2000, 2000); - tx = [tx zeros(1,Nsam-length(tx))]; - tx = tx(1:Nsam); - rx = tx; - - if hf_en - %rx = [zeros(1,path_delay_samples) tx(1:Nsam-path_delay_samples)]; - rx = hf_gain * tx(1:Nsam) .* spread1(1:Nsam); - rx += hf_gain * [zeros(1,path_delay_samples) tx(1:Nsam-path_delay_samples)] .* spread2(1:Nsam); - end - - rx = rx .* exp(j*woffset*(1:Nsam)); - - noise = sqrt(variance)*(0.5*randn(1,Nsam) + j*0.5*randn(1,Nsam)); - rx += noise; - - % some spare samples at end to allow for timing est window - - rx = [rx zeros(1,Nsamperframe)]; - - % pilot based phase est, we use known tx symbols as pilots ---------- - - rx_sym = zeros(Nrp, Nc+2); - 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); - timing_est_log = []; - timing_est = 1; - if timing_en - sample_point = timing_est; - else - sample_point = Ncp/2; - end - - for r=1:Ns:Nrp-Ns - - if timing_en - - % update timing every frame - - if (mod(r-1,Ns) == 0) && (r != 1) && (r != Nrp) - st = (r-1)*(M+Ncp) + 1 - floor(window_width/2) + (timing_est-1); - en = st + Nsamperframe-1 + M+Ncp + window_width-1; - ft_est = coarse_sync(tstates, rx(st:en), rate_fs_pilot_samples); - timing_est += ft_est - ceil(window_width/2); - if verbose - printf("ft_est: %2d timing_est: %2d sample_point: %2d\n", ft_est, timing_est, sample_point); - end - delta_t = [delta_t ft_est - ceil(window_width/2)]; - sample_point = max(timing_est+Ncp/4, sample_point); - sample_point = min(timing_est+Ncp, sample_point); - end - end - - timing_est_log = [timing_est_log timing_est]; - - % down convert at current timing instant - - % previous pilot - - if r > Ns+1 - rr = r-Ns; - st = (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1; - for c=1:Nc+2 - acarrier = rx(st:en) .* conj(W(c,:)); - rx_sym(rr,c) = sum(acarrier); - end - end - - % pilot - this frame - pilot - - for rr=r:r+Ns - - st = (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1; - for c=1:Nc+2 - acarrier = rx(st:en) .* conj(W(c,:)); - rx_sym(rr,c) = sum(acarrier); - end - end - - % next pilot - - if r < Nrp - 2*Ns - rr = r+2*Ns; - st = (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1; - for c=1:Nc+2 - acarrier = rx(st:en) .* conj(W(c,:)); - rx_sym(rr,c) = sum(acarrier); - end - end - - % OK - now estimate and correct phase - - for c=2:Nc+1 - - % 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 % c=2:Nc+1 - end % r=1:Ns:Nrp-Ns - - - % remove pilots to give us just data symbols and demodulate - - rx_bits = []; rx_np = []; - for r=1:Nrp - if mod(r-1,Ns) != 0 - arowofsymbols = rx_corr(r,2:Nc+1); - rx_np = [rx_np arowofsymbols]; - if bps == 1 - arowofbits = real(arowofsymbols) > 0; - end - if bps == 2 - arowofbits = zeros(1,Nc); - for c=1:Nc - arowofbits((c-1)*2+1:c*2) = qpsk_demod(arowofsymbols(c)); - end - end - rx_bits = [rx_bits arowofbits]; - end - end - assert(length(rx_bits) == Nbits); - - - % calculate BER stats as a block, after pilots extracted - - errors = xor(tx_bits, rx_bits); - Nerrs = 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(1:Nsam).*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]); - title('Scatter'); - - if hf_en - figure(4); clf; - subplot(211) - plot(abs(spread1(1:Nsam))); - %hold on; plot(abs(spread2(1:Nsam)),'g'); hold off; - subplot(212) - plot(angle(spread1(1:Nsam))); - title('spread1 amp and phase'); - 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); - title('Phase est'); - -#{ - % todo, work out a way to plot rate Fs hf model phase - if sim_in.hf_en - plot(angle(hf_model(:,2:Nc+1))); - end -#} - - axis([1 Nrp -pi pi]); - - figure(6); clf; - subplot(211) - stem(delta_t) - subplot(212) - plot(timing_est_log); - title('Timing'); - end - - sim_out.ber(nn) = sum(Nerrs)/Nbits; - sim_out.pilot_overhead = 10*log10(Ns/(Ns-1)); - sim_out.M = M; sim_out.Fs = Fs; sim_out.Ncp = Ncp; - sim_out.Nrowsperframe = Nrowsperframe; sim_out.Nsamperframe = Nsamperframe; - end -endfunction - - -function run_single - Ts = 0.016; - sim_in.Tcp = 0.002; - sim_in.Rs = 1/Ts; sim_in.bps = 2; sim_in.Nc = 16; sim_in.Ns = 8; - - %sim_in.Nsec = 3*(sim_in.Ns+1)/sim_in.Rs; % one frame - sim_in.Nsec = 30; - - sim_in.EbNodB = 3; - sim_in.verbose = 1; - sim_in.hf_en = 1; - sim_in.foff_hz = 0; - sim_in.timing_en = 1; - - run_sim(sim_in); -end - - -% Plot BER against Eb/No curves for AWGN and HF - -% Target operating point Eb/No for HF 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 -% -% For AWGN target is 2dB so -1dB op point. - -function run_curves - Ts = 0.010; - sim_in.Rs = 1/Ts; - sim_in.Tcp = 0.002; - sim_in.bps = 2; sim_in.Ns = 8; sim_in.Nc = 8; sim_in.verbose = 0; - sim_in.foff_hz = 0; - - pilot_overhead = (sim_in.Ns-1)/sim_in.Ns; - cp_overhead = Ts/(Ts+sim_in.Tcp); - overhead_dB = -10*log10(pilot_overhead*cp_overhead); - - sim_in.hf_en = 0; - sim_in.Nsec = 30; - sim_in.EbNodB = -3:5; - awgn_EbNodB = sim_in.EbNodB; - - awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNodB/10))); - awgn = run_sim(sim_in); - - sim_in.hf_en = 1; - sim_in.Nsec = 120; - sim_in.EbNodB = 1:8; - - EbNoLin = 10.^(sim_in.EbNodB/10); - hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); - - hf = run_sim(sim_in); - - figure(4); clf; - semilogy(awgn_EbNodB, awgn_theory,'b+-;AWGN theory;'); - hold on; - semilogy(sim_in.EbNodB, hf_theory,'b+-;HF theory;'); - semilogy(awgn_EbNodB+overhead_dB, awgn_theory,'g+-;AWGN lower bound with pilot + CP overhead;'); - semilogy(sim_in.EbNodB+overhead_dB, hf_theory,'g+-;HF lower bound with pilot + CP overhead;'); - semilogy(awgn_EbNodB+overhead_dB, awgn.ber,'r+-;AWGN sim;'); - semilogy(sim_in.EbNodB+overhead_dB, hf.ber,'r+-;HF sim;'); - hold off; - axis([-3 8 1E-2 2E-1]) - xlabel('Eb/No (dB)'); - ylabel('BER'); - grid; grid minor on; - legend('boxoff'); -end - - -% Run an acquisition test, returning vectors of estimation errors - -function [delta_t delta_foff] = acquisition_test(Ntests=10, EbNodB=100, foff_hz=0, hf_en=0, fine_en=0) - - % generate test signal at a given Eb/No and frequency offset - - Ts = 0.016; - sim_in.Tcp = 0.002; - sim_in.Rs = 1/Ts; sim_in.bps = 2; sim_in.Nc = 16; sim_in.Ns = 8; - - sim_in.Nsec = Ntests*(sim_in.Ns+1)/sim_in.Rs; - - sim_in.EbNodB = EbNodB; - sim_in.verbose = 0; - sim_in.hf_en = hf_en; - sim_in.foff_hz = foff_hz; sim_in.timing_en = 0; - - % set up acquistion - - Nsamperframe = states.Nsamperframe = sim_out.Nsamperframe; - states.M = sim_out.M; states.Ncp = sim_out.Ncp; - states.verbose = 0; - states.Fs = sim_out.Fs; - - % test fine or acquisition over test signal - #{ - fine: - start with coarse timing instant - - on each frame est timing a few samples about that point - - update timing instant - - corr: - where is best plcase to sample - - just before end of symbol? - - how long should sequence be? - - add extra? - - aim for last possible moment? - - man I hope IL isn't too big..... - #} - - delta_t = []; delta_t = []; delta_foff = []; - - if fine_en - - window_width = 5; % search +/-2 samples from current timing instant - timing_instant = Nsamperframe+1; % start at correct instant for AWGN - % start at second frame so we can search -2 ... +2 - - while timing_instant < (length(rx) - (Nsamperframe + length(rate_fs_pilot_samples) + window_width)) - st = timing_instant - ceil(window_width/2); - en = st + Nsamperframe-1 + length(rate_fs_pilot_samples) + window_width-1; - [ft_est foff_est] = coarse_sync(states, rx(st:en), rate_fs_pilot_samples); - printf("ft_est: %d timing_instant %d %d\n", ft_est, timing_instant, mod(timing_instant, Nsamperframe)); - timing_instant += Nsamperframe + ft_est - ceil(window_width/2); - delta_t = [delta_ft ft_est - ceil(window_width/2)]; - end - else - % for coarse simulation we just use contant window shifts - - st = 0.5*Nsamperframe; - en = 2.5*Nsamperframe - 1; - ct_target = Nsamperframe/2; - - for w=1:Nsamperframe:length(rx)-3*Nsamperframe - %st = w+0.5*Nsamperframe; en = st+2*Nsamperframe-1; - %[ct_est foff_est] = coarse_sync(states, rx(st:en), rate_fs_pilot_samples); - [ct_est foff_est] = coarse_sync(states, rx(w+st:w+en), rate_fs_pilot_samples); - printf("ct_est: %4d foff_est: %3.1f\n", ct_est, foff_est); - - % valid coarse timing ests are modulo Nsamperframe - - delta_t = [delta_ct ct_est-ct_target]; - delta_foff = [delta_foff (foff_est-foff_hz)]; - end - end - -endfunction - - -% Run some tests for various acquisition conditions. Probability of -% acquistion is what matters, e.g. if it's 50% we can expect sync -% within 2 frames -% P(t)/P(f) P(t)/P(f) -% Eb/No AWGN HF -% +/- 25Hz -1/3 1.0/0.3 0.96/0.3 -% +/- 5Hz -1/3 1.0/0.347 0.96/0.55 -% +/- 25Hz 10/10 1.00/0.92 0.99/0.77 - -function acquisition_histograms(fine_en = 0) - Fs = 8000; - Ntests = 100; - - % allowable tolerance for acquistion - - ftol_hz = 2.0; - ttol_samples = 0.002*Fs; - - % AWGN channel operating point - - [dct dfoff] = acquisition_test(Ntests, -1, 25, 0, fine_en); - - % Probability of acquistion is what matters, e.g. if it's 50% we can - % expect sync within 2 frames - - printf("AWGN P(time offset acq) = %3.2f\n", length(find (abs(dct) < ttol_samples))/length(dct)); - if fine_en == 0 - printf("AWGN P(freq offset acq) = %3.2f\n", length(find (abs(dfoff) < ftol_hz))/length(dfoff)); - end - - figure(1) - hist(dct(find (abs(dct) < ttol_samples))) - if fine_en == 0 - figure(2) - hist(dfoff) - end - - % HF channel operating point - - [dct dfoff] = acquisition_test(Ntests, 3, 25, 1, fine_en); - - printf("HF P(time offset acq) = %3.2f\n", length(find (abs(dct) < ttol_samples))/length(dct)); - if fine_en == 0 - printf("HF P(freq offset acq) = %3.2f\n", length(find (abs(dfoff) < ftol_hz))/length(dfoff)); - end - - figure(3) - hist(dct(find (abs(dct) < ttol_samples))) - if fine_en == 0 - figure(4) - hist(dfoff) - end - -endfunction - - - -% choose simulation to run here ------------------------------------------------------- - -format; -more off; - -run_single -%run_curves -%acquisition_histograms(1) - diff --git a/codec2-dev/octave/bpsk_hf_rs.m b/codec2-dev/octave/bpsk_hf_rs.m deleted file mode 100644 index 8b013fe8..00000000 --- a/codec2-dev/octave/bpsk_hf_rs.m +++ /dev/null @@ -1,786 +0,0 @@ -% bpsk_hf_rs.m -% David Rowe Mar 2017 -% -% Rate Rs BPSK simulation to explore techniques for phase estimation -% over multiple carriers in HF channel - -#{ - TODO: - [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 - [X] test single points - + 1dB IL @ 6dB HF, 0.4 dB @ 2dB AWGN - [X] try linear interpolation - [X] try longer time windows - [X] try combining mod stripping phase est inside frame - [X] curves taking into account pilot losses - [ ] remove border carriers, interpolate edge carrier - [ ] modify for QPSK - [ ] change name -#} - -1; - -% 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) - bit0 = real(symbol*exp(j*pi/4)) < 0; - bit1 = imag(symbol*exp(j*pi/4)) < 0; - two_bits = [bit1 bit0]; -endfunction - - -function sim_out = run_sim(sim_in) - Rs = 100; - - bps = sim_in.bps; - EbNodB = sim_in.EbNodB; - verbose = sim_in.verbose; - hf_en = sim_in.hf_en; - hf_phase = sim_in.hf_phase; - phase_offset = sim_in.phase_offset; - - Ns = sim_in.Ns; % step size for pilots - Nc = sim_in.Nc; % Number of cols, aka number of carriers - - Nbitsperframe = (Ns-1)*Nc*bps; - Nrowsperframe = Nbitsperframe/(Nc*bps); - if verbose - printf("bps:.........: %d\n", bps); - printf("Nbitsperframe: %d\n", Nbitsperframe); - printf("Nrowsperframe: %d\n", Nrowsperframe); - end - - % 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*bps); % 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*bps) == floor(Nbits/(Nc*bps)), "Nbits/(Nc*bps) 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 - - % some typical values, or replace with user supplied - - dopplerSpreadHz = 1.0; path_delay = 1E-3*Rs; - - if isfield(sim_in, "dopplerSpreadHz") - dopplerSpreadHz = sim_in.dopplerSpreadHz; - end - if isfield(sim_in, "path_delay") - path_delay = sim_in.path_delay; - end - printf("Doppler Spread: %3.2f Hz Path Delay: %3.2f symbols\n", dopplerSpreadHz, path_delay); - randn('seed',1); - spread1 = doppler_spread(dopplerSpreadHz, Rs, sim_in.Nsec*Rs*1.1); - spread2 = doppler_spread(dopplerSpreadHz, Rs, sim_in.Nsec*Rs*1.1); - - % 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 - - % 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); - randn('seed',1); - - EsNo = bps * (10 .^ (EbNodB(nn)/10)); - variance = 1/(EsNo/2); - noise = sqrt(variance)*(0.5*randn(Nrp,Nc+2) + j*0.5*randn(Nrp,Nc+2)); - - % generate tx bits - - tx_bits = rand(1,Nbits) > 0.5; - - % map to symbols - - if bps == 1 - tx_symb = 2*tx_bits - 1; - end - if bps == 2 - for s=1:Nbits/bps - tx_symb(s) = qpsk_mod(tx_bits(2*(s-1)+1:2*s)); - end - end - - % place symbols in multi-carrier frame with pilots and boundary carriers - - tx = []; s = 1; - for f=1:Nframes - aframe = zeros(Nrowsperframe,Nc+2); - aframe(1,:) = 1; - for r=1:Nrowsperframe - arowofsymbols = tx_symb(s:s+Nc-1); - s += Nc; - aframe(r+1,2:Nc+1) = arowofsymbols; - end - tx = [tx; aframe]; - end - tx = [tx; ones(1,Nc+2)]; % final row of pilots - [nr nc] = size(tx); - assert(nr == Nrp); - - 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. - - % Note Rs carrier spacing, sample rate is Rs - - 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(r,c) = hf_gain*(spread1(r) + exp(-j*w*path_delay)*spread2(r)); - end - - if hf_phase - rx(r,:) = rx(r,:) .* hf_model(r,:); - else - rx(r,:) = rx(r,:) .* abs(hf_model(r,:)); - end - end - end - - rx += noise; - - % pilot based phase est, we use known tx symbols as pilots ---------- - - rx_corr = rx; - - if sim_in.pilot_phase_est - - % est phase from pilots either side of data symbols - % adjust phase of data symbol - % demodulate and count errors of just data - - 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_rect1 = sum(rx(r,cr)*tx(r,cr)'); - aphase_est_pilot_rect2 = sum(rx(r+Ns,cr)*tx(r+Ns,cr)'); - - % optionally use next step of pilots in past and future - - if sim_in.pilot_wide - if r > Ns+1 - aphase_est_pilot_rect1 += sum(rx(r-Ns,cr)*tx(r-Ns,cr)'); - end - if r < Nrp - 2*Ns - aphase_est_pilot_rect2 += sum(rx(r+2*Ns,cr)*tx(r+2*Ns,cr)'); - end - end - - % correct phase offset using phase estimate - - for rr=r+1:r+Ns-1 - a = b = 1; - if sim_in.pilot_interp - b = (rr-r)/Ns; a = 1 - b; - end - %printf("rr: %d a: %4.3f b: %4.3f\n", rr, a, b); - aphase_est_pilot = angle(a*aphase_est_pilot_rect1 + b*aphase_est_pilot_rect2); - phase_est_pilot_log(rr,c) = aphase_est_pilot; - rx_corr(rr,c) = rx(rr,c) * exp(-j*aphase_est_pilot); - end - - if sim_in.stripped_phase_est - % Optional modulation stripping feed fwd phase estimation, to refine - % pilot-based phase estimate. Doing it after pilot based phase estimation - % means we don't need to deal with ambiguity, which is difficult to handle - % in low SNR channels. - - % Use vector of 7 symbols around current data symbol. We could use a 2D - % window if we can work out how best to correct with pilot-est and avoid - % ambiguities - - for rr=r+1:r+Ns-1 - - % extract a matrix of nearby samples with pilot-based offset removed - - amatrix = rx(max(1,rr-3):min(Nrp,rr+3),c) .* exp(-j*aphase_est_pilot); - - % modulation strip and est phase - - stripped = abs(amatrix) .* exp(j*2*angle(amatrix)); - aphase_est_stripped = angle(sum(sum(stripped)))/2; - phase_est_stripped_log(rr,c) = aphase_est_stripped; - - % correct rx symbols based on both phase ests - - phase_est_log(rr,c) = angle(exp(j*(aphase_est_pilot+aphase_est_stripped))); - rx_corr(rr,c) = rx(rr,c) * exp(-j*phase_est_log(rr,c)); - end - end % sim_in.stripped_phase_est - - end % r=1:Ns:Nrp-Ns - - end % c=2:Nc+1 - end % sim_in.pilot_phase_est - - - if isfield(sim_in, "ml_pd") && sim_in.ml_pd - - % Bill's ML with pilots phase detector, does phase est and demodulation - - rx_bits = []; rx_np = []; - aframeofbits = zeros(Ns-1, Nc); - for r=1:Ns:Nrp-Ns - - % demodulate this frame, ML operates carrier by carrier - - for c=2:Nc+1 - arxcol = rx(r:r+Ns, c); - arxcol(1) = rx(r, c-1) + rx(r, c+1); - arxcol(Ns+1) = rx(r+Ns, c-1) + rx(r+Ns, c+1); - [acolofbits aphase_est] = ml_pd(rot90(arxcol), bps, [1 Ns+1]); - aframeofbits(:,c-1) = xor(acolofbits, ones(1,Ns-1)); - rx_np = [rx_np rot90(arxcol) .* exp(-j*aphase_est)]; - end - - % unpack from frame into linear array of bits - - for rr=1:Ns-1 - rx_bits = [rx_bits aframeofbits(rr,:)]; - end - - end - else - - % remove pilots to give us just data symbols and demodulate - - rx_bits = []; rx_np = []; - for r=1:Nrp - if mod(r-1,Ns) != 0 - arowofsymbols = rx_corr(r,2:Nc+1); - rx_np = [rx_np arowofsymbols]; - if bps == 1 - arowofbits = real(arowofsymbols) > 0; - end - if bps == 2 - arowofbits = zeros(1,Nc); - for c=1:Nc - arowofbits((c-1)*2+1:c*2) = qpsk_demod(arowofsymbols(c)); - end - end - rx_bits = [rx_bits arowofbits]; - end - end - end - %tx_bits - %rx_bits - assert(length(rx_bits) == Nbits); - - %phase_test - %phase_est_log - - % calculate BER stats as a block, after pilots extracted - - errors = xor(tx_bits, rx_bits); - Nerrs = sum(errors); - - 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 hf_en - figure(2); clf; - plot(abs(hf_model(:,2:Nc+1))); - end - - if sim_in.pilot_phase_est - figure(3); 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.stripped_phase_est - plot(phase_est_stripped_log(:,2:Nc+1),'ro', 'markersize', 5); - end - if sim_in.hf_en && sim_in.hf_phase - plot(angle(hf_model(:,2:Nc+1))); - end - if sim_in.phase_test - plot(phase_test(:,2:Nc+1)); - end - axis([1 Nrp -pi pi]); - end - 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 = 7; - sim_in.Ns = 5; - sim_in.Nsec = 240; - sim_in.EbNodB = 1:8; - sim_in.verbose = 0; - sim_in.pilot_phase_est = 0; - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 0; - sim_in.stripped_phase_est = 0; - sim_in.phase_offset = 0; - sim_in.phase_test = 0; - sim_in.hf_en = 1; - sim_in.hf_phase = 0; - - sim_in.Ns = 5; - hf_ref_Ns_5_no_phase = run_sim(sim_in); - sim_in.Ns = 9; - hf_ref_Ns_9_no_phase = run_sim(sim_in); - - sim_in.hf_phase = 1; - sim_in.pilot_phase_est = 1; - - sim_in.Ns = 5; - hf_Ns_5 = run_sim(sim_in); - - sim_in.Ns = 9; - hf_Ns_9 = run_sim(sim_in); - - sim_in.Ns = 17; - hf_Ns_17 = run_sim(sim_in); - - figure(4); clf; - semilogy(sim_in.EbNodB, hf_ref_Ns_5_no_phase.ber,'b+-;Ns=5 HF ref no phase;'); - hold on; - semilogy(sim_in.EbNodB, hf_ref_Ns_9_no_phase.ber,'c+-;Ns=9 HF ref no phase;'); - semilogy(sim_in.EbNodB, hf_Ns_5.ber,'g+--;Ns=5;'); - semilogy(sim_in.EbNodB + hf_Ns_5.pilot_overhead, hf_Ns_5.ber,'go-;Ns=5 with pilot overhead;'); - semilogy(sim_in.EbNodB, hf_Ns_9.ber,'r+--;Ns=9;'); - semilogy(sim_in.EbNodB + hf_Ns_9.pilot_overhead, hf_Ns_9.ber,'ro-;Ns=9 with pilot overhead;'); - semilogy(sim_in.EbNodB, hf_Ns_17.ber,'k+--;Ns=17;'); - semilogy(sim_in.EbNodB + hf_Ns_17.pilot_overhead, hf_Ns_17.ber,'ko-;Ns=17 with pilot overhead;'); - hold off; - axis([1 8 4E-2 2E-1]) - xlabel('Eb/No (dB)'); - ylabel('BER'); - grid; grid minor on; - legend('boxoff'); - title('HF Multipath 1Hz Doppler 1ms delay'); - -end - - -% Generate HF curves for some alternative, experimental methods tested -% during development, such as interpolation, refinements using -% modulation stripping, narrow window. - -function run_curves_hf_alt - sim_in.Nc = 7; - sim_in.Ns = 5; - sim_in.Nsec = 60; - sim_in.EbNodB = 1:8; - sim_in.verbose = 0; - sim_in.pilot_phase_est = 0; - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 0; - sim_in.stripped_phase_est = 0; - sim_in.phase_offset = 0; - sim_in.phase_test = 0; - sim_in.hf_en = 1; - sim_in.hf_phase = 0; - - sim_in.Ns = 9; - hf_ref_Ns_9_no_phase = run_sim(sim_in); - - sim_in.hf_phase = 1; - sim_in.pilot_phase_est = 1; - hf_Ns_9 = run_sim(sim_in); - - sim_in.stripped_phase_est = 1; - hf_Ns_9_stripped = run_sim(sim_in); - - sim_in.stripped_phase_est = 0; - sim_in.pilot_wide = 0; - hf_Ns_9_narrow = run_sim(sim_in); - - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 1; - hf_Ns_9_interp = run_sim(sim_in); - - figure(6); clf; - semilogy(sim_in.EbNodB, hf_ref_Ns_9_no_phase.ber,'c+-;Ns=9 HF ref no phase;'); - hold on; - semilogy(sim_in.EbNodB, hf_Ns_9.ber,'r+--;Ns=9;'); - semilogy(sim_in.EbNodB, hf_Ns_9_stripped.ber,'g+--;Ns=9 stripped refinement;'); - semilogy(sim_in.EbNodB, hf_Ns_9_narrow.ber,'b+--;Ns=9 narrow;'); - semilogy(sim_in.EbNodB, hf_Ns_9_interp.ber,'k+--;Ns=9 interp;'); - hold off; - axis([1 8 4E-2 2E-1]) - xlabel('Eb/No (dB)'); - ylabel('BER'); - grid; grid minor on; - legend('boxoff'); - title('HF Multipath 1Hz Doppler 1ms delay'); - -end - - -% Generate HF curves for fixed Ns but different HF channels. - -function run_curves_hf_channels - sim_in.Nc = 7; - sim_in.Ns = 9; - sim_in.Nsec = 240; - sim_in.EbNodB = 1:8; - sim_in.verbose = 0; - sim_in.pilot_phase_est = 0; - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 0; - sim_in.stripped_phase_est = 0; - sim_in.phase_offset = 0; - sim_in.phase_test = 0; - sim_in.hf_en = 1; - sim_in.hf_phase = 0; - - hf_Ns_9_1hz_1ms_no_phase = run_sim(sim_in); - - sim_in.hf_phase = 1; - sim_in.pilot_phase_est = 1; - hf_Ns_9_1hz_1ms = run_sim(sim_in); - - Rs = 100; - - sim_in.dopplerSpreadHz = 1.0; - sim_in.path_delay = 500E-6*Rs; - hf_Ns_9_1hz_500us = run_sim(sim_in); - - sim_in.dopplerSpreadHz = 1.0; - sim_in.path_delay = 2E-3*Rs; - hf_Ns_9_1hz_2ms = run_sim(sim_in); - - sim_in.dopplerSpreadHz = 2.0; - sim_in.path_delay = 1E-3*Rs; - hf_Ns_9_2hz_1ms = run_sim(sim_in); - - sim_in.dopplerSpreadHz = 2.0; - sim_in.path_delay = 1E-3*Rs; - hf_Ns_9_2hz_2ms = run_sim(sim_in); - - sim_in.dopplerSpreadHz = 2.0; - sim_in.path_delay = 2E-3*Rs; - hf_Ns_9_2hz_2ms = run_sim(sim_in); - - sim_in.dopplerSpreadHz = 4.0; - sim_in.path_delay = 1E-3*Rs; - hf_Ns_9_4hz_1ms = run_sim(sim_in); - - figure(6); clf; - semilogy(sim_in.EbNodB, hf_Ns_9_1hz_1ms_no_phase.ber,'c+-;Ns=9 1Hz 1ms ref no phase;'); - hold on; - semilogy(sim_in.EbNodB, hf_Ns_9_1hz_500us.ber,'k+-;Ns=9 1Hz 500us;'); - semilogy(sim_in.EbNodB, hf_Ns_9_1hz_1ms.ber,'r+-;Ns=9 1Hz 1ms;'); - semilogy(sim_in.EbNodB, hf_Ns_9_1hz_2ms.ber,'bo-;Ns=9 1Hz 2ms;'); - semilogy(sim_in.EbNodB, hf_Ns_9_2hz_1ms.ber,'g+-;Ns=9 2Hz 1ms;'); - semilogy(sim_in.EbNodB, hf_Ns_9_2hz_2ms.ber,'mo-;Ns=9 2Hz 2ms;'); - semilogy(sim_in.EbNodB, hf_Ns_9_4hz_1ms.ber,'c+-;Ns=9 4Hz 1ms;'); - hold off; - axis([1 8 4E-2 2E-1]) - xlabel('Eb/No (dB)'); - ylabel('BER'); - grid; grid minor on; - legend('boxoff'); - title('HF Multipath Ns = 9'); - -end - - -% AWGN curves for BPSK and QPSK. Coded Eb/No operating point is 2dB, -% so raw BER for rate 1/2 will be -1dB - -function run_curves_awgn_bpsk_qpsk - sim_in.Nc = 7; - sim_in.Ns = 7; - sim_in.Nsec = 30; - sim_in.verbose = 0; - sim_in.pilot_phase_est = 0; - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 0; - sim_in.stripped_phase_est = 1; - sim_in.phase_offset = 0; - sim_in.phase_test = 0; - sim_in.hf_en = 0; - sim_in.hf_phase = 0; - - sim_in.EbNodB = -3:5; - - ber_awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNodB/10))); - - sim_in.bps = 1; - awgn_bpsk = run_sim(sim_in); - sim_in.bps = 2; - awgn_qpsk = run_sim(sim_in); - - figure(5); clf; - semilogy(sim_in.EbNodB, ber_awgn_theory,'b+-;AWGN Theory;'); - hold on; - semilogy(sim_in.EbNodB, awgn_bpsk.ber,'g+-;Ns=7 BPSK;'); - semilogy(sim_in.EbNodB + awgn_bpsk.pilot_overhead, awgn_bpsk.ber,'go-;Ns=7 BPSK with pilot overhead;'); - semilogy(sim_in.EbNodB, awgn_qpsk.ber,'r+-;Ns=7 QPSK;'); - semilogy(sim_in.EbNodB + awgn_qpsk.pilot_overhead, awgn_qpsk.ber,'ro-;Ns=7 QPSK with pilot overhead;'); - hold off; - axis([-3 5 4E-3 2E-1]) - xlabel('Eb/No (dB)'); - ylabel('BER'); - grid; grid minor on; - legend('boxoff'); - title('AWGN'); -end - - -% HF multipath curves for BPSK and QPSK. Coded operating point is about 3dB - -function run_curves_hf_bpsk_qpsk - sim_in.Nc = 7; - sim_in.Ns = 7; - sim_in.Nsec = 120; - sim_in.verbose = 0; - sim_in.pilot_phase_est = 1; - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 0; - sim_in.stripped_phase_est = 0; - sim_in.phase_offset = 0; - sim_in.phase_test = 0; - sim_in.hf_en = 1; - sim_in.hf_phase = 1; - - sim_in.EbNodB = 1:8; - - EbNoLin = 10.^(sim_in.EbNodB/10); - hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); - - sim_in.bps = 1; - hf_bpsk = run_sim(sim_in); - sim_in.bps = 2; - hf_qpsk = run_sim(sim_in); - - figure(5); clf; - semilogy(sim_in.EbNodB, hf_theory,'b+-;HF Theory;'); - hold on; - semilogy(sim_in.EbNodB, hf_bpsk.ber,'g+-;Ns=7 BPSK;'); - semilogy(sim_in.EbNodB + hf_bpsk.pilot_overhead, hf_bpsk.ber,'go-;Ns=7 BPSK with pilot overhead;'); - semilogy(sim_in.EbNodB, hf_qpsk.ber,'r+-;Ns=7 QPSK;'); - semilogy(sim_in.EbNodB + hf_qpsk.pilot_overhead, hf_qpsk.ber,'ro-;Ns=7 QPSK with pilot overhead;'); - hold off; - axis([1 8 4E-3 2E-1]) - xlabel('Eb/No (dB)'); - ylabel('BER'); - grid; grid minor on; - legend('boxoff'); - title('HF Multipath'); -end - - -% AWGN curves for BPSK using 3 carrier 2D matrix pilot and ML pilot - -function run_curves_awgn_ml - sim_in.bps = 1; - sim_in.Nc = 7; - sim_in.Ns = 7; - sim_in.Nsec = 10; - sim_in.verbose = 0; - sim_in.pilot_phase_est = 1; - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 0; - sim_in.stripped_phase_est = 1; - sim_in.phase_offset = 0; - sim_in.phase_test = 0; - sim_in.hf_en = 0; - sim_in.hf_phase = 0; - - sim_in.EbNodB = -3:5; - - ber_awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNodB/10))); - - awgn_2d = run_sim(sim_in); - sim_in.pilot_phase_est = 0; - sim_in.ml_pd = 1; - awgn_ml = run_sim(sim_in); - - figure(5); clf; - semilogy(sim_in.EbNodB, ber_awgn_theory,'b+-;AWGN Theory;'); - hold on; - semilogy(sim_in.EbNodB, awgn_2d.ber,'g+-;Ns=7 3 carrier pilot BPSK;'); - semilogy(sim_in.EbNodB, awgn_ml.ber,'ro-;Ns=7 ML pilot BPSK;'); - hold off; - axis([-3 5 4E-3 5E-1]) - xlabel('Eb/No (dB)'); - ylabel('BER'); - grid; grid minor on; - legend('boxoff'); - title('AWGN'); -end - - -% HF multipath curves for ML - -function run_curves_hf_ml - sim_in.bps = 1; - sim_in.Nc = 7; - sim_in.Ns = 14; - sim_in.Nsec = 120; - sim_in.verbose = 0; - sim_in.pilot_phase_est = 1; - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 0; - sim_in.stripped_phase_est = 0; - sim_in.phase_offset = 0; - sim_in.phase_test = 0; - sim_in.hf_en = 1; - sim_in.hf_phase = 1; - - sim_in.EbNodB = 1:8; - - EbNoLin = 10.^(sim_in.EbNodB/10); - hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); - - hf_2d = run_sim(sim_in); - sim_in.pilot_phase_est = 0; - sim_in.ml_pd = 1; - hf_ml = run_sim(sim_in); - - figure(7); clf; - semilogy(sim_in.EbNodB, hf_theory,'b+-;HF Theory;'); - hold on; - semilogy(sim_in.EbNodB, hf_2d.ber,'g+-;Ns=7 3 carrier pilot BPSK;'); - semilogy(sim_in.EbNodB, hf_ml.ber,'ro-;Ns=7 ML pilot BPSK;'); - hold off; - axis([1 8 4E-3 2E-1]) - xlabel('Eb/No (dB)'); - ylabel('BER'); - grid; grid minor on; - legend('boxoff'); - title('HF Multipath'); -end - - - -function run_single - sim_in.bps = 1; - sim_in.Nsec = 30; - sim_in.Nc = 3; - sim_in.Ns = 9; - sim_in.EbNodB = -1; - sim_in.verbose = 1; - sim_in.pilot_phase_est = 1; - sim_in.pilot_wide = 1; - sim_in.pilot_interp = 0; - sim_in.stripped_phase_est = 0; - sim_in.phase_offset = 0; - sim_in.phase_test = 0; - sim_in.hf_en = 0; - sim_in.hf_phase = 0; - sim_in.ml_pd = 1; - - run_sim(sim_in); -end - - -format; -more off; - -%run_single -%run_curves_hf_bpsk_qpsk -%run_curves_hf_channels -run_curves_hf_ml - - - - diff --git a/codec2-dev/octave/ofdm_fs.m b/codec2-dev/octave/ofdm_fs.m new file mode 100644 index 00000000..c265adc7 --- /dev/null +++ b/codec2-dev/octave/ofdm_fs.m @@ -0,0 +1,771 @@ +% ofdm_fs.m +% David Rowe Mar 2017 +% +% Rate Fs BPSK/QPSK OFDM simulation, rate Fs verison of ofdm_rs.m with +% OFDM based up and down conversion. + +#{ + TODO: + [X] strip back experimental stuff to just features we need + [X] ZOH/integrator + [X] OFDM up and down conversion + [X] rate Fs HF model and HF results + [X] add QPSK + [X] add CP + [X] fine timing estimator and sample clock offset tracking + [X] acquisition coarse timing & freq offset estimation + [ ] adjust waveform parameters for real world + [ ] Nsec run time take into account CP + [ ] plot phase ests + [ ] handle border carriers + [ ] start with phantom carriers + + but unhappy with 1800Hz bandwidth + [ ] also try interpolation or just single row + [ ] compute SNR and PAPR + [ ] SSB bandpass filtering +#} + +1; + +% 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) + bit0 = real(symbol*exp(j*pi/4)) < 0; + bit1 = imag(symbol*exp(j*pi/4)) < 0; + two_bits = [bit1 bit0]; +endfunction + + +% Correlates the OFDM pilot symbol samples with a window of received +% samples to determine the most likely timing offset. Combines two +% frames pilots so we need at least Nsamperframe+M+Ncp samples in rx. +% Also determines frequency offset at maximimum correlation. Can be +% used for acquisition (coarse timing a freq offset), and fine timing + +function [t_est foff_est] = coarse_sync(states, rx, rate_fs_pilot_samples) + Nsamperframe = states.Nsamperframe; Fs = states.Fs; + Npsam = length(rate_fs_pilot_samples); + verbose = states.verbose; + + Ncorr = length(rx) - (Nsamperframe+Npsam) + 1; + assert(Ncorr > 0); + corr = zeros(1,Ncorr); + for i=1:Ncorr + corr(i) = abs(rx(i:i+Npsam-1) * rate_fs_pilot_samples'); + corr(i) += abs(rx(i+Nsamperframe:i+Nsamperframe+Npsam-1) * rate_fs_pilot_samples'); + end + + [mx t_est] = max(abs(corr)); + + C = abs(fft(rx(t_est:t_est+Npsam-1) .* conj(rate_fs_pilot_samples), Fs)); + C += abs(fft(rx(t_est+Nsamperframe:t_est+Nsamperframe+Npsam-1) .* conj(rate_fs_pilot_samples), Fs)); + + fmax = 30; + [mx_pos foff_est_pos] = max(C(1:fmax)); + [mx_neg foff_est_neg] = max(C(Fs-fmax+1:Fs)); + + if mx_pos > mx_neg + foff_est = foff_est_pos - 1; + else + foff_est = foff_est_neg - fmax - 1; + end + + if verbose > 1 + %printf("t_est: %d\n", t_est); + figure(7); clf; + plot(abs(corr)) + figure(8) + plot(C) + axis([0 200 0 0.4]) + end +endfunction + + +% init function +% load function +% default est on, but way to optionally disable + +#{ + Frame has Ns-1 data symbols between pilots, e.g. for Ns=3: + + PPP + DDD + DDD + PPP +#} + +function states = ofdm_init(bps, Rs, Tcp, Ns, Nc) + states.Fs = 8000; + states.bps = bps; + states.Rs = Rs; + states.Tcp = Tcp; + states.Ns = Ns; % step size for pilots + states.Nc = Nc; % Number of cols, aka number of carriers + states.M = states.Fs/Rs; + states.Ncp = Tcp*states.Fs; + states.Nbitsperframe = (Ns-1)*Nc*bps; + states.Nrowsperframe = states.Nbitsperframe/(Nc*bps); + states.Nsamperframe = (states.Nrowsperframe+1)*(states.M+states.Ncp); + + % generate same pilots each time + + rand('seed',1); + states.pilots = 1 - 2*(rand(1,Nc+2) > 0.5); +endfunction + + +function [sim_out rate_fs_pilot_samples rx] = run_sim(sim_in) + + % set up core modem constants + + states = ofdm_init(sim_in.bps, sim_in.Rs, sim_in.Tcp, sim_in.Ns, sim_in.Nc); + ofdm_load_const; + + % simulation parameters and flags + + woffset = 2*pi*sim_in.foff_hz/Fs; + EbNodB = sim_in.EbNodB; + verbose = states.verbose = sim_in.verbose; + hf_en = sim_in.hf_en; + timing_en = sim_in.timing_en; + foff_est_en = sim_in.foff_est_en; + phase_est_en = sim_in.phase_est_en; + + if verbose + printf("Rs:.........: %4.2f\n", Rs); + printf("M:..........: %d\n", M); + printf("bps:.........: %d\n", bps); + printf("Nbitsperframe: %d\n", Nbitsperframe); + printf("Nrowsperframe: %d\n", Nrowsperframe); + end + + % 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 + + Nrows = sim_in.Nsec*Rs; + Nframes = floor((Nrows-1)/Ns); + Nbits = Nframes * Nbitsperframe; % number of payload data bits + + Nr = Nbits/(Nc*bps); % Number of data rows to get Nbits total + + % double check if Nbits fit neatly into carriers + + assert(Nbits/(Nc*bps) == floor(Nbits/(Nc*bps)), "Nbits/(Nc*bps) must be an integer"); + + Nrp = Nr + Nframes + 1; % number of rows once pilots inserted + % extra row of pilots at end + + 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); + printf("Nframes: %d\n", Nframes); + printf("Nrp....: %d (number of rows including pilots)\n", Nrp); + end + + % set up HF model --------------------------------------------------------------- + + if hf_en + + % some typical values, or replace with user supplied + + dopplerSpreadHz = 1.0; path_delay_ms = 1; + + if isfield(sim_in, "dopplerSpreadHz") + dopplerSpreadHz = sim_in.dopplerSpreadHz; + end + if isfield(sim_in, "path_delay_ms") + path_delay_ms = sim_in.path_delay_ms; + end + path_delay_samples = path_delay_ms*Fs/1000; + printf("Doppler Spread: %3.2f Hz Path Delay: %3.2f ms %d samples\n", dopplerSpreadHz, path_delay_ms, path_delay_samples); + + % generate same fading pattern for every run + + randn('seed',1); + + spread1 = doppler_spread(dopplerSpreadHz, Fs, (sim_in.Nsec*(M+Ncp)/M+0.2)*Fs); + spread2 = doppler_spread(dopplerSpreadHz, Fs, (sim_in.Nsec*(M+Ncp)/M+0.2)*Fs); + + % sometimes doppler_spread() doesn't return exactly the number of samples we need + + assert(length(spread1) >= Nrp*M, "not enough doppler spreading samples"); + assert(length(spread2) >= Nrp*M, "not enough doppler spreading samples"); + + hf_gain = 1.0/sqrt(var(spread1)+var(spread2)); + % printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1)); + end + + % init timing est states + + tstates.Nsamperframe = Nsamperframe; tstates.M = M; tstates.Ncp = Ncp; + tstates.verbose = 0; tstates.Fs = Fs; + delta_t = []; + window_width = 11; % search window_width/2 from current timing instant + + % simulate for each Eb/No point ------------------------------------ + + for nn=1:length(EbNodB) + rand('seed',1); + randn('seed',1); + + EbNo = bps * (10 .^ (EbNodB(nn)/10)); + variance = 1/(M*EbNo/2); + + % generate tx bits + + tx_bits = rand(1,Nbits) > 0.5; + + % map to symbols in linear array + + if bps == 1 + tx_sym_lin = 2*tx_bits - 1; + end + if bps == 2 + for s=1:Nbits/bps + tx_sym_lin(s) = qpsk_mod(tx_bits(2*(s-1)+1:2*s)); + end + end + + % place symbols in multi-carrier frame with pilots and boundary carriers + + tx_sym = []; s = 1; + for f=1:Nframes + aframe = zeros(Nrowsperframe,Nc+2); + aframe(1,:) = pilots; + for r=1:Nrowsperframe + arowofsymbols = tx_sym_lin(s:s+Nc-1); + s += Nc; + aframe(r+1,2:Nc+1) = arowofsymbols; + end + tx_sym = [tx_sym; aframe]; + end + tx_sym = [tx_sym; pilots]; % final row of pilots + [nr nc] = size(tx_sym); + assert(nr == Nrp); + + % OFDM up conversion and upsampling to rate Fs --------------------- + + w = (0:Nc+1)*2*pi*Rs/Fs; + W = zeros(Nc+2,M); + for c=1:Nc+2 + W(c,:) = exp(j*w(c)*(0:M-1)); + end + + Nsam = Nrp*(M+Ncp); + tx = []; + + % OFDM upconvert symbol by symbol so we can add CP + + for r=1:Nrp + asymbol = tx_sym(r,:) * W/M; + asymbol_cp = [asymbol(M-Ncp+1:M) asymbol]; + tx = [tx asymbol_cp]; + end + assert(length(tx) == Nsam); + + % OFDM symbol of pilots is used for coarse timing and freq during acquisition, and fine timing + + rate_fs_pilot_samples = tx(1:M+Ncp); + + % channel simulation --------------------------------------------------------------- + + if isfield(sim_in, "sample_clock_offset_ppm") + + if sim_in.sample_clock_offset_ppm + timebase = floor(abs(1E6/sim_in.sample_clock_offset_ppm)); + if sim_in.sample_clock_offset_ppm > 0 + tx = resample(tx, timebase+1, timebase); + else + tx = resample(tx, timebase, timebase+1); + end + + % make sure length is correct for rest of simulation + + tx = [tx zeros(1,Nsam-length(tx))]; + tx = tx(1:Nsam); + end + end + + rx = tx; + + if hf_en + %rx = [zeros(1,path_delay_samples) tx(1:Nsam-path_delay_samples)]; + rx = hf_gain * tx(1:Nsam) .* spread1(1:Nsam); + rx += hf_gain * [zeros(1,path_delay_samples) tx(1:Nsam-path_delay_samples)] .* spread2(1:Nsam); + end + + rx = rx .* exp(j*woffset*(1:Nsam)); + + noise = sqrt(variance)*(0.5*randn(1,Nsam) + j*0.5*randn(1,Nsam)); + rx += noise; + + % some spare samples at end to allow for timing est window + + rx = [rx zeros(1,Nsamperframe)]; + + % pilot based phase est, we use known tx symbols as pilots ---------- + + rx_sym = zeros(Nrp, Nc+2); + 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); + timing_est_log = []; + timing_est = 1; + if timing_en + sample_point = timing_est; + else + sample_point = Ncp; + end + foff_est_hz_log = []; + foff_est_hz = 0; + foff_est_gain = 0.1; + Nerrs_log = []; + + for r=1:Ns:Nrp-Ns + + % printf("symbol r: %d\n", r); + + woff_est = 2*pi*foff_est_hz/Fs; + + if timing_en + + % update timing every frame + + if (mod(r-1,Ns) == 0) && (r != 1) && (r != Nrp) + st = (r-1)*(M+Ncp) + 1 - floor(window_width/2) + (timing_est-1); + en = st + Nsamperframe-1 + M+Ncp + window_width-1; + ft_est = coarse_sync(states, rx(st:en) .* exp(-j*woff_est*(st:en)), rate_fs_pilot_samples); + timing_est += ft_est - ceil(window_width/2); + if verbose > 1 + printf(" ft_est: %2d timing_est: %2d sample_point: %2d\n", ft_est, timing_est, sample_point); + %printf(" r : %d st: %d en: %d\n", r, st, en); + end + delta_t = [delta_t ft_est - ceil(window_width/2)]; + sample_point = max(timing_est+Ncp/4, sample_point); + sample_point = min(timing_est+Ncp, sample_point); + end + end + + timing_est_log = [timing_est_log timing_est]; + + % down convert at current timing instant + + % previous pilot + + if r >= Ns+1 + rr = r-Ns; + st = (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1; + for c=1:Nc+2 + acarrier = rx(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); + rx_sym(rr,c) = sum(acarrier); + end + %printf(" rr: %d st: %d en: %d\n", rr, st, en); + end + + % pilot - this frame - pilot + + for rr=r:r+Ns + st = (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1; + for c=1:Nc+2 + acarrier = rx(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); + rx_sym(rr,c) = sum(acarrier); + end + %printf(" rr: %d st: %d en: %d\n", rr, st, en); + end + + % next pilot + + if r <= Nrp - 2*Ns + rr = r+2*Ns; + st = (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1; + for c=1:Nc+2 + acarrier = rx(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); + rx_sym(rr,c) = sum(acarrier); + end + % printf(" rr: %d st: %d en: %d\n", rr, st, en); + end + + % est freq err based on all carriers + + if foff_est_en + freq_err_rect = sum(rx_sym(r,:))' * sum(rx_sym(r+Ns,:)); + freq_err_hz = angle(freq_err_rect)*Rs/(2*pi*Ns); + foff_est_hz += foff_est_gain*freq_err_hz; + end + foff_est_hz_log = [foff_est_hz_log foff_est_hz]; + + % OK - now estimate and correct phase + + for c=2:Nc+1 + + % 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; + if phase_est_en + rx_corr(rr,c) = rx_sym(rr,c) * exp(-j*aphase_est_pilot); + else + rx_corr(rr,c) = rx_sym(rr,c); + end + end + + end % c=2:Nc+1 + end % r=1:Ns:Nrp-Ns + + + % remove pilots to give us just data symbols and demodulate + + rx_bits = []; rx_np = []; + for r=1:Nrp + if mod(r-1,Ns) != 0 + arowofsymbols = rx_corr(r,2:Nc+1); + rx_np = [rx_np arowofsymbols]; + if bps == 1 + arowofbits = real(arowofsymbols) > 0; + end + if bps == 2 + arowofbits = zeros(1,Nc); + for c=1:Nc + arowofbits((c-1)*2+1:c*2) = qpsk_demod(arowofsymbols(c)); + end + end + rx_bits = [rx_bits arowofbits]; + end + end + assert(length(rx_bits) == Nbits); + + + % calculate BER stats as a block, after pilots extracted + + errors = xor(tx_bits, rx_bits); + Nerrs = sum(errors); + for f=1:Nframes + st = (f-1)*Nbitsperframe+1; en = st + Nbitsperframe-1; + Nerrs_log(f) = sum(xor(tx_bits(st:en), rx_bits(st:en))); + end + + printf("EbNodB: %3.2f BER: %5.4f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/Nbits, Nbits, Nerrs); + + if verbose + + figure(1); clf; + plot(rx_np,'+'); + axis([-2 2 -2 2]); + title('Scatter'); + + figure(2); clf; + plot(phase_est_log(:,2:Nc+1),'+', 'markersize', 10); + hold on; + plot(phase_est_pilot_log(:,2:Nc+1),'g+', 'markersize', 5); + title('Phase est'); + axis([1 Nrp -pi pi]); + + figure(3); clf; + subplot(211) + stem(delta_t) + title('Fine Timing'); + subplot(212) + plot(timing_est_log); + + figure(4); clf; + plot(foff_est_hz_log) + axis([1 max(Nframes,2) -3 3]); + title('Fine Freq'); + + figure(5); clf; + plot(Nerrs_log); + +#{ + figure(2) + Tx = abs(fft(tx(1:Nsam).*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]) +#} + +#{ + if hf_en + figure(4); clf; + subplot(211) + plot(abs(spread1(1:Nsam))); + %hold on; plot(abs(spread2(1:Nsam)),'g'); hold off; + subplot(212) + plot(angle(spread1(1:Nsam))); + title('spread1 amp and phase'); + end +#} + +#{ + % todo, work out a way to plot rate Fs hf model phase + if sim_in.hf_en + plot(angle(hf_model(:,2:Nc+1))); + end +#} + + + end + + sim_out.ber(nn) = sum(Nerrs)/Nbits; + sim_out.pilot_overhead = 10*log10(Ns/(Ns-1)); + sim_out.M = M; sim_out.Fs = Fs; sim_out.Ncp = Ncp; + sim_out.Nrowsperframe = Nrowsperframe; sim_out.Nsamperframe = Nsamperframe; + end +endfunction + + +function run_single + Ts = 0.018; + sim_in.Tcp = 0.002; + sim_in.Rs = 1/Ts; sim_in.bps = 2; sim_in.Nc = 16; sim_in.Ns = 8; + + %sim_in.Nsec = 1*(sim_in.Ns+1)/sim_in.Rs; % one frame + sim_in.Nsec = 10; + + sim_in.EbNodB = 6; + sim_in.verbose = 1; + sim_in.hf_en = 1; + sim_in.foff_hz = 0; + sim_in.timing_en = 1; + sim_in.sample_clock_offset_ppm = 0; + sim_in.foff_est_en = 1; + sim_in.phase_est_en = 1; + + run_sim(sim_in); +end + + +% Plot BER against Eb/No curves for AWGN and HF + +% Target operating point Eb/No for HF 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 +% +% For AWGN target is 2dB so -1dB op point. + +function run_curves + Ts = 0.010; + sim_in.Rs = 1/Ts; + sim_in.Tcp = 0.002; + sim_in.bps = 2; sim_in.Ns = 8; sim_in.Nc = 8; sim_in.verbose = 0; + sim_in.foff_hz = 0; + + pilot_overhead = (sim_in.Ns-1)/sim_in.Ns; + cp_overhead = Ts/(Ts+sim_in.Tcp); + overhead_dB = -10*log10(pilot_overhead*cp_overhead); + + sim_in.hf_en = 0; + sim_in.Nsec = 30; + sim_in.EbNodB = -3:5; + awgn_EbNodB = sim_in.EbNodB; + + awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNodB/10))); + awgn = run_sim(sim_in); + + sim_in.hf_en = 1; + sim_in.Nsec = 120; + sim_in.EbNodB = 1:8; + + EbNoLin = 10.^(sim_in.EbNodB/10); + hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); + + hf = run_sim(sim_in); + + figure(4); clf; + semilogy(awgn_EbNodB, awgn_theory,'b+-;AWGN theory;'); + hold on; + semilogy(sim_in.EbNodB, hf_theory,'b+-;HF theory;'); + semilogy(awgn_EbNodB+overhead_dB, awgn_theory,'g+-;AWGN lower bound with pilot + CP overhead;'); + semilogy(sim_in.EbNodB+overhead_dB, hf_theory,'g+-;HF lower bound with pilot + CP overhead;'); + semilogy(awgn_EbNodB+overhead_dB, awgn.ber,'r+-;AWGN sim;'); + semilogy(sim_in.EbNodB+overhead_dB, hf.ber,'r+-;HF sim;'); + hold off; + axis([-3 8 1E-2 2E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); +end + + +% Run an acquisition test, returning vectors of estimation errors + +function [delta_t delta_foff] = acquisition_test(Ntests=10, EbNodB=100, foff_hz=0, hf_en=0, fine_en=0) + + % generate test signal at a given Eb/No and frequency offset + + Ts = 0.016; + sim_in.Tcp = 0.002; + sim_in.Rs = 1/Ts; sim_in.bps = 2; sim_in.Nc = 16; sim_in.Ns = 8; + + sim_in.Nsec = Ntests*(sim_in.Ns+1)/sim_in.Rs; + + sim_in.EbNodB = EbNodB; + sim_in.verbose = 0; + sim_in.hf_en = hf_en; + sim_in.foff_hz = foff_hz; sim_in.timing_en = 0; + + % set up acquistion + + Nsamperframe = states.Nsamperframe = sim_out.Nsamperframe; + states.M = sim_out.M; states.Ncp = sim_out.Ncp; + states.verbose = 0; + states.Fs = sim_out.Fs; + + % test fine or acquisition over test signal + #{ + fine: - start with coarse timing instant + - on each frame est timing a few samples about that point + - update timing instant + + corr: - where is best plcase to sample + - just before end of symbol? + - how long should sequence be? + - add extra? + - aim for last possible moment? + - man I hope IL isn't too big..... + #} + + delta_t = []; delta_t = []; delta_foff = []; + + if fine_en + + window_width = 5; % search +/-2 samples from current timing instant + timing_instant = Nsamperframe+1; % start at correct instant for AWGN + % start at second frame so we can search -2 ... +2 + + while timing_instant < (length(rx) - (Nsamperframe + length(rate_fs_pilot_samples) + window_width)) + st = timing_instant - ceil(window_width/2); + en = st + Nsamperframe-1 + length(rate_fs_pilot_samples) + window_width-1; + [ft_est foff_est] = coarse_sync(states, rx(st:en), rate_fs_pilot_samples); + printf("ft_est: %d timing_instant %d %d\n", ft_est, timing_instant, mod(timing_instant, Nsamperframe)); + timing_instant += Nsamperframe + ft_est - ceil(window_width/2); + delta_t = [delta_ft ft_est - ceil(window_width/2)]; + end + else + % for coarse simulation we just use contant window shifts + + st = 0.5*Nsamperframe; + en = 2.5*Nsamperframe - 1; + ct_target = Nsamperframe/2; + + for w=1:Nsamperframe:length(rx)-3*Nsamperframe + %st = w+0.5*Nsamperframe; en = st+2*Nsamperframe-1; + %[ct_est foff_est] = coarse_sync(states, rx(st:en), rate_fs_pilot_samples); + [ct_est foff_est] = coarse_sync(states, rx(w+st:w+en), rate_fs_pilot_samples); + printf("ct_est: %4d foff_est: %3.1f\n", ct_est, foff_est); + + % valid coarse timing ests are modulo Nsamperframe + + delta_t = [delta_ct ct_est-ct_target]; + delta_foff = [delta_foff (foff_est-foff_hz)]; + end + end + +endfunction + + +% Run some tests for various acquisition conditions. Probability of +% acquistion is what matters, e.g. if it's 50% we can expect sync +% within 2 frames +% P(t)/P(f) P(t)/P(f) +% Eb/No AWGN HF +% +/- 25Hz -1/3 1.0/0.3 0.96/0.3 +% +/- 5Hz -1/3 1.0/0.347 0.96/0.55 +% +/- 25Hz 10/10 1.00/0.92 0.99/0.77 + +function acquisition_histograms(fine_en = 0) + Fs = 8000; + Ntests = 100; + + % allowable tolerance for acquistion + + ftol_hz = 2.0; + ttol_samples = 0.002*Fs; + + % AWGN channel operating point + + [dct dfoff] = acquisition_test(Ntests, -1, 25, 0, fine_en); + + % Probability of acquistion is what matters, e.g. if it's 50% we can + % expect sync within 2 frames + + printf("AWGN P(time offset acq) = %3.2f\n", length(find (abs(dct) < ttol_samples))/length(dct)); + if fine_en == 0 + printf("AWGN P(freq offset acq) = %3.2f\n", length(find (abs(dfoff) < ftol_hz))/length(dfoff)); + end + + figure(1) + hist(dct(find (abs(dct) < ttol_samples))) + if fine_en == 0 + figure(2) + hist(dfoff) + end + + % HF channel operating point + + [dct dfoff] = acquisition_test(Ntests, 3, 25, 1, fine_en); + + printf("HF P(time offset acq) = %3.2f\n", length(find (abs(dct) < ttol_samples))/length(dct)); + if fine_en == 0 + printf("HF P(freq offset acq) = %3.2f\n", length(find (abs(dfoff) < ftol_hz))/length(dfoff)); + end + + figure(3) + hist(dct(find (abs(dct) < ttol_samples))) + if fine_en == 0 + figure(4) + hist(dfoff) + end + +endfunction + + + +% choose simulation to run here ------------------------------------------------------- + +format; +more off; + +run_single +%run_curves +%acquisition_histograms(1) + diff --git a/codec2-dev/octave/ofdm_load_const.m b/codec2-dev/octave/ofdm_load_const.m new file mode 100644 index 00000000..8b355623 --- /dev/null +++ b/codec2-dev/octave/ofdm_load_const.m @@ -0,0 +1,15 @@ +% make like C #define for ofdm modem + +Fs = states.Fs; +bps = states.bps; +Rs = states.Rs; +Tcp = states.Tcp; +Ns = states.Ns; +Nc = states.Nc; +M = states.M; +Ncp = states.Ncp; +bps = sim_in.bps; +Nbitsperframe = states.Nbitsperframe; +Nrowsperframe = states.Nrowsperframe; +Nsamperframe = states.Nsamperframe; +pilots = states.pilots; diff --git a/codec2-dev/octave/ofdm_rs.m b/codec2-dev/octave/ofdm_rs.m new file mode 100644 index 00000000..9b667310 --- /dev/null +++ b/codec2-dev/octave/ofdm_rs.m @@ -0,0 +1,788 @@ +% ofdm_rs.m +% David Rowe Mar 2017 +% +% Rate Rs BPSK/QPSK simulation to explore techniques for +% phase estimation over multiple carriers in HF channel. Rate +% Rs version of ofdm_fs.m + +#{ + TODO: + [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 + [X] test single points + + 1dB IL @ 6dB HF, 0.4 dB @ 2dB AWGN + [X] try linear interpolation + [X] try longer time windows + [X] try combining mod stripping phase est inside frame + [X] curves taking into account pilot losses + [ ] remove border carriers, interpolate edge carrier + [X] modify for QPSK +#} + +1; + +% 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) + bit0 = real(symbol*exp(j*pi/4)) < 0; + bit1 = imag(symbol*exp(j*pi/4)) < 0; + two_bits = [bit1 bit0]; +endfunction + + +function sim_out = run_sim(sim_in) + Rs = 100; + + bps = sim_in.bps; + EbNodB = sim_in.EbNodB; + verbose = sim_in.verbose; + hf_en = sim_in.hf_en; + hf_phase = sim_in.hf_phase; + phase_offset = sim_in.phase_offset; + + Ns = sim_in.Ns; % step size for pilots + Nc = sim_in.Nc; % Number of cols, aka number of carriers + + Nbitsperframe = (Ns-1)*Nc*bps; + Nrowsperframe = Nbitsperframe/(Nc*bps); + if verbose + printf("bps:.........: %d\n", bps); + printf("Nbitsperframe: %d\n", Nbitsperframe); + printf("Nrowsperframe: %d\n", Nrowsperframe); + end + + % 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*bps); % 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*bps) == floor(Nbits/(Nc*bps)), "Nbits/(Nc*bps) 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 + + % some typical values, or replace with user supplied + + dopplerSpreadHz = 1.0; path_delay = 1E-3*Rs; + + if isfield(sim_in, "dopplerSpreadHz") + dopplerSpreadHz = sim_in.dopplerSpreadHz; + end + if isfield(sim_in, "path_delay") + path_delay = sim_in.path_delay; + end + printf("Doppler Spread: %3.2f Hz Path Delay: %3.2f symbols\n", dopplerSpreadHz, path_delay); + randn('seed',1); + spread1 = doppler_spread(dopplerSpreadHz, Rs, sim_in.Nsec*Rs*1.1); + spread2 = doppler_spread(dopplerSpreadHz, Rs, sim_in.Nsec*Rs*1.1); + + % 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 + + % 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); + randn('seed',1); + + EsNo = bps * (10 .^ (EbNodB(nn)/10)); + variance = 1/(EsNo/2); + noise = sqrt(variance)*(0.5*randn(Nrp,Nc+2) + j*0.5*randn(Nrp,Nc+2)); + + % generate tx bits + + tx_bits = rand(1,Nbits) > 0.5; + + % map to symbols + + if bps == 1 + tx_symb = 2*tx_bits - 1; + end + if bps == 2 + for s=1:Nbits/bps + tx_symb(s) = qpsk_mod(tx_bits(2*(s-1)+1:2*s)); + end + end + + % place symbols in multi-carrier frame with pilots and boundary carriers + + tx = []; s = 1; + for f=1:Nframes + aframe = zeros(Nrowsperframe,Nc+2); + aframe(1,:) = 1; + for r=1:Nrowsperframe + arowofsymbols = tx_symb(s:s+Nc-1); + s += Nc; + aframe(r+1,2:Nc+1) = arowofsymbols; + end + tx = [tx; aframe]; + end + tx = [tx; ones(1,Nc+2)]; % final row of pilots + [nr nc] = size(tx); + assert(nr == Nrp); + + 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. + + % Note Rs carrier spacing, sample rate is Rs + + 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(r,c) = hf_gain*(spread1(r) + exp(-j*w*path_delay)*spread2(r)); + end + + if hf_phase + rx(r,:) = rx(r,:) .* hf_model(r,:); + else + rx(r,:) = rx(r,:) .* abs(hf_model(r,:)); + end + end + end + + rx += noise; + + % pilot based phase est, we use known tx symbols as pilots ---------- + + rx_corr = rx; + + if sim_in.pilot_phase_est + + % est phase from pilots either side of data symbols + % adjust phase of data symbol + % demodulate and count errors of just data + + 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_rect1 = sum(rx(r,cr)*tx(r,cr)'); + aphase_est_pilot_rect2 = sum(rx(r+Ns,cr)*tx(r+Ns,cr)'); + + % optionally use next step of pilots in past and future + + if sim_in.pilot_wide + if r > Ns+1 + aphase_est_pilot_rect1 += sum(rx(r-Ns,cr)*tx(r-Ns,cr)'); + end + if r < Nrp - 2*Ns + aphase_est_pilot_rect2 += sum(rx(r+2*Ns,cr)*tx(r+2*Ns,cr)'); + end + end + + % correct phase offset using phase estimate + + for rr=r+1:r+Ns-1 + a = b = 1; + if sim_in.pilot_interp + b = (rr-r)/Ns; a = 1 - b; + end + %printf("rr: %d a: %4.3f b: %4.3f\n", rr, a, b); + aphase_est_pilot = angle(a*aphase_est_pilot_rect1 + b*aphase_est_pilot_rect2); + phase_est_pilot_log(rr,c) = aphase_est_pilot; + rx_corr(rr,c) = rx(rr,c) * exp(-j*aphase_est_pilot); + end + + if sim_in.stripped_phase_est + % Optional modulation stripping feed fwd phase estimation, to refine + % pilot-based phase estimate. Doing it after pilot based phase estimation + % means we don't need to deal with ambiguity, which is difficult to handle + % in low SNR channels. + + % Use vector of 7 symbols around current data symbol. We could use a 2D + % window if we can work out how best to correct with pilot-est and avoid + % ambiguities + + for rr=r+1:r+Ns-1 + + % extract a matrix of nearby samples with pilot-based offset removed + + amatrix = rx(max(1,rr-3):min(Nrp,rr+3),c) .* exp(-j*aphase_est_pilot); + + % modulation strip and est phase + + stripped = abs(amatrix) .* exp(j*2*angle(amatrix)); + aphase_est_stripped = angle(sum(sum(stripped)))/2; + phase_est_stripped_log(rr,c) = aphase_est_stripped; + + % correct rx symbols based on both phase ests + + phase_est_log(rr,c) = angle(exp(j*(aphase_est_pilot+aphase_est_stripped))); + rx_corr(rr,c) = rx(rr,c) * exp(-j*phase_est_log(rr,c)); + end + end % sim_in.stripped_phase_est + + end % r=1:Ns:Nrp-Ns + + end % c=2:Nc+1 + end % sim_in.pilot_phase_est + + + if isfield(sim_in, "ml_pd") && sim_in.ml_pd + + % Bill's ML with pilots phase detector, does phase est and demodulation + + rx_bits = []; rx_np = []; + aframeofbits = zeros(Ns-1, Nc); + for r=1:Ns:Nrp-Ns + + % demodulate this frame, ML operates carrier by carrier + + for c=2:Nc+1 + arxcol = rx(r:r+Ns, c); + arxcol(1) = rx(r, c-1) + rx(r, c+1); + arxcol(Ns+1) = rx(r+Ns, c-1) + rx(r+Ns, c+1); + [acolofbits aphase_est] = ml_pd(rot90(arxcol), bps, [1 Ns+1]); + aframeofbits(:,c-1) = xor(acolofbits, ones(1,Ns-1)); + rx_np = [rx_np rot90(arxcol) .* exp(-j*aphase_est)]; + end + + % unpack from frame into linear array of bits + + for rr=1:Ns-1 + rx_bits = [rx_bits aframeofbits(rr,:)]; + end + + end + else + + % remove pilots to give us just data symbols and demodulate + + rx_bits = []; rx_np = []; + for r=1:Nrp + if mod(r-1,Ns) != 0 + arowofsymbols = rx_corr(r,2:Nc+1); + rx_np = [rx_np arowofsymbols]; + if bps == 1 + arowofbits = real(arowofsymbols) > 0; + end + if bps == 2 + arowofbits = zeros(1,Nc); + for c=1:Nc + arowofbits((c-1)*2+1:c*2) = qpsk_demod(arowofsymbols(c)); + end + end + rx_bits = [rx_bits arowofbits]; + end + end + end + %tx_bits + %rx_bits + assert(length(rx_bits) == Nbits); + + %phase_test + %phase_est_log + + % calculate BER stats as a block, after pilots extracted + + errors = xor(tx_bits, rx_bits); + Nerrs = sum(errors); + + printf("EbNodB: %3.2f BER: %5.4f 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 hf_en + figure(2); clf; + plot(abs(hf_model(:,2:Nc+1))); + end + + if sim_in.pilot_phase_est + figure(3); 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.stripped_phase_est + plot(phase_est_stripped_log(:,2:Nc+1),'ro', 'markersize', 5); + end + if sim_in.hf_en && sim_in.hf_phase + plot(angle(hf_model(:,2:Nc+1))); + end + if sim_in.phase_test + plot(phase_test(:,2:Nc+1)); + end + axis([1 Nrp -pi pi]); + end + 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 = 7; + sim_in.Ns = 5; + sim_in.Nsec = 240; + sim_in.EbNodB = 1:8; + sim_in.verbose = 0; + sim_in.pilot_phase_est = 0; + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 0; + sim_in.stripped_phase_est = 0; + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 1; + sim_in.hf_phase = 0; + + sim_in.Ns = 5; + hf_ref_Ns_5_no_phase = run_sim(sim_in); + sim_in.Ns = 9; + hf_ref_Ns_9_no_phase = run_sim(sim_in); + + sim_in.hf_phase = 1; + sim_in.pilot_phase_est = 1; + + sim_in.Ns = 5; + hf_Ns_5 = run_sim(sim_in); + + sim_in.Ns = 9; + hf_Ns_9 = run_sim(sim_in); + + sim_in.Ns = 17; + hf_Ns_17 = run_sim(sim_in); + + figure(4); clf; + semilogy(sim_in.EbNodB, hf_ref_Ns_5_no_phase.ber,'b+-;Ns=5 HF ref no phase;'); + hold on; + semilogy(sim_in.EbNodB, hf_ref_Ns_9_no_phase.ber,'c+-;Ns=9 HF ref no phase;'); + semilogy(sim_in.EbNodB, hf_Ns_5.ber,'g+--;Ns=5;'); + semilogy(sim_in.EbNodB + hf_Ns_5.pilot_overhead, hf_Ns_5.ber,'go-;Ns=5 with pilot overhead;'); + semilogy(sim_in.EbNodB, hf_Ns_9.ber,'r+--;Ns=9;'); + semilogy(sim_in.EbNodB + hf_Ns_9.pilot_overhead, hf_Ns_9.ber,'ro-;Ns=9 with pilot overhead;'); + semilogy(sim_in.EbNodB, hf_Ns_17.ber,'k+--;Ns=17;'); + semilogy(sim_in.EbNodB + hf_Ns_17.pilot_overhead, hf_Ns_17.ber,'ko-;Ns=17 with pilot overhead;'); + hold off; + axis([1 8 4E-2 2E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); + title('HF Multipath 1Hz Doppler 1ms delay'); + +end + + +% Generate HF curves for some alternative, experimental methods tested +% during development, such as interpolation, refinements using +% modulation stripping, narrow window. + +function run_curves_hf_alt + sim_in.Nc = 7; + sim_in.Ns = 5; + sim_in.Nsec = 60; + sim_in.EbNodB = 1:8; + sim_in.verbose = 0; + sim_in.pilot_phase_est = 0; + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 0; + sim_in.stripped_phase_est = 0; + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 1; + sim_in.hf_phase = 0; + + sim_in.Ns = 9; + hf_ref_Ns_9_no_phase = run_sim(sim_in); + + sim_in.hf_phase = 1; + sim_in.pilot_phase_est = 1; + hf_Ns_9 = run_sim(sim_in); + + sim_in.stripped_phase_est = 1; + hf_Ns_9_stripped = run_sim(sim_in); + + sim_in.stripped_phase_est = 0; + sim_in.pilot_wide = 0; + hf_Ns_9_narrow = run_sim(sim_in); + + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 1; + hf_Ns_9_interp = run_sim(sim_in); + + figure(6); clf; + semilogy(sim_in.EbNodB, hf_ref_Ns_9_no_phase.ber,'c+-;Ns=9 HF ref no phase;'); + hold on; + semilogy(sim_in.EbNodB, hf_Ns_9.ber,'r+--;Ns=9;'); + semilogy(sim_in.EbNodB, hf_Ns_9_stripped.ber,'g+--;Ns=9 stripped refinement;'); + semilogy(sim_in.EbNodB, hf_Ns_9_narrow.ber,'b+--;Ns=9 narrow;'); + semilogy(sim_in.EbNodB, hf_Ns_9_interp.ber,'k+--;Ns=9 interp;'); + hold off; + axis([1 8 4E-2 2E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); + title('HF Multipath 1Hz Doppler 1ms delay'); + +end + + +% Generate HF curves for fixed Ns but different HF channels. + +function run_curves_hf_channels + sim_in.Nc = 7; + sim_in.Ns = 9; + sim_in.Nsec = 240; + sim_in.EbNodB = 1:8; + sim_in.verbose = 0; + sim_in.pilot_phase_est = 0; + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 0; + sim_in.stripped_phase_est = 0; + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 1; + sim_in.hf_phase = 0; + + hf_Ns_9_1hz_1ms_no_phase = run_sim(sim_in); + + sim_in.hf_phase = 1; + sim_in.pilot_phase_est = 1; + hf_Ns_9_1hz_1ms = run_sim(sim_in); + + Rs = 100; + + sim_in.dopplerSpreadHz = 1.0; + sim_in.path_delay = 500E-6*Rs; + hf_Ns_9_1hz_500us = run_sim(sim_in); + + sim_in.dopplerSpreadHz = 1.0; + sim_in.path_delay = 2E-3*Rs; + hf_Ns_9_1hz_2ms = run_sim(sim_in); + + sim_in.dopplerSpreadHz = 2.0; + sim_in.path_delay = 1E-3*Rs; + hf_Ns_9_2hz_1ms = run_sim(sim_in); + + sim_in.dopplerSpreadHz = 2.0; + sim_in.path_delay = 1E-3*Rs; + hf_Ns_9_2hz_2ms = run_sim(sim_in); + + sim_in.dopplerSpreadHz = 2.0; + sim_in.path_delay = 2E-3*Rs; + hf_Ns_9_2hz_2ms = run_sim(sim_in); + + sim_in.dopplerSpreadHz = 4.0; + sim_in.path_delay = 1E-3*Rs; + hf_Ns_9_4hz_1ms = run_sim(sim_in); + + figure(6); clf; + semilogy(sim_in.EbNodB, hf_Ns_9_1hz_1ms_no_phase.ber,'c+-;Ns=9 1Hz 1ms ref no phase;'); + hold on; + semilogy(sim_in.EbNodB, hf_Ns_9_1hz_500us.ber,'k+-;Ns=9 1Hz 500us;'); + semilogy(sim_in.EbNodB, hf_Ns_9_1hz_1ms.ber,'r+-;Ns=9 1Hz 1ms;'); + semilogy(sim_in.EbNodB, hf_Ns_9_1hz_2ms.ber,'bo-;Ns=9 1Hz 2ms;'); + semilogy(sim_in.EbNodB, hf_Ns_9_2hz_1ms.ber,'g+-;Ns=9 2Hz 1ms;'); + semilogy(sim_in.EbNodB, hf_Ns_9_2hz_2ms.ber,'mo-;Ns=9 2Hz 2ms;'); + semilogy(sim_in.EbNodB, hf_Ns_9_4hz_1ms.ber,'c+-;Ns=9 4Hz 1ms;'); + hold off; + axis([1 8 4E-2 2E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); + title('HF Multipath Ns = 9'); + +end + + +% AWGN curves for BPSK and QPSK. Coded Eb/No operating point is 2dB, +% so raw BER for rate 1/2 will be -1dB + +function run_curves_awgn_bpsk_qpsk + sim_in.Nc = 7; + sim_in.Ns = 7; + sim_in.Nsec = 30; + sim_in.verbose = 0; + sim_in.pilot_phase_est = 0; + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 0; + sim_in.stripped_phase_est = 1; + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 0; + sim_in.hf_phase = 0; + + sim_in.EbNodB = -3:5; + + ber_awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNodB/10))); + + sim_in.bps = 1; + awgn_bpsk = run_sim(sim_in); + sim_in.bps = 2; + awgn_qpsk = run_sim(sim_in); + + figure(5); clf; + semilogy(sim_in.EbNodB, ber_awgn_theory,'b+-;AWGN Theory;'); + hold on; + semilogy(sim_in.EbNodB, awgn_bpsk.ber,'g+-;Ns=7 BPSK;'); + semilogy(sim_in.EbNodB + awgn_bpsk.pilot_overhead, awgn_bpsk.ber,'go-;Ns=7 BPSK with pilot overhead;'); + semilogy(sim_in.EbNodB, awgn_qpsk.ber,'r+-;Ns=7 QPSK;'); + semilogy(sim_in.EbNodB + awgn_qpsk.pilot_overhead, awgn_qpsk.ber,'ro-;Ns=7 QPSK with pilot overhead;'); + hold off; + axis([-3 5 4E-3 2E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); + title('AWGN'); +end + + +% HF multipath curves for BPSK and QPSK. Coded operating point is about 3dB + +function run_curves_hf_bpsk_qpsk + sim_in.Nc = 7; + sim_in.Ns = 7; + sim_in.Nsec = 120; + sim_in.verbose = 0; + sim_in.pilot_phase_est = 1; + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 0; + sim_in.stripped_phase_est = 0; + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 1; + sim_in.hf_phase = 1; + + sim_in.EbNodB = 1:8; + + EbNoLin = 10.^(sim_in.EbNodB/10); + hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); + + sim_in.bps = 1; + hf_bpsk = run_sim(sim_in); + sim_in.bps = 2; + hf_qpsk = run_sim(sim_in); + + figure(5); clf; + semilogy(sim_in.EbNodB, hf_theory,'b+-;HF Theory;'); + hold on; + semilogy(sim_in.EbNodB, hf_bpsk.ber,'g+-;Ns=7 BPSK;'); + semilogy(sim_in.EbNodB + hf_bpsk.pilot_overhead, hf_bpsk.ber,'go-;Ns=7 BPSK with pilot overhead;'); + semilogy(sim_in.EbNodB, hf_qpsk.ber,'r+-;Ns=7 QPSK;'); + semilogy(sim_in.EbNodB + hf_qpsk.pilot_overhead, hf_qpsk.ber,'ro-;Ns=7 QPSK with pilot overhead;'); + hold off; + axis([1 8 4E-3 2E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); + title('HF Multipath'); +end + + +% AWGN curves for BPSK using 3 carrier 2D matrix pilot and ML pilot + +function run_curves_awgn_ml + sim_in.bps = 1; + sim_in.Nc = 7; + sim_in.Ns = 7; + sim_in.Nsec = 10; + sim_in.verbose = 0; + sim_in.pilot_phase_est = 1; + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 0; + sim_in.stripped_phase_est = 1; + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 0; + sim_in.hf_phase = 0; + + sim_in.EbNodB = -3:5; + + ber_awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNodB/10))); + + awgn_2d = run_sim(sim_in); + sim_in.pilot_phase_est = 0; + sim_in.ml_pd = 1; + awgn_ml = run_sim(sim_in); + + figure(5); clf; + semilogy(sim_in.EbNodB, ber_awgn_theory,'b+-;AWGN Theory;'); + hold on; + semilogy(sim_in.EbNodB, awgn_2d.ber,'g+-;Ns=7 3 carrier pilot BPSK;'); + semilogy(sim_in.EbNodB, awgn_ml.ber,'ro-;Ns=7 ML pilot BPSK;'); + hold off; + axis([-3 5 4E-3 5E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); + title('AWGN'); +end + + +% HF multipath curves for ML + +function run_curves_hf_ml + sim_in.bps = 1; + sim_in.Nc = 7; + sim_in.Ns = 14; + sim_in.Nsec = 120; + sim_in.verbose = 0; + sim_in.pilot_phase_est = 1; + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 0; + sim_in.stripped_phase_est = 0; + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 1; + sim_in.hf_phase = 1; + + sim_in.EbNodB = 1:8; + + EbNoLin = 10.^(sim_in.EbNodB/10); + hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); + + hf_2d = run_sim(sim_in); + sim_in.pilot_phase_est = 0; + sim_in.ml_pd = 1; + hf_ml = run_sim(sim_in); + + figure(7); clf; + semilogy(sim_in.EbNodB, hf_theory,'b+-;HF Theory;'); + hold on; + semilogy(sim_in.EbNodB, hf_2d.ber,'g+-;Ns=7 3 carrier pilot BPSK;'); + semilogy(sim_in.EbNodB, hf_ml.ber,'ro-;Ns=7 ML pilot BPSK;'); + hold off; + axis([1 8 4E-3 2E-1]) + xlabel('Eb/No (dB)'); + ylabel('BER'); + grid; grid minor on; + legend('boxoff'); + title('HF Multipath'); +end + + + +function run_single + sim_in.bps = 2; + sim_in.Nsec = 30; + sim_in.Nc = 16; + sim_in.Ns = 8; + sim_in.EbNodB = 4; + sim_in.verbose = 1; + + sim_in.pilot_phase_est = 0; + sim_in.pilot_wide = 1; + sim_in.pilot_interp = 0; + sim_in.stripped_phase_est = 0; + sim_in.ml_pd = 0; + + sim_in.phase_offset = 0; + sim_in.phase_test = 0; + sim_in.hf_en = 0; + sim_in.hf_phase = 1; + + run_sim(sim_in); +end + + +format; +more off; + +run_single +%run_curves_hf_bpsk_qpsk +%run_curves_hf_channels +%run_curves_hf_ml + + + +