From be9acf2659ac2aee76cee4bb18c20b5c6f2dc72c Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sat, 29 Apr 2017 22:40:11 +0000 Subject: [PATCH] split rate fs ofdm simulation into multiple files git-svn-id: https://svn.code.sf.net/p/freetel/code@3111 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/{ofdm_fs.m => ofdm_dev.m} | 402 +------------------- codec2-dev/octave/ofdm_lib.m | 373 ++++++++++++++++++ 2 files changed, 384 insertions(+), 391 deletions(-) rename codec2-dev/octave/{ofdm_fs.m => ofdm_dev.m} (57%) create mode 100644 codec2-dev/octave/ofdm_lib.m diff --git a/codec2-dev/octave/ofdm_fs.m b/codec2-dev/octave/ofdm_dev.m similarity index 57% rename from codec2-dev/octave/ofdm_fs.m rename to codec2-dev/octave/ofdm_dev.m index 821f9842..2cfed8af 100644 --- a/codec2-dev/octave/ofdm_fs.m +++ b/codec2-dev/octave/ofdm_dev.m @@ -1,391 +1,15 @@ -% ofdm_fs.m -% David Rowe Mar 2017 +% ofdm_dev.m +% David Rowe April 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 - - -#{ - 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); - - % carrier tables for up and down conversion - - w = (0:Nc+1)*2*pi*Rs/states.Fs; - W = zeros(Nc+2,states.M); - for c=1:Nc+2 - W(c,:) = exp(j*w(c)*(0:states.M-1)); - end - states.w = w; - states.W = W; - - % fine timing search +/- window_width/2 from current timing instant - - states.window_width = 11; - - % Receive buffer: D P DDD P DDD P DDD P D - % ^ - % also see ofdm_demod() ... - - states.Nrxbuf = 3*states.Nsamperframe+states.M+states.Ncp + 2*(states.M + states.Ncp); - states.rxbuf = zeros(1, states.Nrxbuf); - - % default settings on a bunch of options and states - - states.verbose = 0; - states.timing_en = 1; - states.foff_est_en = 1; - states.phase_est_en = 1; - - states.foff_est_gain = 0.1; - states.foff_est_hz = 0; - states.sample_point = states.timing_est = 1; - states.nin = states.Nsamperframe; - -endfunction - - -% --------------------------- -% Modulates one frame of bits -% --------------------------- - -function tx = ofdm_mod(states, tx_bits) - ofdm_load_const; - - assert(length(tx_bits) == Nbitsperframe); - - % map to symbols in linear array - - if bps == 1 - tx_sym_lin = 2*tx_bits - 1; - end - if bps == 2 - for s=1:Nbitsperframe/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; - aframe = zeros(Ns,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]; - - % OFDM upconvert symbol by symbol so we can add CP - - tx = []; - for r=1:Ns - asymbol = tx_sym(r,:) * W/M; - asymbol_cp = [asymbol(M-Ncp+1:M) asymbol]; - tx = [tx asymbol_cp]; - end -endfunction - - -% ----------------------------- -% Demodulates one frame of bits -% ----------------------------- - -#{ - - For phase estimation we need to maintain buffer of 3 frames plus - one pilot, so we have 4 pilots total. '^' is the start of current - frame that we are demodulating. - - P DDD P DDD P DDD P - ^ - - Then add one symbol either side to account for movement in - sampling instant due to sample clock differences: - - D P DDD P DDD P DDD P D - ^ -#} - -function [rx_bits states aphase_est_pilot_log rx_np] = ofdm_demod(states, rxbuf_in) - ofdm_load_const; - - % extra states that are st up at run time rather than init time - - timing_est = states.timing_est; - timing_en = states.timing_en; - foff_est_hz = states.foff_est_hz; - foff_est_gain = states.foff_est_gain; - foff_est_en = states.foff_est_en; - sample_point = states.sample_point; - rate_fs_pilot_samples = states.rate_fs_pilot_samples; - verbose = states.verbose; - phase_est_en = states.phase_est_en; - - % insert latest input samples into rxbuf - - rxbuf(1:Nrxbuf-states.nin) = rxbuf(states.nin+1:Nrxbuf); - rxbuf(Nrxbuf-states.nin+1:Nrxbuf) = rxbuf_in; - - woff_est = 2*pi*foff_est_hz/Fs; - - % update timing estimate -------------------------------------------------- - - delta_t = sample_point = 0; - if timing_en - % update timing at start of every frame - - st = M+Ncp + Nsamperframe + 1 - floor(window_width/2) + (timing_est-1); - en = st + Nsamperframe-1 + M+Ncp + window_width-1; - - ft_est = coarse_sync(states, rxbuf(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); - end - - % Black magic to keep sample_point inside cyclic prefix. Or something like that. - - 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 - - % down convert at current timing instant---------------------------------- - - % todo: this cld be more efficent, as pilot r becomes r-Ns on next frame - - rx_sym = zeros(1+Ns+1+1, Nc+2); - - % previous pilot - - st = M+Ncp + Nsamperframe + (-Ns)*(M+Ncp) + 1 + sample_point; en = st + M - 1; - - for c=1:Nc+2 - acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); - rx_sym(1,c) = sum(acarrier); - end - - % pilot - this frame - pilot - - for rr=1:Ns+1 - st = M+Ncp + Nsamperframe + (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1; - for c=1:Nc+2 - acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); - rx_sym(rr+1,c) = sum(acarrier); - end - end - - % next pilot - - st = M+Ncp + Nsamperframe + (2*Ns)*(M+Ncp) + 1 + sample_point; en = st + M - 1; - for c=1:Nc+2 - acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); - rx_sym(Ns+3,c) = sum(acarrier); - end - - % est freq err based on all carriers ------------------------------------ - - if foff_est_en - freq_err_rect = sum(rx_sym(2,:))' * sum(rx_sym(2+Ns,:)); - freq_err_hz = angle(freq_err_rect)*Rs/(2*pi*Ns); - foff_est_hz += foff_est_gain*freq_err_hz; - end - - % OK - now estimate and correct phase ---------------------------------- - - aphase_est_pilot = 10*ones(1,Nc+2); - 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(2,cr)*pilots(cr)') + sum(rx_sym(2+Ns,cr)*pilots(cr)'); - - % use next step of pilots in past and future - - aphase_est_pilot_rect += sum(rx_sym(1,cr)*pilots(cr)'); - aphase_est_pilot_rect += sum(rx_sym(2+Ns+1,cr)*pilots(cr)'); - - aphase_est_pilot(c) = angle(aphase_est_pilot_rect); - end - - % correct phase offset using phase estimate, and demodulate - % bits, separate loop as it runs across cols (carriers) to get - % frame bit ordering correct - - aphase_est_pilot_log = []; - rx_bits = []; rx_np = []; - for rr=1:Ns-1 - for c=2:Nc+1 - if phase_est_en - rx_corr = rx_sym(rr+2,c) * exp(-j*aphase_est_pilot(c)); - else - rx_corr = rx_sym(rr+2,c); - end - rx_np = [rx_np rx_corr]; - if bps == 1 - abit = real(rx_corr) > 0; - end - if bps == 2 - abit = qpsk_demod(rx_corr); - end - rx_bits = [rx_bits abit]; - end % c=2:Nc+1 - aphase_est_pilot_log = [aphase_est_pilot_log; aphase_est_pilot]; - end - - % Adjust nin to take care of sample clock offset - - nin = Nsamperframe; - if timing_en - thresh = (M+Ncp)/8; - tshift = (M+Ncp)/4; - if timing_est > thresh - nin = Nsamperframe+tshift; - timing_est -= tshift; - sample_point -= tshift; - end - if timing_est < -thresh - nin = Nsamperframe-tshift; - timing_est += tshift; - sample_point += tshift; - end - end - - states.rxbuf = rxbuf; - states.nin = nin; - states.timing_est = timing_est; - states.sample_point = sample_point; - states.delta_t = delta_t; - states.foff_est_hz = foff_est_hz; -endfunction +% Simulations used for development and testing of Rate Fs BPSK/QPSK +% OFDM modem. +ofdm_lib; #{ TODO: - [ ] Some states need warm rest at the start of each simulation point - [ ] rate_fs_pilot_samples generated in init - [ ] move buffer shift into demod + [ ] compute SNR and PAPR + [ ] SSB bandpass filtering [ ] way to simulate aquisition and demod [ ] testframe based #} @@ -505,11 +129,6 @@ function [sim_out rate_fs_pilot_samples rx] = run_sim(sim_in) tx = [tx tx(st:en)]; assert(length(tx) == Nsam); - % OFDM symbol of pilots is used for coarse timing and freq during acquisition, and fine timing - % TODO: put this in init code - - states.rate_fs_pilot_samples = tx(1:M+Ncp); - % channel simulation --------------------------------------------------------------- if isfield(sim_in, "sample_clock_offset_ppm") @@ -634,7 +253,7 @@ function [sim_out rate_fs_pilot_samples rx] = run_sim(sim_in) title('Fine Freq'); figure(5); clf; - plot(Nerrs log); + plot(Nerrs_log); #{ figure(2) @@ -884,8 +503,9 @@ function acquisition_histograms(fine_en = 0) endfunction - -% choose simulation to run here ------------------------------------------------------- +% --------------------------------------------------------- +% choose simulation to run here +% --------------------------------------------------------- format; more off; diff --git a/codec2-dev/octave/ofdm_lib.m b/codec2-dev/octave/ofdm_lib.m new file mode 100644 index 00000000..1004ddbb --- /dev/null +++ b/codec2-dev/octave/ofdm_lib.m @@ -0,0 +1,373 @@ +% ofdm_lib.m +% David Rowe Mar 2017 +% + +% Library of functions that implement a BPSK/QPSK OFDM modem. Rate Fs +% verison of ofdm_rs.m with OFDM based up and down conversion, and all +% those nasty real-world details like fine freq, timing. + +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 + + +%------------------------------------------------------------- +% ofdm_init +%------------------------------------------------------------- + +#{ + 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); + + % carrier tables for up and down conversion + + w = (0:Nc+1)*2*pi*Rs/states.Fs; + W = zeros(Nc+2,states.M); + for c=1:Nc+2 + W(c,:) = exp(j*w(c)*(0:states.M-1)); + end + states.w = w; + states.W = W; + + % fine timing search +/- window_width/2 from current timing instant + + states.window_width = 11; + + % Receive buffer: D P DDD P DDD P DDD P D + % ^ + % also see ofdm_demod() ... + + states.Nrxbuf = 3*states.Nsamperframe+states.M+states.Ncp + 2*(states.M + states.Ncp); + states.rxbuf = zeros(1, states.Nrxbuf); + + % default settings on a bunch of options and states + + states.verbose = 0; + states.timing_en = 1; + states.foff_est_en = 1; + states.phase_est_en = 1; + + states.foff_est_gain = 0.1; + states.foff_est_hz = 0; + states.sample_point = states.timing_est = 1; + states.nin = states.Nsamperframe; + + % generate OFDM pilot symbol, used for timing and freq offset est + + rate_fs_pilot_samples = states.pilots * W/states.M; + states.rate_fs_pilot_samples = [rate_fs_pilot_samples(states.M-states.Ncp+1:states.M) rate_fs_pilot_samples]; + +endfunction + + +% -------------------------------------- +% ofdm_mod - modulates one frame of bits +% -------------------------------------- + +function tx = ofdm_mod(states, tx_bits) + ofdm_load_const; + + assert(length(tx_bits) == Nbitsperframe); + + % map to symbols in linear array + + if bps == 1 + tx_sym_lin = 2*tx_bits - 1; + end + if bps == 2 + for s=1:Nbitsperframe/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; + aframe = zeros(Ns,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]; + + % OFDM upconvert symbol by symbol so we can add CP + + tx = []; + for r=1:Ns + asymbol = tx_sym(r,:) * W/M; + asymbol_cp = [asymbol(M-Ncp+1:M) asymbol]; + tx = [tx asymbol_cp]; + end +endfunction + + +% ------------------------------------------ +% ofdm_demod - Demodulates one frame of bits +% ------------------------------------------ + +#{ + + For phase estimation we need to maintain buffer of 3 frames plus + one pilot, so we have 4 pilots total. '^' is the start of current + frame that we are demodulating. + + P DDD P DDD P DDD P + ^ + + Then add one symbol either side to account for movement in + sampling instant due to sample clock differences: + + D P DDD P DDD P DDD P D + ^ +#} + +function [rx_bits states aphase_est_pilot_log rx_np] = ofdm_demod(states, rxbuf_in) + ofdm_load_const; + + % extra states that are st up at run time rather than init time + + timing_est = states.timing_est; + timing_en = states.timing_en; + foff_est_hz = states.foff_est_hz; + foff_est_gain = states.foff_est_gain; + foff_est_en = states.foff_est_en; + sample_point = states.sample_point; + rate_fs_pilot_samples = states.rate_fs_pilot_samples; + verbose = states.verbose; + phase_est_en = states.phase_est_en; + + % insert latest input samples into rxbuf + + rxbuf(1:Nrxbuf-states.nin) = rxbuf(states.nin+1:Nrxbuf); + rxbuf(Nrxbuf-states.nin+1:Nrxbuf) = rxbuf_in; + + woff_est = 2*pi*foff_est_hz/Fs; + + % update timing estimate -------------------------------------------------- + + delta_t = sample_point = 0; + if timing_en + % update timing at start of every frame + + st = M+Ncp + Nsamperframe + 1 - floor(window_width/2) + (timing_est-1); + en = st + Nsamperframe-1 + M+Ncp + window_width-1; + + ft_est = coarse_sync(states, rxbuf(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); + end + + % Black magic to keep sample_point inside cyclic prefix. Or something like that. + + 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 + + % down convert at current timing instant---------------------------------- + + % todo: this cld be more efficent, as pilot r becomes r-Ns on next frame + + rx_sym = zeros(1+Ns+1+1, Nc+2); + + % previous pilot + + st = M+Ncp + Nsamperframe + (-Ns)*(M+Ncp) + 1 + sample_point; en = st + M - 1; + + for c=1:Nc+2 + acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); + rx_sym(1,c) = sum(acarrier); + end + + % pilot - this frame - pilot + + for rr=1:Ns+1 + st = M+Ncp + Nsamperframe + (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1; + for c=1:Nc+2 + acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); + rx_sym(rr+1,c) = sum(acarrier); + end + end + + % next pilot + + st = M+Ncp + Nsamperframe + (2*Ns)*(M+Ncp) + 1 + sample_point; en = st + M - 1; + for c=1:Nc+2 + acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:)); + rx_sym(Ns+3,c) = sum(acarrier); + end + + % est freq err based on all carriers ------------------------------------ + + if foff_est_en + freq_err_rect = sum(rx_sym(2,:))' * sum(rx_sym(2+Ns,:)); + freq_err_hz = angle(freq_err_rect)*Rs/(2*pi*Ns); + foff_est_hz += foff_est_gain*freq_err_hz; + end + + % OK - now estimate and correct phase ---------------------------------- + + aphase_est_pilot = 10*ones(1,Nc+2); + 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(2,cr)*pilots(cr)') + sum(rx_sym(2+Ns,cr)*pilots(cr)'); + + % use next step of pilots in past and future + + aphase_est_pilot_rect += sum(rx_sym(1,cr)*pilots(cr)'); + aphase_est_pilot_rect += sum(rx_sym(2+Ns+1,cr)*pilots(cr)'); + + aphase_est_pilot(c) = angle(aphase_est_pilot_rect); + end + + % correct phase offset using phase estimate, and demodulate + % bits, separate loop as it runs across cols (carriers) to get + % frame bit ordering correct + + aphase_est_pilot_log = []; + rx_bits = []; rx_np = []; + for rr=1:Ns-1 + for c=2:Nc+1 + if phase_est_en + rx_corr = rx_sym(rr+2,c) * exp(-j*aphase_est_pilot(c)); + else + rx_corr = rx_sym(rr+2,c); + end + rx_np = [rx_np rx_corr]; + if bps == 1 + abit = real(rx_corr) > 0; + end + if bps == 2 + abit = qpsk_demod(rx_corr); + end + rx_bits = [rx_bits abit]; + end % c=2:Nc+1 + aphase_est_pilot_log = [aphase_est_pilot_log; aphase_est_pilot]; + end + + % Adjust nin to take care of sample clock offset + + nin = Nsamperframe; + if timing_en + thresh = (M+Ncp)/8; + tshift = (M+Ncp)/4; + if timing_est > thresh + nin = Nsamperframe+tshift; + timing_est -= tshift; + sample_point -= tshift; + end + if timing_est < -thresh + nin = Nsamperframe-tshift; + timing_est += tshift; + sample_point += tshift; + end + end + + states.rxbuf = rxbuf; + states.nin = nin; + states.timing_est = timing_est; + states.sample_point = sample_point; + states.delta_t = delta_t; + states.foff_est_hz = foff_est_hz; +endfunction + + -- 2.25.1