diff options
Diffstat (limited to 'octave/ofdm_lib.m')
| -rw-r--r-- | octave/ofdm_lib.m | 1284 |
1 files changed, 1284 insertions, 0 deletions
diff --git a/octave/ofdm_lib.m b/octave/ofdm_lib.m new file mode 100644 index 0000000..e4ae35a --- /dev/null +++ b/octave/ofdm_lib.m @@ -0,0 +1,1284 @@ +% ofdm_lib.m +% David Rowe Mar 2017 + +#{ + Library of functions that implement a PSK OFDM modem. +#} + +1; +qam16; +esno_est; +ofdm_mode; +ofdm_state; +ofdm_helper; + +%------------------------------------------------------------- +% ofdm_init +%------------------------------------------------------------- + +#{ + + Modem frame has a pilot every Ns symbols. There are Ns-1 data + symbols between every pilot. + + e.g. for Ns=4, Nc=6: + + |-Nc-| Time + DDDDDD | + PPPPPPPP --- | + DDDDDD | | + DDDDDD Ns | + DDDDDD | | + PPPPPPPP --- \|/ + DDDDDD | | + + Freq------------------> + + In this figure, time flows down, freq across. +#} + +function states = ofdm_init(config) + Rs = config.Rs; + Tcp = config.Tcp; + Ns = config.Ns; + Nc = config.Nc; + bps = config.bps; + Np = config.Np; + Ntxtbits = config.Ntxtbits; + Nuwbits = config.Nuwbits; + ftwindow_width = config.ftwindow_width; + timing_mx_thresh = config.timing_mx_thresh; + tx_uw = config.tx_uw; + bad_uw_errors = config.bad_uw_errors; + amp_scale = config.amp_scale; + amp_est_mode = config.amp_est_mode; + EsNo_est_all_symbols = config.EsNo_est_all_symbols; + EsNodB = config.EsNodB; + state_machine = config.state_machine; + edge_pilots = config.edge_pilots; + clip_gain1 = config.clip_gain1; + clip_gain2 = config.clip_gain2; + foff_limiter = config.foff_limiter; + txbpf_width_Hz = config.txbpf_width_Hz; + data_mode = config.data_mode; + + states.Fs = 8000; + states.bps = bps; + states.Rs = Rs; + states.Tcp = Tcp; + states.Ns = Ns; % one pilot every Ns symbols, e.g. Ns=4, ...PDDDPDDDP... + states.Nc = Nc; % Number of carriers + states.M = states.Fs/Rs; % oversampling rate + states.Ncp = Tcp*states.Fs; + states.Nbitsperframe = (Ns-1)*Nc*bps; % total bits in all data symbols in modem frame + states.Nsampersymbol = states.M+states.Ncp; % number of samples in a single symbol + states.Nsamperframe = Ns*states.Nsampersymbol; % number of samples in a modem frame + states.qam16 = [ + 1 + j, 1 + j*3, 3 + j, 3 + j*3; + 1 - j, 1 - j*3, 3 - j, 3 - j*3; + -1 + j, -1 + j*3, -3 + j, -3 + j*3; + -1 - j, -1 - j*3, -3 - j, -3 - j*3]/3; + rms = sqrt(states.qam16(:)'*states.qam16(:)/16);% set average Es to 1 + states.qam16 /= rms; + states.qam16 *= exp(-j*pi/4); % same rotation as QPSK constellation + states.Np = Np; % number of modem frames per packet. In some modes we want + % the total packet of data to span multiple modem frames, e.g. HF data + % and/or when the FEC codeword is larger than the one + % modem frame. In other modes (e.g. 700D/2020) Np=1, ie the modem frame + % is the same length as the packet/FEC codeword. + states.Nbitsperpacket = Np*states.Nbitsperframe; + states.Tpacket = Np*Ns*(Tcp+1/Rs); % time for one packet in ms + + states.Ntxtbits = Ntxtbits; % reserved bits/frame for auxiliary text information. Uncoded/unprotected so may + % be of limited use going forward, consider setting to 0 + states.Nuwbits = Nuwbits; + + % some basic sanity checks + assert(floor(states.M) == states.M); + + % UW symbol placement. + % Note we need to fill each UW symbols with bits. The LDPC decoder + % works on symbols so we can't break up any symbols into UW/FEC + % encoded bits. + + states.uw_ind = states.uw_ind_sym = []; + + % lets see if all UW syms will fit in frame + Nuwsyms = states.Nuwbits/bps; + Ndatasymsperframe = (Ns-1)*Nc; + states.spread_uw = 0; + if states.spread_uw + uw_step = 1.8*floor(states.Nbitsperpacket/states.Nuwbits); + else + uw_step = Nc+1; % default step for UW sym placement + end + last_sym = floor(Nuwsyms*uw_step/bps+1); + if last_sym > states.Np*Ndatasymsperframe + uw_step = Nc-1; % try a different step + end + last_sym = floor(Nuwsyms*uw_step/bps+1); + assert(last_sym <= states.Np*Ndatasymsperframe); % we still can't fit them all + + % Place UW symbols in frame + for i=1:Nuwsyms + ind_sym = floor(i*uw_step/bps+1); + % printf("%d sym: %d\n",i, ind_sym); + states.uw_ind_sym = [states.uw_ind_sym ind_sym]; % symbol index + for b=bps-1:-1:0 + states.uw_ind = [states.uw_ind bps*ind_sym-b]; % bit index + end + end + + % how many of the first few frames have UW symbols in them + Nsymsperframe = states.Nbitsperframe/states.bps; + states.Nuwframes = ceil(states.uw_ind_sym(end)/Nsymsperframe); + + states.tx_uw = tx_uw; + assert(length(states.tx_uw) == states.Nuwbits); + tx_uw_syms = []; + for b=1:bps:states.Nuwbits + if bps == 2 tx_uw_syms = [tx_uw_syms qpsk_mod(states.tx_uw(b:b+1))]; end + if bps == 4 tx_uw_syms = [tx_uw_syms qam16_mod(states.qam16, states.tx_uw(b:b+bps-1))]; end + end + states.tx_uw_syms = tx_uw_syms; + + % if the UW has this many errors it is "bad", the binomal cdf can be used to + % set this with the ofdm_determine_bad_uw_errors() function below + % + % Nuw=12; plot(0:Nuw, binocdf(0:Nuw,Nuw,0.05)); hold on; plot(binocdf(0:Nuw,Nuw,0.5)); hold off; + states.bad_uw_errors = bad_uw_errors; + + states.ofdm_peak = 16384; + % use this to scale tx output to 16 bit short to a peak value of 16384. Adjusted by experiment + states.amp_scale = amp_scale; + % when using the clipping, this is the manual gain value. Adjusted by experiment, trade off between + % increased average power and BER + states.clip_gain1 = clip_gain1; + states.clip_gain2 = clip_gain2; + states.txbpf_width_Hz = txbpf_width_Hz; + + % this is used to scale inputs to LDPC decoder to make it amplitude indep + states.mean_amp = 0; + + % use a fixed EsNo for LDPC decoder, this seems to work OK and avoid another estimator + states.EsNodB = EsNodB; + + % generate same BPSK pilots each time + rand('seed',1); + states.pilots = 1 - 2*(rand(1,Nc+2) > 0.5); + %printf("number of pilots total: %d\n", length(states.pilots)); + + % If set, place pilots at carrier 1 and Nc+2 to support low bandwidth phase est over grid + % of 12 pilot_samples. Used for 700D and 2020 + states.edge_pilots = edge_pilots; + if states.edge_pilots == 0 + states.pilots(1) = 0; + states.pilots(Nc+2) = 0; + end + + % carrier tables for up and down conversion + states.fcentre = fcentre = 1500; + alower = fcentre - Rs * (Nc/2); % approx frequency of lowest carrier + Nlower = round(alower / Rs) - 1; % round this to nearest integer multiple from 0Hz to keep DFT happy + %printf(" fcentre: %f alower: %f alower/Rs: %f Nlower: %d\n", fcentre, alower, alower/Rs, Nlower); + w = (Nlower:Nlower+Nc+1)*2*pi/(states.Fs/Rs); + 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, + % set this to roughly twice the maximum delay spread + states.ftwindow_width = ftwindow_width; + + % magic number we adjust by experiment (see ofdm_dev.m acquisition tests, blog post on 700D sync) + states.timing_mx_thresh = timing_mx_thresh; + + % Receive buffer: rxbufst + D P DDD P DDD P DDD P D + % ^ + % nominal start of current modem frame + + if length(data_mode) + Nrxbufhistory = (states.Np+2)*states.Nsamperframe; % extra storage at start of rxbuf to allow us to step back in time + else + Nrxbufhistory = 0; + end + states.rxbufst = Nrxbufhistory; % start of rxbuf window used for demod of current rx frame + states.Nrxbufhistory = Nrxbufhistory; + + % D P DDD P DDD P DDD P D + states.Nrxbufmin = states.Nsampersymbol + 3*states.Nsamperframe + states.Nsampersymbol + states.Nsampersymbol; + states.Nrxbuf = Nrxbufhistory + states.Nrxbufmin; + 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.phase_est_bandwidth = "high"; + states.dpsk = 0; + states.amp_est_mode = amp_est_mode; + + states.foff_est_gain = 0.1; + states.foff_limiter = foff_limiter; + states.foff_est_hz = 0; + states.sample_point = states.timing_est = 1; + states.nin = states.Nsamperframe; + states.timing_valid = 0; + states.timing_mx = 0; + states.coarse_foff_est_hz = 0; + + states.foff_metric = 0; + + % generate OFDM pilot symbol, used for timing and freq offset est + + rate_fs_pilot_samples = states.pilots * W/states.M; + + % During tuning it was found that not including the cyc prefix in + % rate_fs_pilot_samples produced better fest results + + %states.rate_fs_pilot_samples = [rate_fs_pilot_samples(states.M-states.Ncp+1:states.M) rate_fs_pilot_samples]; + states.rate_fs_pilot_samples = [zeros(1,states.Ncp) rate_fs_pilot_samples]; + + % pre-compute a constant used to detect valid modem frames + + Npsam = length(states.rate_fs_pilot_samples); + states.timing_norm = Npsam*(states.rate_fs_pilot_samples * states.rate_fs_pilot_samples'); + % printf("timing_norm: %f\n", states.timing_norm) + + % sync state machine + + states.sync_state = states.last_sync_state = 'search'; + states.uw_errors = 0; + states.sync_counter = 0; + states.frame_count = 0; % number of frames we have been in sync + states.sync_start = 0; + states.sync_end = 0; + states.modem_frame = 0; % keep track of how many frames received in packet + states.state_machine = state_machine; % mode specific state machine + states.packetsperburst = 0; % for OFDM data modes, how many packets before we reset state machine + states.postambledetectoren = strcmp(data_mode,"burst"); + states.npre = states.npost = 0; % counters for logging + + % LDPC code is optionally enabled + + states.rate = 1.0; + states.ldpc_en = 0; + + % init some output states for logging + + states.rx_sym = zeros(1+Ns+1+1, Nc+2); + + % Es/No (SNR) est states + + states.EsNo_est_all_symbols = EsNo_est_all_symbols; + states.clock_offset_est = 0; + + % pre-amble for data modes + states.data_mode = data_mode; + if length(states.data_mode) + states.tx_preamble = ofdm_generate_preamble(states, 2); + states.tx_postamble = ofdm_generate_preamble(states, 3); + end + + % automated tests + test_qam16_mod_demod(states.qam16); + test_assemble_disassemble(states); +endfunction + +% 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 out = freq_shift(in, foff, Fs) + foff_rect = exp(j*2*pi*foff/Fs); + foff_phase_rect = exp(j*0); + + for r=1:length(in) + foff_phase_rect *= foff_rect; + out(r) = in(r)*foff_phase_rect; + end +endfunction + + +% ----------------------------------------------------------------- +% ofdm_mod - modulates a complete packet (one or more modem frames) +% ---------------------------------------------------------------- + +function tx = ofdm_mod(states, tx_bits) + ofdm_load_const; + assert(length(tx_bits) == Nbitsperpacket); + + % map to symbols in linear array + + if bps == 1 + tx_sym_lin = 2*tx_bits - 1; + end + if bps == 2 + for s=1:Nbitsperpacket/bps + tx_sym_lin(s) = qpsk_mod(tx_bits(2*(s-1)+1:2*s)); + end + end + if bps == 4 + for s=1:Nbitsperpacket/bps + tx_sym_lin(s) = qam16_mod(states.qam16,tx_bits(4*(s-1)+1:4*s)); + end + end + + tx = ofdm_txframe(states, tx_sym_lin); +endfunction + + +% ---------------------------------------------- +% ofdm_txframe - modulates one packet of symbols +% ---------------------------------------------- + +function tx = ofdm_txframe(states, tx_sym_lin) + ofdm_load_const; + assert(length(tx_sym_lin) == Nbitsperpacket/bps); + + % place data symbols in multi-carrier frame with pilots and boundary carriers + + s = 1; tx_frame = zeros(Np*Ns,Nc+2); + for r=1:Np*Ns + if mod(r-1,Ns) == 0 + % row of pilots + tx_frame(r,:) = pilots; + else + % row of data symbols + arowofsymbols = tx_sym_lin(s:s+Nc-1); + tx_frame(r,2:Nc+1) = arowofsymbols; + s += Nc; + if states.dpsk + tx_frame(r,2:Nc+1) = tx_frame(r,2:Nc+1) .* tx_frame(r-1,2:Nc+1); + end + end + end + % make sure we use all the symbols + assert((s-1) == length(tx_sym_lin)); + + % OFDM upconvert symbol by symbol so we can add CP + + tx = []; + for r=1:Ns*Np + asymbol = tx_frame(r,:) * W/M; + asymbol_cp = [asymbol(M-Ncp+1:M) asymbol]; + tx = [tx asymbol_cp]; + end +endfunction + + +% ----------------------------------------------------------- +% est_timing +% ----------------------------------------------------------- + +#{ + Correlates known samples (for example pilots or a preamble) with a window of received + samples to determine the most likely timing offset. Optionally combines + known samples from two frames (e.g. pilots at start of this and next frame) + so we need at least Nsamperframe+M+Ncp samples in rx. + + Can be used for acquisition (coarse timing), and fine timing. Tends + to break down when freq offset approaches +/- symbol rate (e.g +/- + 25 Hz for 700D). +#} + +function [t_est timing_valid timing_mx av_level] = est_timing(states, rx, known_samples, step, dual=1) + ofdm_load_const; + Npsam = length(known_samples); + + Ncorr = length(rx) - (Nsamperframe+Npsam); + corr = zeros(1,Ncorr); + %printf("Npsam: %d M+Ncp: %d Ncorr: %d Nsamperframe: %d step: %d\n", Npsam, M+Ncp, Ncorr, Nsamperframe, step); + + % normalise correlation so we can compare to a threshold across varying input levels + + av_level = 2*sqrt(states.timing_norm*(rx*rx')/length(rx)) + 1E-12; + + % correlate with pilots at start and (optionally) end of frame to determine timing offset + + for i=1:step:Ncorr + rx1 = rx(i:i+Npsam-1); + corr_st = rx1 * known_samples'; + corr_en = 0; + if dual + % for the streaming voice modes we also correlate with pilot samples at start of next frame + rx2 = rx(i+Nsamperframe:i+Nsamperframe+Npsam-1); + corr_en = rx2 * known_samples'; + end + corr(i) = (abs(corr_st) + abs(corr_en))/av_level; + end + + [timing_mx t_est] = max(abs(corr)); + % only declare timing valid if there are enough samples in rxbuf to demodulate a frame + timing_valid = (abs(rx(t_est)) > 0) && (timing_mx > timing_mx_thresh); + + if verbose > 1 + printf(" av_level: %5.4f mx: %4.3f timing_est: %4d timing_valid: %d\n", av_level, timing_mx, t_est, timing_valid); + end + if verbose > 2 + figure(10); clf; + subplot(211); plot(rx) + subplot(212); plot(corr) + figure(11); clf; plot(real(known_samples)); + end + +endfunction + + +% ----------------------------------------------------------- +% est_freq_offset_known_corr +% ----------------------------------------------------------- + +#{ + Determines frequency offset at current timing estimate, used for + coarse freq offset estimation during streaming mode acquisition. +#} + +function foff_est = est_freq_offset_known_corr(states, rx, known_samples, t_est, dual=1) + ofdm_load_const; + Npsam = length(known_samples); + + % extract pilot samples from either end of frame + rx1 = rx(t_est:t_est+Npsam-1); rx2 = rx(t_est+Nsamperframe:t_est+Nsamperframe+Npsam-1); + + % "mix" these down (correlate) with 0 Hz offset pilot samples + corr_st = rx1 .* conj(known_samples); + if dual + corr_en = rx2 .* conj(known_samples); + end + + % sample sum of DFT magnitude of correlated signals at each freq offset and look for peak + st = -20; en = 20; foff_est = 0; Cabs_max = 0; + + for f=st:en + w = 2*pi*f/Fs; + C_st = corr_st * exp(j*w*(0:Npsam-1))'; + C_en = 0; + if dual + C_en = corr_en * exp(j*w*(0:Npsam-1))'; + end + Cabs = abs(C_st) + abs(C_en); + %printf("f: %4.1f Cabs: %f Cmax: %f\n", f, Cabs, Cabs_max); + if Cabs > Cabs_max + Cabs_max = Cabs; + foff_est = f; + end + end + + if states.verbose > 1 + printf(" foff_est: %f\n", foff_est); + end + +endfunction + + +% Joint estimation used for data mode burst acquistion + +function [t_est foff_est timing_mx] = est_timing_and_freq(states, rx, known_samples, tstep, fmin, fmax, fstep) + ofdm_load_const; + Npsam = length(known_samples); + + Ncorr = length(rx) - Npsam + 1; + corr = zeros(1,Ncorr); + + % set up matrix of freq shifted known samples for correlation with received signal. Each row + % is the known samples shifted by a different freq offset + + M = []; + for afcoarse=fmin:fstep:fmax + w = 2*pi*afcoarse/Fs; + wvec = exp(j*w*(0:Npsam-1)); + M = [M; known_samples .* wvec]; + end + + % At each timing position, correlate with known samples at all possible freq offsets. Result + % is a column vector for each timing offset. Each matrix cell is a freq,timing coordinate + + corr = []; + for t=1:tstep:Ncorr + rx1 = rx(t:t+Npsam-1); + col = M * rx1'; + corr = [corr, col]; + end + + % best timing offset is the col with the global max of the corr matrix + max_col = max(abs(corr)); + [mx mx_col] = max(max_col); + t_est = (mx_col-1)*tstep; + + % obtain normalised real number for timing mx + mag1 = known_samples*known_samples'; + mag2 = rx(t_est+1:t_est+Npsam)*rx(t_est+1:t_est+Npsam)'; + timing_mx = mx*mx'/(mag1*mag2+1E-12); + + % determine frequency offset for row where max is located + [tmp freq_row] = max(corr(:,mx_col)); + foff_est = fmin + fstep*(freq_row-1); + + if verbose > 1 + printf(" t_est: %d timing:mx: %f foff_est: %f\n", t_est, timing_mx, foff_est); + end + if verbose > 2 + figure(10); clf; + subplot(211); plot(rx) + subplot(212); plot(corr) + figure(11); clf; plot(real(known_samples)); + end + +endfunction + + +% streaming mode acquistion, used mainly for voice modes + +function [timing_valid states] = ofdm_sync_search_stream(states) + ofdm_load_const; + + st = rxbufst + M+Ncp + Nsamperframe + 1; en = st + 2*Nsamperframe + M+Ncp - 1; + + % Attempt coarse timing estimate (i.e. detect start of frame) at a range of frequency offsets + + timing_mx = 0; fcoarse = 0; timing_valid = 0; ct_est = 1; + for afcoarse=-40:40:40 + % vector of local oscillator samples to shift input vector + % these could be computed on the fly to save memory, or pre-computed in flash at tables as they are static + + if afcoarse != 0 + w = 2*pi*afcoarse/Fs; + wvec = exp(-j*w*(0:2*Nsamperframe+M+Ncp-1)); + + % choose best timing offset metric at this freq offset + [act_est atiming_valid atiming_mx] = est_timing(states, wvec .* states.rxbuf(st:en), states.rate_fs_pilot_samples, 2); + else + % exp(-j*0) is just 1 when afcoarse is 0 + [act_est atiming_valid atiming_mx] = est_timing(states, states.rxbuf(st:en), states.rate_fs_pilot_samples, 2); + end + + %printf("afcoarse: %f atiming_mx: %f\n", afcoarse, atiming_mx); + + if atiming_mx > timing_mx + ct_est = act_est; + timing_valid = atiming_valid; + timing_mx = atiming_mx; + fcoarse = afcoarse; + end + end + + % refine freq est within -/+ 20 Hz window + + if fcoarse != 0 + w = 2*pi*fcoarse/Fs; + wvec = exp(-j*w*(0:2*Nsamperframe+M+Ncp-1)); + foff_est = est_freq_offset_known_corr(states, wvec .* states.rxbuf(st:en), states.rate_fs_pilot_samples, ct_est); + foff_est += fcoarse; + else + % exp(-j*0) is just 1 when fcoarse is 0 + foff_est = est_freq_offset_known_corr(states, states.rxbuf(st:en), states.rate_fs_pilot_samples, ct_est); + end + + if verbose + printf(" ct_est: %4d mx: %3.2f coarse_foff: %5.1f timing_valid: %d", ct_est, timing_mx, foff_est, timing_valid); + end + + if timing_valid + states.nin = ct_est - 1; + else + states.nin = Nsamperframe; + end + + states.timing_valid = timing_valid; + states.timing_mx = timing_mx; + states.coarse_foff_est_hz = foff_est; + states.sample_point = states.timing_est = 1; +endfunction + + +% two stage acquisition detector for burst mode + +function results = burst_acquisition_detector(states, rx, n, known_sequence) + ofdm_load_const; + + % initial search over coarse grid + tstep = 4; fstep = 5; + [ct_est foff_est timing_mx] = est_timing_and_freq(states, rx(n:n+2*Nsamperframe-1), known_sequence, + tstep, fmin = -50, fmax = 50, fstep); + % refine estimate over finer grid + fmin = foff_est - ceil(fstep/2); fmax = foff_est + ceil(fstep/2); + fine_st = max(1, n + ct_est - tstep/2); fine_en = fine_st + Nsamperframe + tstep - 1; + [ct_est foff_est timing_mx] = est_timing_and_freq(states, rx(fine_st:fine_en), known_sequence, 1, fmin, fmax, 1); + % refer ct_est to nominal start of frame rx_buf(n) + ct_est += fine_st - n; + results.ct_est = ct_est; results.foff_est = foff_est; results.timing_mx = timing_mx; +end + + +% Burst mode acquisition ------------------------------------------ + +function [timing_valid states] = ofdm_sync_search_burst(states) + ofdm_load_const; + + pre_post = ""; + st = rxbufst + M+Ncp + Nsamperframe + 1; en = st + 2*Nsamperframe - 1; + pre = burst_acquisition_detector(states, states.rxbuf, st, states.tx_preamble); + if states.postambledetectoren + post = burst_acquisition_detector(states, states.rxbuf, st, states.tx_postamble); + end + + if isfield(states,"postambletest") pre.timing_mx = 0; end % force ignore preamble to test postamble + + if (states.postambledetectoren == 0) || (pre.timing_mx > post.timing_mx) + timing_mx = pre.timing_mx; ct_est = pre.ct_est; foff_est = pre.foff_est; + pre_post = "pre"; + else + timing_mx = post.timing_mx; ct_est = post.ct_est; foff_est = post.foff_est; + pre_post = "post"; + end + timing_valid = timing_mx > timing_mx_thresh; + + if timing_valid + % potential candidate found .... + + % calculate number of samples we need on next buffer to get into sync + if strcmp(pre_post, "post") + states.nin = 0; + % printf("\n rxbufst: %d ", states.rxbufst); + states.rxbufst -= states.Np*states.Nsamperframe; % backup to first modem frame in packet + states.rxbufst += ct_est - 1; + states.npost++; + % printf("%d\n", states.rxbufst); + else + % ct_est is start of preamble, so advance past that to start of first modem frame + states.nin = Nsamperframe + ct_est - 1; + states.npre++; + end + else + states.nin = Nsamperframe; + end + + states.ct_est = ct_est; + states.timing_valid = timing_valid; + states.timing_mx = timing_mx; + states.sample_point = states.timing_est = 1; + states.foff_est_hz = foff_est; + + if verbose + printf(" ct_est: %4d nin: %4d mx: %3.2f foff_est: %5.1f timing_valid: %d %4s", + ct_est, states.nin, timing_mx, foff_est, timing_valid, pre_post); + end +endfunction + + +% ---------------------------------------------------------------------------------- +% ofdm_sync_search - attempts to find coarse sync parameters for modem initial sync +% ---------------------------------------------------------------------------------- + +function [timing_valid states] = ofdm_sync_search(states, rxbuf_in) + ofdm_load_const; + + % update rxbuf so it is primed for when we have to call ofdm_demod() + + states.rxbuf(1:Nrxbuf-states.nin) = states.rxbuf(states.nin+1:Nrxbuf); + states.rxbuf(Nrxbuf-states.nin+1:Nrxbuf) = rxbuf_in; + + if strcmp(states.data_mode, "burst") + [timing_valid states] = ofdm_sync_search_burst(states); + else + [timing_valid states] = ofdm_sync_search_stream(states); + 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 + ^ + + Returns: + rx_bits - (hard decoded/raw/uncoded) demodulated data bits from packet + aphase_est - phase est for each data symbol + rx_np - output data symbols after phase correction + rx_amp - amplitude estimates for each symbol +#} + +function [states rx_bits achannel_est_rect_log rx_np rx_amp] = ofdm_demod(states, rxbuf_in) + ofdm_load_const; + + % insert latest input samples into rxbuf + + rxbuf(1:Nrxbuf-states.nin) = rxbuf(states.nin+1:Nrxbuf); + rxbuf(Nrxbuf-states.nin+1:Nrxbuf) = rxbuf_in; + + % get latest freq offset estimate + + woff_est = 2*pi*foff_est_hz/Fs; + + % update timing estimate -------------------------------------------------- + + delta_t = coarse_foff_est_hz = timing_valid = timing_mx = 0; + if timing_en + % update timing at start of every frame + + % search for timing in a window centered on timing_est, the window will typically be around 2Ncp wide as we could + % get a shift of +Ncp or -Ncp if we swing from one delay extreme to another + st = rxbufst + M+Ncp + Nsamperframe + 1 - floor(ftwindow_width/2) + (timing_est-1); + en = st + Nsamperframe-1 + M+Ncp + ftwindow_width-1; + + [ft_est timing_valid timing_mx] = est_timing(states, rxbuf(st:en) .* exp(-j*woff_est*(st:en)), rate_fs_pilot_samples, 1); + % printf(" timing_est: %d ft_est: %d timing_valid: %d timing_mx: %d\n", timing_est, ft_est, timing_valid, timing_mx); + + % if we are in a deep fade timing_valid will not be asserted as ft_est will be garbage, so we don't + % adjust timing est, just freewheel for now + if timing_valid + + % adjust timing_est based on ft_est + timing_est = timing_est + ft_est - ceil(ftwindow_width/2); + + % Track the ideal sampling point, which is Ncp for a multipath signal whose delay varies between 0 and Ncp. The + % timing est will be bouncing back and forth due to multipath so we may need to use the upper or lower limit of + % the timing est to track the ideal sample_point. A good way to explore this algorithm is to disable the feedback + % loop for nin adjustment below, and look at the plots from ofdm_rx with +ve and -ve sample clock offsets + % (sox can be used to resample). The "4" constants are small guard bands so we don't stumble outside of the CP + % due to noise. + + delta_t = ft_est - ceil(ftwindow_width/2); % just used for plotting + sample_point = max(timing_est+4, sample_point); % we are at max timing est, so sample point just above + sample_point = min(timing_est+Ncp-4, sample_point); % we are at min timing_est, so sample point Ncp above + end + + if verbose > 1 + printf(" ft_est: %2d mx: %3.2f coarse_foff: %4.1f foff: %4.1f\n", ft_est, timing_mx, coarse_foff_est_hz, foff_est_hz); + end + + end + + % down convert at current timing instant---------------------------------- + + rx_sym = zeros(1+Ns+1+1, Nc+2); + + % previous pilot + + st = rxbufst + 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 = rxbufst + 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 = rxbufst + 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,:)); + + % prevent instability in atan(im/re) when real part near 0 + + freq_err_rect += 1E-6; + + %printf("freq_err_rect: %f %f angle: %f\n", real(freq_err_rect), imag(freq_err_rect), angle(freq_err_rect)); + freq_err_hz = angle(freq_err_rect)*Rs/(2*pi*Ns); + if states.foff_limiter + freq_err_hz = max(freq_err_hz,-1); + freq_err_hz = min(freq_err_hz, 1); + end + foff_est_hz = foff_est_hz + foff_est_gain*freq_err_hz; + end + + % OK - now channel for each carrier and correct phase ---------------------------------- + + achannel_est_rect = zeros(1,Nc+2); + aamp_est_pilot = zeros(1,Nc+2); + for c=2:Nc+1 + + % estimate channel for this carrier using an average of 12 pilots + % in a rect 2D window centred on this carrier + + % PPP <-- frame-1 + % --- + % PPP <-- you are here + % DDD + % DDD + % PPP <-- frame+1 + % --- + % PPP <-- frame+2 + + if isfield(states, "phase_est_bandwidth") + phase_est_bandwidth = states.phase_est_bandwidth; + else + phase_est_bandwidth = "low"; + end + + if strcmp(phase_est_bandwidth, "high") + % Only use pilots at start and end of this frame to track quickly changes in phase + % present. Useful for initial sync where freq offset est may be a bit off, and + % for high Doppler channels. As less pilots are averaged, low SNR performance + % will be poorer. + achannel_est_rect(c) = rx_sym(2,c)*pilots(c)'; % frame + achannel_est_rect(c) += rx_sym(2+Ns,c)*pilots(c)'; % frame+1 + aamp_est_pilot(c) = abs(rx_sym(2,c)) + abs(rx_sym(2+Ns,c)); + elseif strcmp(phase_est_bandwidth, "low") + % Average over a bunch of pilots in adjacent carriers, and past and future frames, good + % low SNR performance, but will fall over with high Doppler or freq offset. + cr = c-1:c+1; + achannel_est_rect(c) = rx_sym(2,cr)*pilots(cr)'; % frame + achannel_est_rect(c) += rx_sym(2+Ns,cr)*pilots(cr)'; % frame+1 + aamp_est_pilot(c) = sum(abs(rx_sym(2,cr))); + aamp_est_pilot(c) += sum(abs(rx_sym(2+Ns,cr))); + + % use next step of pilots in past and future + + achannel_est_rect(c) += rx_sym(1,cr)*pilots(cr)'; % frame-1 + achannel_est_rect(c) += rx_sym(2+Ns+1,cr)*pilots(cr)'; % frame+2 + aamp_est_pilot(c) += sum(abs(rx_sym(1,cr))); + aamp_est_pilot(c) += sum(abs(rx_sym(2+Ns+1,cr))); + end + end + + % pilots are estimated over 12 pilot symbols, so find average + + if strcmp(phase_est_bandwidth, "high") + achannel_est_rect /= 2; + aamp_est_pilot /= 2; + elseif strcmp(phase_est_bandwidth, "low") + achannel_est_rect /= 12; + aamp_est_pilot /= 12; + end + + aphase_est_pilot = angle(achannel_est_rect); + if states.amp_est_mode == 0 + % legacy 700D/2020 ampl estimator for compatibility with current C code + aamp_est_pilot = abs(achannel_est_rect); + end + achannel_est_rect = aamp_est_pilot.*exp(j*aphase_est_pilot); + + % correct phase offset using phase estimate, and demodulate + % bits, separate loop as it runs across cols (carriers) to get + % frame bit ordering correct + + rx_bits = []; rx_np = []; rx_amp = []; achannel_est_rect_log = []; + for rr=1:Ns-1 + for c=2:Nc+1 + if phase_est_en + if states.dpsk + rx_corr = rx_sym(rr+2,c) * rx_sym(rr+1,c)'; + else + rx_corr = rx_sym(rr+2,c) * exp(-j*aphase_est_pilot(c)); + end + else + rx_corr = rx_sym(rr+2,c); + end + + rx_np = [rx_np rx_corr]; + rx_amp = [rx_amp aamp_est_pilot(c)]; + + % hard decision demod + if bps == 1 abit = real(rx_corr) > 0; end + if bps == 2 abit = qpsk_demod(rx_corr); end + if bps == 4 abit = qam16_demod(states.qam16, rx_corr, max(1E-12,aamp_est_pilot(c))); end + rx_bits = [rx_bits abit]; + end % c=2:Nc+1 + achannel_est_rect_log = [achannel_est_rect_log; achannel_est_rect(2:Nc+1)]; + end + + % Adjust nin to take care of sample clock offset. When debugong or exploring how timing loop works + % it's a good idea to comment out this code to "open the loop". + + nin = Nsamperframe; + if timing_en && timing_valid + states.clock_offset_est = 0.9*states.clock_offset_est + 0.1*abs(states.timing_est - timing_est)/Nsamperframe; + 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 + + % use internal rxbuf samples if they are available + rxbufst_next = rxbufst + nin; + %printf("\nrxbufst: %d rxbufst_next: %d nin: %d Nrxbufmin: %d rqd: %d Nrxbuf: %d\n", + % rxbufst, rxbufst_next, nin, Nrxbufmin, rxbufst_next + Nrxbufmin, Nrxbuf); + if rxbufst_next + Nrxbufmin <= Nrxbuf + % printf("Can maybe use rxbufst!\n"); + rxbufst = rxbufst_next; + nin = 0; + end + + % maintain mean amp estimate for LDPC decoder + states.mean_amp = 0.9*states.mean_amp + 0.1*mean(rx_amp); + + states.rx_sym = rx_sym; + states.rxbuf = rxbuf; + states.nin = nin; + states.rxbufst = rxbufst; + states.timing_valid = timing_valid; + states.timing_mx = timing_mx; + states.timing_est = timing_est; + states.sample_point = sample_point; + states.delta_t = delta_t; + states.foff_est_hz = foff_est_hz; + states.coarse_foff_est_hz = coarse_foff_est_hz; % just used for tofdm +endfunction + + +function SNR3kdB = snr_from_esno(states, EsNodB) + ofdm_load_const; + + % We integrate over M samples to get the received symbols. Additional signal power + % is used for the cyclic prefix samples. + cyclic_power = 10*log10((Ncp+M)/M); + % Es is the energy for each symbol. To get signal power lets + % multiply by symbols/second, and calculate noise power in 3000 Hz. + SNR3kdB = EsNodB + 10*log10(Nc*Rs/3000) + cyclic_power; +endfunction + +% ---------------------------------------------------------------------------------- +% assemble_modem_packet - assemble modem packet from UW, payload, and txt bits +% ---------------------------------------------------------------------------------- + +function modem_frame = assemble_modem_packet(states, payload_bits, txt_bits) + ofdm_load_const; + + # Due to the operation of the FEC encoder or interleaver, Tx data + # usually comes in "packet size" chunks, so assembly operates on an + # entire packet (multiple modem frames if Np>1) + + p = 1; u = 1; + modem_frame = zeros(1,Nbitsperpacket); + + for b=1:Nbitsperpacket-Ntxtbits; + if (u <= Nuwbits) && (b == uw_ind(u)) + modem_frame(b) = tx_uw(u++); + else + modem_frame(b) = payload_bits(p++); + end + end + t = 1; + for b=Nbitsperpacket-Ntxtbits+1:Nbitsperpacket + modem_frame(b) = txt_bits(t++); + end + assert(u == (Nuwbits+1)); + assert(p = (length(payload_bits)+1)); +endfunction + + +% ---------------------------------------------------------------------------------- +% assemble_modem_packet_symbols - assemble modem packet from UW, payload, and txt bits +% ---------------------------------------------------------------------------------- + +function modem_frame = assemble_modem_packet_symbols(states, payload_syms, txt_syms) + ofdm_load_const; + + Nsymsperpacket = Nbitsperpacket/bps; + Nuwsyms = Nuwbits/bps; + Ntxtsyms = Ntxtbits/bps; + modem_frame = zeros(1,Nsymsperpacket); + p = 1; u = 1; + + for s=1:Nsymsperpacket-Ntxtsyms; + if (u <= Nuwsyms) && (s == uw_ind_sym(u)) + modem_frame(s) = states.tx_uw_syms(u++); + else + modem_frame(s) = payload_syms(p++); + end + end + t = 1; + for s=Nsymsperpacket-Ntxtsyms+1:Nsymsperpacket + modem_frame(s) = txt_syms(t++); + end + assert(u == (Nuwsyms+1)); + assert(p = (length(payload_syms)+1)); +endfunction + + +% ------------------------------------------------------------------------------------------------ +% extract_uw - extract just the UW from the first few frames of a packet, to check UW +% during acquisition +% ------------------------------------------------------------------------------------------------- + +function rx_uw = extract_uw(states, rx_syms, rx_amps) + ofdm_load_const; + + Nsymsperframe = Nbitsperframe/bps; + assert(length(rx_syms) == Nuwframes*Nsymsperframe); + Nuwsyms = Nuwbits/bps; + rx_uw_syms = zeros(1,Nuwsyms); + rx_uw_amps = zeros(1,Nuwsyms); + u = 1; + + for s=1:Nuwframes*Nsymsperframe + if (u <= Nuwsyms) && (s == uw_ind_sym(u)) + rx_uw_syms(u) = rx_syms(s); + rx_uw_amps(u) = rx_amps(s); + u++; + end + end + assert(u == (Nuwsyms+1)); + + % now demodulate UW bits + rx_uw = zeros(1,Nuwbits); + + for s=1:Nuwsyms + if bps == 2 + rx_uw(bps*(s-1)+1:bps*s) = qpsk_demod(rx_uw_syms(s)); + elseif bps == 4 + rx_uw(bps*(s-1)+1:bps*s) = qam16_demod(states.qam16,rx_uw_syms(s), max(1E-12,rx_amps(s))); + end + end +endfunction + + +% ------------------------------------------------------------------------------------------------ +% disassemble_modem_packet - extract UW, txt bits, and payload symbols from a packet of symbols +% ------------------------------------------------------------------------------------------------- + +function [rx_uw payload_syms payload_amps txt_bits] = disassemble_modem_packet(states, modem_frame_syms, modem_frame_amps) + ofdm_load_const; + + Nsymsperpacket = Nbitsperpacket/bps; + Nuwsyms = Nuwbits/bps; + Ntxtsyms = Ntxtbits/bps; + payload_syms = zeros(1,Nsymsperpacket-Nuwsyms-Ntxtsyms); + payload_amps = zeros(1,Nsymsperpacket-Nuwsyms-Ntxtsyms); + rx_uw_syms = zeros(1,Nuwsyms); + rx_uw_amps = zeros(1,Nuwsyms); + txt_syms = zeros(1,Ntxtsyms); + p = 1; u = 1; + + for s=1:Nsymsperpacket-Ntxtsyms; + if (u <= Nuwsyms) && (s == uw_ind_sym(u)) + rx_uw_syms(u) = modem_frame_syms(s); + rx_uw_amps(u) = modem_frame_amps(s); + u++; + else + payload_syms(p) = modem_frame_syms(s); + payload_amps(p++) = modem_frame_amps(s); + end + end + t = 1; + for s=Nsymsperpacket-Ntxtsyms+1:Nsymsperpacket + txt_syms(t++) = modem_frame_syms(s); + end + assert(u == (Nuwsyms+1)); + assert(p = (Nsymsperpacket+1)); + + % now demodulate UW and txt bits + + rx_uw = zeros(1,Nuwbits); + txt_bits = zeros(1,Ntxtbits); + + for s=1:Nuwsyms + if bps == 2 + rx_uw(bps*(s-1)+1:bps*s) = qpsk_demod(rx_uw_syms(s)); + elseif bps == 4 + rx_uw(bps*(s-1)+1:bps*s) = qam16_demod(states.qam16,rx_uw_syms(s),rx_uw_amps(s)); + end + end + for s=1:Ntxtsyms + txt_bits(2*s-1:2*s) = qpsk_demod(txt_syms(s)); + end + +endfunction + + +%----------------------------------------------------------------------- +% ofdm_rand - a psuedo-random number generator that we can implement +% in C with identical results to Octave. Returns an unsigned +% int between 0 and 32767 +%----------------------------------------------------------------------- + +function r = ofdm_rand(n, seed=1) + r = zeros(1,n); + for i=1:n + seed = mod(1103515245 * seed + 12345, 32768); + r(i) = seed; + end +endfunction + + +% build a single modem frame preamble vector for reliable single frame acquisition +% on data modes +function tx_preamble = ofdm_generate_preamble(states, seed=2) + tmp_states = states; + % tweak local copy of states so we can generate a 1 modem-frame packet + tmp_states.Np = 1; tmp_states.Nbitsperpacket = tmp_states.Nbitsperframe; + preamble_bits = ofdm_rand(tmp_states.Nbitsperframe, seed) > 16384; + tx_preamble = ofdm_mod(tmp_states, preamble_bits); +endfunction + + +% ------------------------------------------------------------------------------ +% Handle FEC encoding/decoding +% ------------------------------------------------------------------------------ + +function [frame_bits bits_per_frame] = fec_encode(states, code_param, mode, payload_bits) + ofdm_load_const; + if code_param.data_bits_per_frame != code_param.ldpc_data_bits_per_frame + % optionally lower the code rate by "one stuffing" - setting Nunused data bits to 1 + Nunused = code_param.ldpc_data_bits_per_frame - code_param.data_bits_per_frame; + frame_bits = LdpcEncode([payload_bits ones(1,Nunused)], code_param.H_rows, code_param.P_matrix); + % remove unused data bits from codeword, as they are known to the receiver and don't need to be transmitted + frame_bits = [ frame_bits(1:code_param.data_bits_per_frame) frame_bits(code_param.ldpc_data_bits_per_frame+1:end) ]; + else + frame_bits = LdpcEncode(payload_bits, code_param.H_rows, code_param.P_matrix); + end + bits_per_frame = length(frame_bits); + +endfunction + +function [rx_bits paritychecks] = fec_decode(states, code_param, ... + payload_syms_de, payload_amps_de, ... + mean_amp, EsNo) + ofdm_load_const; + % note ldpc_dec() handles optional lower code rate zero-stuffing + [rx_codeword paritychecks] = ldpc_dec(code_param, mx_iter=100, demod=0, dec=0, ... + payload_syms_de/mean_amp, EsNo, + payload_amps_de/mean_amp); + rx_bits = rx_codeword(1:code_param.data_bits_per_frame); +endfunction + + +function [tx nclipped] = ofdm_clip(states, tx, threshold_level, plot_en=0) + ofdm_load_const; + tx_ = tx; + ind = find(abs(tx) > threshold_level); + nclipped = length(ind); + tx(ind) = threshold_level*exp(j*angle(tx(ind))); + if plot_en + figure(2); clf; plot(abs(tx_(1:5*M))); hold on; plot(abs(tx(1:5*M))); hold off; + endif +end + +% two stage Hilbert clipper to improve PAPR +function tx = ofdm_hilbert_clipper(states, tx, tx_clip_en) + tx *= states.amp_scale; + + % optional compressor to improve PAPR + + nclipped = 0; + if tx_clip_en + if states.verbose + printf("%f %f\n", states.clip_gain1, states.clip_gain2); + end + [tx nclipped] = ofdm_clip(states, tx*states.clip_gain1, states.ofdm_peak); + + cutoff_norm = states.txbpf_width_Hz/states.Fs; + w_centre = mean(states.w); centre_norm = w_centre/(2*pi); + tx = ofdm_complex_bandpass_filter(cutoff_norm, centre_norm,100,tx); + + % filter messes up peak levels use this to get us back to approx 16384 + tx *= states.clip_gain2; + end + + % Hilbert Clipper 2 - remove any really low probability outliers after clipping/filtering + % even on vanilla Tx + [tx tmp] = ofdm_clip(states, tx, states.ofdm_peak); + + % note this is PAPR of complex signal, PAPR of real signal will be 3dB-ish larger + peak = max(abs(tx)); RMS = sqrt(mean(abs(tx).^2)); + cpapr = 10*log10((peak.^2)/(RMS.^2)); + + if states.verbose + printf("Peak: %4.2f RMS: %5.2f CPAPR: %4.2f clipped: %5.2f%%\n", + peak, RMS, cpapr, nclipped*100/length(tx)); + end +endfunction + + +% Complex bandpass filter built from low pass prototype as per src/filter.c, +% cutoff_freq and center_freq are normalised such that cutoff_freq = 0.5 is Fs/2 +function out = ofdm_complex_bandpass_filter(cutoff_freq,center_freq,n_coeffs,in) + lowpass_coeff = fir1(n_coeffs-1, cutoff_freq); + k = (0:n_coeffs-1); + bandpass_coeff = lowpass_coeff .* exp(j*2*pi*center_freq*k); + out = filter(bandpass_coeff,1,in); +endfunction + + +% Complex bandpass filter for Rx - just used on the very low SNR modes to help +% with acquisition +function [rx delay_samples] = ofdm_rx_filter(states, mode, rx) + delay_samples = 0; + if strcmp(mode,"datac4") || strcmp(mode,"datac13") || strcmp(mode,"datac14") + w_centre = mean(states.w); centre_norm = w_centre/(2*pi); + n_coeffs = 100; + cutoff_Hz = 400; cutoff_norm = cutoff_Hz/states.Fs; + rx = ofdm_complex_bandpass_filter(cutoff_norm,centre_norm,n_coeffs,rx); + delay_samples = n_coeffs/2; + end +endfunction + + +% returns an unpacked CRC16 (array of 16 bits) calculated from an array of unpacked bits +function unpacked_crc16 = crc16_unpacked(unpacked_bits) + % pack into bytes + mod(length(unpacked_bits),8); + assert(mod(length(unpacked_bits),8) == 0); + nbytes = length(unpacked_bits)/8; + mask = 2 .^ (7:-1:0); + for i=1:nbytes + st = (i-1)*8 + 1; en = st+7; + bytes(i) = sum(mask .* unpacked_bits(st:en)); + end + crc16_hex = crc16(bytes); + crc16_dec = [hex2dec(crc16_hex(1:2)) hex2dec(crc16_hex(3:4)) ]; + unpacked_crc16 = []; + for b=1:length(crc16_dec) + unpacked_crc16 = [unpacked_crc16 bitand(crc16_dec(b), mask) > 0]; + end +endfunction |