+++ /dev/null
-% 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)
-
+++ /dev/null
-% 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
-
-
-
-
--- /dev/null
+% 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)
+
--- /dev/null
+% 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;
--- /dev/null
+% 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
+
+
+
+
projects / freetel-svn-tracking.git / commitdiff