From: drowe67 Date: Sun, 18 Mar 2018 20:14:20 +0000 (+0000) Subject: renamed, as about to build simpler ofdtm_tx and ofdem_rx to test OFDM mdoem only X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=908329556443bd77d4e158194ae3140dc3d2cf5e;p=freetel-svn-tracking.git renamed, as about to build simpler ofdtm_tx and ofdem_rx to test OFDM mdoem only git-svn-id: https://svn.code.sf.net/p/freetel/code@3416 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/ofdm_ldpc_rx.m b/codec2-dev/octave/ofdm_ldpc_rx.m new file mode 100644 index 00000000..98476c5a --- /dev/null +++ b/codec2-dev/octave/ofdm_ldpc_rx.m @@ -0,0 +1,295 @@ +% ofdm_ldpc_rx.m +% David Rowe April 2017 +% +% OFDM file based rx, with LDPC and interleaver + +#{ + TODO: + [ ] proper EsNo estimation + [ ] some sort of real time GUI display to watch signal evolving + [ ] est SNR or Eb/No of recieved signal + [ ] way to fall out of sync +#} + +function ofdm_rx(filename, interleave_frames = 1, error_pattern_filename) + ofdm_lib; + ldpc; + gp_interleaver; + more off; + + % init modem + + Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8; + states = ofdm_init(bps, Rs, Tcp, Ns, Nc); + ofdm_load_const; + states.verbose = 1; + + % Set up LDPC code + + mod_order = 4; bps = 2; modulation = 'QPSK'; mapping = 'gray'; + demod_type = 0; decoder_type = 0; max_iterations = 100; + + EsNo = 10; % TODO: fixme + + init_cml('/home/david/Desktop/cml/'); + load HRA_112_112.txt + [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping); + assert(Nbitsperframe == code_param.code_bits_per_frame); + + % load real samples from file + + Ascale= 2E5*1.1491; + frx=fopen(filename,"rb"); rx = 2*fread(frx, Inf, "short")/4E5; fclose(frx); + Nsam = length(rx); Nframes = floor(Nsam/Nsamperframe); + prx = 1; + + % buffers for interleaved frames + + Nsymbolsperframe = code_param.code_bits_per_frame/bps + Nsymbolsperinterleavedframe = interleave_frames*Nsymbolsperframe + rx_np = zeros(1, Nsymbolsperinterleavedframe); + rx_amp = zeros(1, Nsymbolsperinterleavedframe); + + % OK generate tx frame for BER calcs + + rand('seed', 100); + tx_bits = []; tx_codewords = []; + for f=1:interleave_frames + atx_bits = round(rand(1,code_param.data_bits_per_frame)); + tx_bits = [tx_bits atx_bits]; + acodeword = LdpcEncode(atx_bits, code_param.H_rows, code_param.P_matrix); + tx_codewords = [tx_codewords acodeword]; + end + + % used for rx frame sync on interleaved symbols - we demod the + % entire interleaved frame to raw bits + + tx_symbols = []; + for s=1:Nsymbolsperinterleavedframe + tx_symbols = [tx_symbols qpsk_mod( tx_codewords(2*(s-1)+1:2*s) )]; + end + + tx_symbols = gp_interleave(tx_symbols); + + tx_bits_raw = []; + for s=1:Nsymbolsperinterleavedframe + tx_bits_raw = [tx_bits_raw qpsk_demod(tx_symbols(s))]; + end + + % init logs and BER stats + + rx_bits = []; rx_np_log = []; timing_est_log = []; delta_t_log = []; foff_est_hz_log = []; + phase_est_pilot_log = []; + Terrs = Tbits = Terrs_coded = Tbits_coded = Tpackets = Tpacketerrs = 0; + Nbitspervocframe = 28; + Nerrs_coded_log = Nerrs_log = []; + error_positions = []; + + % 'prime' rx buf to get correct coarse timing (for now) + + prx = 1; + nin = Nsamperframe+2*(M+Ncp); + states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx(prx:nin); + prx += nin; + + state = 'searching'; frame_count = 0; + + % main loop ---------------------------------------------------------------- + + for f=1:Nframes + + % insert samples at end of buffer, set to zero if no samples + % available to disable phase estimation on future pilots on last + % frame of simulation + + lnew = min(Nsam-prx,states.nin); + rxbuf_in = zeros(1,states.nin); + + if lnew + rxbuf_in(1:lnew) = rx(prx:prx+lnew-1); + end + prx += states.nin; + + [rx_bits_raw states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in); + frame_count++; + + % update sliding windows of rx-ed symbols and symbol amplitudes + + rx_np(1:Nsymbolsperinterleavedframe-Nsymbolsperframe) = rx_np(Nsymbolsperframe+1:Nsymbolsperinterleavedframe); + rx_np(Nsymbolsperinterleavedframe-Nsymbolsperframe+1:Nsymbolsperinterleavedframe) = arx_np; + rx_amp(1:Nsymbolsperinterleavedframe-Nsymbolsperframe) = rx_amp(Nsymbolsperframe+1:Nsymbolsperinterleavedframe); + rx_amp(Nsymbolsperinterleavedframe-Nsymbolsperframe+1:Nsymbolsperinterleavedframe) = arx_amp; + + printf("f: %d state: %s frame_count: %d\n", f, state, frame_count); + + % If looking for sync: check raw BER on frame just received + % against all possible positions in the interleaver frame. + + % iterate state machine ------------------------------------ + + next_state = state; + if strcmp(state,'searching') + + % If looking for sync: check raw BER on frame just received + % against all possible positions in the interleaver frame. + + for ff=1:interleave_frames + st = (ff-1)*Nbitsperframe+1; en = st + Nbitsperframe - 1; + errors = xor(tx_bits_raw(st:en), rx_bits_raw); + Nerrs = sum(errors); + printf(" ff: %d Nerrs: %d\n", ff, Nerrs); + if Nerrs/Nbitsperframe < 0.1 + next_state = 'synced'; + % make sure we get an interleave frame with correct freq offset + % note this introduces a lot of delay, a better idea would be to + % run demod again from interleave_frames back with now-known freq offset + frame_count = ff - interleave_frames; + end + end + end + + if strcmp(state,'synced') + if Nerrs/Nbitsperframe > 0.2 + %next_state = 'searching'; + end + end + + state = next_state; + + if strcmp(state,'searching') + + % still searching? Attempt coarse timing estimate (i.e. detect start of frame) + + st = M+Ncp + Nsamperframe + 1; en = st + 2*Nsamperframe; + [ct_est foff_est] = coarse_sync(states, states.rxbuf(st:en), states.rate_fs_pilot_samples); + if states.verbose + printf(" Nerrs: %d ct_est: %4d foff_est: %3.1f\n", Nerrs, ct_est, foff_est); + end + + % calculate number of samples we need on next buffer to get into sync + + states.nin = Nsamperframe + ct_est - 1; + + % reset modem states + + states.sample_point = states.timing_est = 1; + states.foff_est_hz = foff_est; + end + + if strcmp(state,'synced') + + % we are in sync so log states + + rx_np_log = [rx_np_log arx_np]; + timing_est_log = [timing_est_log states.timing_est]; + delta_t_log = [delta_t_log states.delta_t]; + foff_est_hz_log = [foff_est_hz_log states.foff_est_hz]; + phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log]; + end + + if strcmp(state,'synced') && (frame_count == interleave_frames) + + % de-interleave QPSK symbols + + arx_np = gp_deinterleave(rx_np); + arx_amp = gp_deinterleave(rx_amp); + + % measure uncoded bit errors per modem frame + + rx_bits_raw = []; + for s=1:Nsymbolsperinterleavedframe + rx_bits_raw = [rx_bits_raw qpsk_demod(rx_np(s))]; + end + for ff=1:interleave_frames + st = (ff-1)*Nbitsperframe+1; en = st+Nbitsperframe-1; + errors = xor(tx_bits_raw(st:en), rx_bits_raw(st:en)); + Nerrs = sum(errors); + if Nerrs/Nbitsperframe < 0.2 + Terrs += Nerrs; + Nerrs_log = [Nerrs_log Nerrs]; + Tbits += Nbitsperframe; + end + end + + % LDPC decode + + rx_bits = []; + for ff=1:interleave_frames + st = (ff-1)*Nbitsperframe/bps+1; en = st + Nbitsperframe/bps - 1; + rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, arx_np(st:en), min(EsNo,30), arx_amp(st:en)); + rx_bits = [rx_bits rx_codeword(1:code_param.data_bits_per_frame)]; + end + + errors_coded = xor(tx_bits, rx_bits); + Nerrs_coded = sum(errors_coded); + if Nerrs/Nbitsperframe < 0.2 + Terrs_coded += Nerrs_coded; + Tbits_coded += code_param.data_bits_per_frame*interleave_frames; + error_positions = [error_positions errors_coded]; + + % measure packet errors based on Codec 2 packet size + + Nvocframes = floor(code_param.data_bits_per_frame*interleave_frames/Nbitspervocframe); + for fv=1:Nvocframes + st = (fv-1)*Nbitspervocframe + 1; + en = st + Nbitspervocframe - 1; + Nvocpacketerrs = sum(xor(tx_bits(st:en), rx_bits(st:en))); + if Nvocpacketerrs + Tpacketerrs++; + end + Tpackets++; + Nerrs_coded_log = [Nerrs_coded_log Nvocpacketerrs]; + end + end + printf(" Nerrs_coded: %d\n", Nerrs_coded); + + frame_count = 0; + end + end + + printf("Coded BER: %5.4f Tbits: %5d Terrs: %5d PER: %5.4f Tpacketerrs: %5d Tpackets: %5d\n", + Terrs_coded/Tbits_coded, Tbits_coded, Terrs_coded, Tpacketerrs/Tpackets, Tpacketerrs, Tpackets); + printf("Raw BER..: %5.4f Tbits: %5d Terrs: %5d\n", Terrs/Tbits, Tbits, Terrs); + + figure(1); clf; + plot(rx_np_log,'+'); + mx = max(abs(rx_np_log)) + axis([-mx mx -mx mx]); + title('Scatter'); + + figure(2); clf; + plot(phase_est_pilot_log(:,2:Nc+1),'g+', 'markersize', 5); + title('Phase est'); + axis([1 length(phase_est_pilot_log) -pi pi]); + + figure(3); clf; + subplot(211) + stem(delta_t_log) + title('delta t'); + subplot(212) + plot(timing_est_log); + title('timing est'); + + figure(4); clf; + plot(foff_est_hz_log) + mx = max(abs(foff_est_hz_log)); + axis([1 max(Nframes,2) -mx mx]); + title('Fine Freq'); + ylabel('Hz') + + figure(5); clf; + subplot(211) + stem(Nerrs_log); + title('Uncoded errrors/modem frame') + axis([1 length(Nerrs_log) 0 Nbitsperframe*rate/2]); + subplot(212) + stem(Nerrs_coded_log); + title('Coded errrors/vocoder frame') + axis([1 length(Nerrs_coded_log) 0 Nbitspervocframe/2]); + + if nargin == 3 + fep = fopen(error_pattern_filename, "wb"); + fwrite(fep, error_positions, "short"); + fclose(fep); + end +endfunction diff --git a/codec2-dev/octave/ofdm_ldpc_tx.m b/codec2-dev/octave/ofdm_ldpc_tx.m new file mode 100644 index 00000000..05bf92bf --- /dev/null +++ b/codec2-dev/octave/ofdm_ldpc_tx.m @@ -0,0 +1,189 @@ +% ofdm_ldpc_tx.m +% David Rowe April 2017 +% +% File based ofdm tx with LDPC encoding and interleaver. Generates a +% file of ofdm samples, including optional channel simulation. + +#{ + Examples: + + i) 4 frame interleaver, 10 seconds, AWGN channel at Eb/No=3dB + + octave:4> ofdm_tx('awgn_ebno_3dB_700d.raw',4, 10,3); + + ii) 4 frame interleaver, 10 seconds, HF channel at Eb/No=6dB + + ofdm_tx('hf_ebno_6dB_700d.raw', 4, 10, 6, 'hf'); +#} + + +#{ + TODO: + [ ] measure and report raw and coded BER + [ ] maybe 10s worth of frames, sync up to any one automatically + + or start with maybe 10 frames + + measure BER match on each one + [ ] model clipping/PA compression + [ ] sample clock offsets + [ ] compare with same SNR from pathsim + [ ] How to pack arbitrary frames into ofdm frame and codec 2 frames + + integer number of ofdm frames? + + how to sync at other end + +#} + +function ofdm_tx(filename, Nsec, interleave_frames = 1, EbNodB=100, channel='awgn', freq_offset_Hz=0) + ofdm_lib; + ldpc; + gp_interleaver; + + % init modem + + Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8; + states = ofdm_init(bps, Rs, Tcp, Ns, Nc); + ofdm_load_const; + + % Set up LDPC code + + mod_order = 4; bps = 2; modulation = 'QPSK'; mapping = 'gray'; + + init_cml('/home/david/Desktop/cml/'); + load HRA_112_112.txt + [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping); + assert(Nbitsperframe == code_param.code_bits_per_frame); + + % Generate fixed test frame of tx bits and run OFDM modulator + + Nrows = Nsec*Rs; + Nframes = floor((Nrows-1)/Ns); + + % Adjust Nframes so we have an integer number of interleaver frames + % in simulation + + Nframes = interleave_frames*round(Nframes/interleave_frames); + + % OK generate an interleaver frame of codewords from random data bits + + rand('seed', 100); + tx_bits = tx_symbols = []; + for f=1:interleave_frames + atx_bits = round(rand(1,code_param.data_bits_per_frame)); + tx_bits = [tx_bits atx_bits]; + codeword = LdpcEncode(atx_bits, code_param.H_rows, code_param.P_matrix); + for b=1:2:Nbitsperframe + tx_symbols = [tx_symbols qpsk_mod(codeword(b:b+1))]; + end + end + + tx_symbols = gp_interleave(tx_symbols); + + atx = []; + for f=1:interleave_frames + st = (f-1)*Nbitsperframe/bps+1; en = st + Nbitsperframe/bps-1; + atx = [atx ofdm_txframe(states, tx_symbols(st:en))]; + end + + tx = []; + for f=1:Nframes/interleave_frames + tx = [tx atx]; + end + + % not very pretty way to process analog signals with exactly the same channel + % todo: work out a cleaner way + + analog_hack = 0; rx_filter = 0; + if analog_hack || rx_filter + + % simulated SSB tx filter + + [b, a] = cheby1(4, 3, [600, 2600]/(Fs/2)); + h = freqz(b,a,(600:2600)/(Fs/(2*pi))); + filt_gain = (2600-600)/sum(abs(h) .^ 2); % ensures power after filter == before filter + end + + if analog_hack + % load analog signal and convert to complex + + s = load_raw('../raw/ve9qrp_10s.raw')'; + tx = hilbert(s); + + % ssb tx filter + + tx = filter(b,a,sqrt(filt_gain)*tx); + + % normalise power to same as ofdm tx + + nom_tx_pwr = 2/(Ns*(M*M)) + Nc/(M*M); + tx_pwr = var(tx); + tx *= sqrt(nom_tx_pwr/tx_pwr); + end + + Nsam = length(tx); + + % channel simulation + + EsNo = rate * bps * (10 .^ (EbNodB/10)); + variance = 1/(M*EsNo/2); + woffset = 2*pi*freq_offset_Hz/Fs; + + SNRdB = EbNodB + 10*log10(700/3000); + printf("EbNo: %3.1f dB SNR(3k) est: %3.1f dB foff: %3.1fHz ", EbNodB, SNRdB, freq_offset_Hz); + + % set up HF model --------------------------------------------------------------- + + if strcmp(channel, 'hf') + randn('seed',1); + + % some typical values, or replace with user supplied + + dopplerSpreadHz = 1; path_delay_ms = 1; + + 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, (Nsec*(M+Ncp)/M)*Fs*1.1); + spread2 = doppler_spread(dopplerSpreadHz, Fs, (Nsec*(M+Ncp)/M)*Fs*1.1); + + % sometimes doppler_spread() doesn't return exactly the number of samples we need + + assert(length(spread1) >= Nsam, "not enough doppler spreading samples"); + assert(length(spread2) >= Nsam, "not enough doppler spreading samples"); + end + + rx = tx; + + if strcmp(channel, 'hf') + rx = tx(1:Nsam) .* spread1(1:Nsam); + rx += [zeros(1,path_delay_samples) tx(1:Nsam-path_delay_samples)] .* spread2(1:Nsam); + + % normalise rx power to same as tx + + nom_rx_pwr = 2/(Ns*(M*M)) + Nc/(M*M); + rx_pwr = var(rx); + rx *= sqrt(nom_rx_pwr/rx_pwr); + end + + rx = rx .* exp(j*woffset*(1:Nsam)); + + % note variance/2 as we are using real() operator, mumble, + % reflection of -ve freq to +ve, mumble, hand wave + + noise = sqrt(variance/2)*0.5*randn(1,Nsam); + rx = real(rx) + noise; + printf("measured SNR: %3.2f dB\n", 10*log10(var(real(tx))/var(noise))+10*log10(4000) - 10*log10(3000)); + + if rx_filter + % ssb rx filter + rx = filter(b,a,sqrt(filt_gain)*rx); + end + + % adjusted by experiment to match rms power of early test signals + + Ascale = 2E5*1.1491; + + frx=fopen(filename,"wb"); fwrite(frx, Ascale*rx, "short"); fclose(frx); +endfunction diff --git a/codec2-dev/octave/ofdm_rx.m b/codec2-dev/octave/ofdm_rx.m deleted file mode 100644 index d108181b..00000000 --- a/codec2-dev/octave/ofdm_rx.m +++ /dev/null @@ -1,295 +0,0 @@ -% ofdm_rx.m -% David Rowe April 2017 -% -% OFDM file based rx. - -#{ - TODO: - [ ] proper EsNo estimation - [ ] some sort of real time GUI display to watch signal evolving - [ ] est SNR or Eb/No of recieved signal - [ ] way to fall out of sync -#} - -function ofdm_rx(filename, interleave_frames = 1, error_pattern_filename) - ofdm_lib; - ldpc; - gp_interleaver; - more off; - - % init modem - - Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8; - states = ofdm_init(bps, Rs, Tcp, Ns, Nc); - ofdm_load_const; - states.verbose = 1; - - % Set up LDPC code - - mod_order = 4; bps = 2; modulation = 'QPSK'; mapping = 'gray'; - demod_type = 0; decoder_type = 0; max_iterations = 100; - - EsNo = 10; % TODO: fixme - - init_cml('/home/david/Desktop/cml/'); - load HRA_112_112.txt - [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping); - assert(Nbitsperframe == code_param.code_bits_per_frame); - - % load real samples from file - - Ascale= 2E5*1.1491; - frx=fopen(filename,"rb"); rx = 2*fread(frx, Inf, "short")/4E5; fclose(frx); - Nsam = length(rx); Nframes = floor(Nsam/Nsamperframe); - prx = 1; - - % buffers for interleaved frames - - Nsymbolsperframe = code_param.code_bits_per_frame/bps - Nsymbolsperinterleavedframe = interleave_frames*Nsymbolsperframe - rx_np = zeros(1, Nsymbolsperinterleavedframe); - rx_amp = zeros(1, Nsymbolsperinterleavedframe); - - % OK generate tx frame for BER calcs - - rand('seed', 100); - tx_bits = []; tx_codewords = []; - for f=1:interleave_frames - atx_bits = round(rand(1,code_param.data_bits_per_frame)); - tx_bits = [tx_bits atx_bits]; - acodeword = LdpcEncode(atx_bits, code_param.H_rows, code_param.P_matrix); - tx_codewords = [tx_codewords acodeword]; - end - - % used for rx frame sync on interleaved symbols - we demod the - % entire interleaved frame to raw bits - - tx_symbols = []; - for s=1:Nsymbolsperinterleavedframe - tx_symbols = [tx_symbols qpsk_mod( tx_codewords(2*(s-1)+1:2*s) )]; - end - - tx_symbols = gp_interleave(tx_symbols); - - tx_bits_raw = []; - for s=1:Nsymbolsperinterleavedframe - tx_bits_raw = [tx_bits_raw qpsk_demod(tx_symbols(s))]; - end - - % init logs and BER stats - - rx_bits = []; rx_np_log = []; timing_est_log = []; delta_t_log = []; foff_est_hz_log = []; - phase_est_pilot_log = []; - Terrs = Tbits = Terrs_coded = Tbits_coded = Tpackets = Tpacketerrs = 0; - Nbitspervocframe = 28; - Nerrs_coded_log = Nerrs_log = []; - error_positions = []; - - % 'prime' rx buf to get correct coarse timing (for now) - - prx = 1; - nin = Nsamperframe+2*(M+Ncp); - states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx(prx:nin); - prx += nin; - - state = 'searching'; frame_count = 0; - - % main loop ---------------------------------------------------------------- - - for f=1:Nframes - - % insert samples at end of buffer, set to zero if no samples - % available to disable phase estimation on future pilots on last - % frame of simulation - - lnew = min(Nsam-prx,states.nin); - rxbuf_in = zeros(1,states.nin); - - if lnew - rxbuf_in(1:lnew) = rx(prx:prx+lnew-1); - end - prx += states.nin; - - [rx_bits_raw states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in); - frame_count++; - - % update sliding windows of rx-ed symbols and symbol amplitudes - - rx_np(1:Nsymbolsperinterleavedframe-Nsymbolsperframe) = rx_np(Nsymbolsperframe+1:Nsymbolsperinterleavedframe); - rx_np(Nsymbolsperinterleavedframe-Nsymbolsperframe+1:Nsymbolsperinterleavedframe) = arx_np; - rx_amp(1:Nsymbolsperinterleavedframe-Nsymbolsperframe) = rx_amp(Nsymbolsperframe+1:Nsymbolsperinterleavedframe); - rx_amp(Nsymbolsperinterleavedframe-Nsymbolsperframe+1:Nsymbolsperinterleavedframe) = arx_amp; - - printf("f: %d state: %s frame_count: %d\n", f, state, frame_count); - - % If looking for sync: check raw BER on frame just received - % against all possible positions in the interleaver frame. - - % iterate state machine ------------------------------------ - - next_state = state; - if strcmp(state,'searching') - - % If looking for sync: check raw BER on frame just received - % against all possible positions in the interleaver frame. - - for ff=1:interleave_frames - st = (ff-1)*Nbitsperframe+1; en = st + Nbitsperframe - 1; - errors = xor(tx_bits_raw(st:en), rx_bits_raw); - Nerrs = sum(errors); - printf(" ff: %d Nerrs: %d\n", ff, Nerrs); - if Nerrs/Nbitsperframe < 0.1 - next_state = 'synced'; - % make sure we get an interleave frame with correct freq offset - % note this introduces a lot of delay, a better idea would be to - % run demod again from interleave_frames back with now-known freq offset - frame_count = ff - interleave_frames; - end - end - end - - if strcmp(state,'synced') - if Nerrs/Nbitsperframe > 0.2 - %next_state = 'searching'; - end - end - - state = next_state; - - if strcmp(state,'searching') - - % still searching? Attempt coarse timing estimate (i.e. detect start of frame) - - st = M+Ncp + Nsamperframe + 1; en = st + 2*Nsamperframe; - [ct_est foff_est] = coarse_sync(states, states.rxbuf(st:en), states.rate_fs_pilot_samples); - if states.verbose - printf(" Nerrs: %d ct_est: %4d foff_est: %3.1f\n", Nerrs, ct_est, foff_est); - end - - % calculate number of samples we need on next buffer to get into sync - - states.nin = Nsamperframe + ct_est - 1; - - % reset modem states - - states.sample_point = states.timing_est = 1; - states.foff_est_hz = foff_est; - end - - if strcmp(state,'synced') - - % we are in sync so log states - - rx_np_log = [rx_np_log arx_np]; - timing_est_log = [timing_est_log states.timing_est]; - delta_t_log = [delta_t_log states.delta_t]; - foff_est_hz_log = [foff_est_hz_log states.foff_est_hz]; - phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log]; - end - - if strcmp(state,'synced') && (frame_count == interleave_frames) - - % de-interleave QPSK symbols - - arx_np = gp_deinterleave(rx_np); - arx_amp = gp_deinterleave(rx_amp); - - % measure uncoded bit errors per modem frame - - rx_bits_raw = []; - for s=1:Nsymbolsperinterleavedframe - rx_bits_raw = [rx_bits_raw qpsk_demod(rx_np(s))]; - end - for ff=1:interleave_frames - st = (ff-1)*Nbitsperframe+1; en = st+Nbitsperframe-1; - errors = xor(tx_bits_raw(st:en), rx_bits_raw(st:en)); - Nerrs = sum(errors); - if Nerrs/Nbitsperframe < 0.2 - Terrs += Nerrs; - Nerrs_log = [Nerrs_log Nerrs]; - Tbits += Nbitsperframe; - end - end - - % LDPC decode - - rx_bits = []; - for ff=1:interleave_frames - st = (ff-1)*Nbitsperframe/bps+1; en = st + Nbitsperframe/bps - 1; - rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, arx_np(st:en), min(EsNo,30), arx_amp(st:en)); - rx_bits = [rx_bits rx_codeword(1:code_param.data_bits_per_frame)]; - end - - errors_coded = xor(tx_bits, rx_bits); - Nerrs_coded = sum(errors_coded); - if Nerrs/Nbitsperframe < 0.2 - Terrs_coded += Nerrs_coded; - Tbits_coded += code_param.data_bits_per_frame*interleave_frames; - error_positions = [error_positions errors_coded]; - - % measure packet errors based on Codec 2 packet size - - Nvocframes = floor(code_param.data_bits_per_frame*interleave_frames/Nbitspervocframe); - for fv=1:Nvocframes - st = (fv-1)*Nbitspervocframe + 1; - en = st + Nbitspervocframe - 1; - Nvocpacketerrs = sum(xor(tx_bits(st:en), rx_bits(st:en))); - if Nvocpacketerrs - Tpacketerrs++; - end - Tpackets++; - Nerrs_coded_log = [Nerrs_coded_log Nvocpacketerrs]; - end - end - printf(" Nerrs_coded: %d\n", Nerrs_coded); - - frame_count = 0; - end - end - - printf("Coded BER: %5.4f Tbits: %5d Terrs: %5d PER: %5.4f Tpacketerrs: %5d Tpackets: %5d\n", - Terrs_coded/Tbits_coded, Tbits_coded, Terrs_coded, Tpacketerrs/Tpackets, Tpacketerrs, Tpackets); - printf("Raw BER..: %5.4f Tbits: %5d Terrs: %5d\n", Terrs/Tbits, Tbits, Terrs); - - figure(1); clf; - plot(rx_np_log,'+'); - mx = max(abs(rx_np_log)) - axis([-mx mx -mx mx]); - title('Scatter'); - - figure(2); clf; - plot(phase_est_pilot_log(:,2:Nc+1),'g+', 'markersize', 5); - title('Phase est'); - axis([1 length(phase_est_pilot_log) -pi pi]); - - figure(3); clf; - subplot(211) - stem(delta_t_log) - title('delta t'); - subplot(212) - plot(timing_est_log); - title('timing est'); - - figure(4); clf; - plot(foff_est_hz_log) - mx = max(abs(foff_est_hz_log)); - axis([1 max(Nframes,2) -mx mx]); - title('Fine Freq'); - ylabel('Hz') - - figure(5); clf; - subplot(211) - stem(Nerrs_log); - title('Uncoded errrors/modem frame') - axis([1 length(Nerrs_log) 0 Nbitsperframe*rate/2]); - subplot(212) - stem(Nerrs_coded_log); - title('Coded errrors/vocoder frame') - axis([1 length(Nerrs_coded_log) 0 Nbitspervocframe/2]); - - if nargin == 3 - fep = fopen(error_pattern_filename, "wb"); - fwrite(fep, error_positions, "short"); - fclose(fep); - end -endfunction diff --git a/codec2-dev/octave/ofdm_tx.m b/codec2-dev/octave/ofdm_tx.m deleted file mode 100644 index e057998f..00000000 --- a/codec2-dev/octave/ofdm_tx.m +++ /dev/null @@ -1,189 +0,0 @@ -% ofdm_tx.m -% David Rowe April 2017 -% -% File based ofdm tx. Generate a file of ofdm samples, inclduing -% optional channel simulation. - -#{ - Examples: - - i) 4 frame interleaver, 10 seconds, AWGN channel at Eb/No=3dB - - octave:4> ofdm_tx('awgn_ebno_3dB_700d.raw',4, 10,3); - - ii) 4 frame interleaver, 10 seconds, HF channel at Eb/No=6dB - - ofdm_tx('hf_ebno_6dB_700d.raw', 4, 10, 6, 'hf'); -#} - - -#{ - TODO: - [ ] measure and report raw and coded BER - [ ] maybe 10s worth of frames, sync up to any one automatically - + or start with maybe 10 frames - + measure BER match on each one - [ ] model clipping/PA compression - [ ] sample clock offsets - [ ] compare with same SNR from pathsim - [ ] How to pack arbitrary frames into ofdm frame and codec 2 frames - + integer number of ofdm frames? - + how to sync at other end - -#} - -function ofdm_tx(filename, Nsec, interleave_frames = 1, EbNodB=100, channel='awgn', freq_offset_Hz=0) - ofdm_lib; - ldpc; - gp_interleaver; - - % init modem - - Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8; - states = ofdm_init(bps, Rs, Tcp, Ns, Nc); - ofdm_load_const; - - % Set up LDPC code - - mod_order = 4; bps = 2; modulation = 'QPSK'; mapping = 'gray'; - - init_cml('/home/david/Desktop/cml/'); - load HRA_112_112.txt - [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping); - assert(Nbitsperframe == code_param.code_bits_per_frame); - - % Generate fixed test frame of tx bits and run OFDM modulator - - Nrows = Nsec*Rs; - Nframes = floor((Nrows-1)/Ns); - - % Adjust Nframes so we have an integer number of interleaver frames - % in simulation - - Nframes = interleave_frames*round(Nframes/interleave_frames); - - % OK generate an interleaver frame of codewords from random data bits - - rand('seed', 100); - tx_bits = tx_symbols = []; - for f=1:interleave_frames - atx_bits = round(rand(1,code_param.data_bits_per_frame)); - tx_bits = [tx_bits atx_bits]; - codeword = LdpcEncode(atx_bits, code_param.H_rows, code_param.P_matrix); - for b=1:2:Nbitsperframe - tx_symbols = [tx_symbols qpsk_mod(codeword(b:b+1))]; - end - end - - tx_symbols = gp_interleave(tx_symbols); - - atx = []; - for f=1:interleave_frames - st = (f-1)*Nbitsperframe/bps+1; en = st + Nbitsperframe/bps-1; - atx = [atx ofdm_txframe(states, tx_symbols(st:en))]; - end - - tx = []; - for f=1:Nframes/interleave_frames - tx = [tx atx]; - end - - % not very pretty way to process analog signals with exactly the same channel - % todo: work out a cleaner way - - analog_hack = 0; rx_filter = 0; - if analog_hack || rx_filter - - % simulated SSB tx filter - - [b, a] = cheby1(4, 3, [600, 2600]/(Fs/2)); - h = freqz(b,a,(600:2600)/(Fs/(2*pi))); - filt_gain = (2600-600)/sum(abs(h) .^ 2); % ensures power after filter == before filter - end - - if analog_hack - % load analog signal and convert to complex - - s = load_raw('../raw/ve9qrp_10s.raw')'; - tx = hilbert(s); - - % ssb tx filter - - tx = filter(b,a,sqrt(filt_gain)*tx); - - % normalise power to same as ofdm tx - - nom_tx_pwr = 2/(Ns*(M*M)) + Nc/(M*M); - tx_pwr = var(tx); - tx *= sqrt(nom_tx_pwr/tx_pwr); - end - - Nsam = length(tx); - - % channel simulation - - EsNo = rate * bps * (10 .^ (EbNodB/10)); - variance = 1/(M*EsNo/2); - woffset = 2*pi*freq_offset_Hz/Fs; - - SNRdB = EbNodB + 10*log10(700/3000); - printf("EbNo: %3.1f dB SNR(3k) est: %3.1f dB foff: %3.1fHz ", EbNodB, SNRdB, freq_offset_Hz); - - % set up HF model --------------------------------------------------------------- - - if strcmp(channel, 'hf') - randn('seed',1); - - % some typical values, or replace with user supplied - - dopplerSpreadHz = 1; path_delay_ms = 1; - - 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, (Nsec*(M+Ncp)/M)*Fs*1.1); - spread2 = doppler_spread(dopplerSpreadHz, Fs, (Nsec*(M+Ncp)/M)*Fs*1.1); - - % sometimes doppler_spread() doesn't return exactly the number of samples we need - - assert(length(spread1) >= Nsam, "not enough doppler spreading samples"); - assert(length(spread2) >= Nsam, "not enough doppler spreading samples"); - end - - rx = tx; - - if strcmp(channel, 'hf') - rx = tx(1:Nsam) .* spread1(1:Nsam); - rx += [zeros(1,path_delay_samples) tx(1:Nsam-path_delay_samples)] .* spread2(1:Nsam); - - % normalise rx power to same as tx - - nom_rx_pwr = 2/(Ns*(M*M)) + Nc/(M*M); - rx_pwr = var(rx); - rx *= sqrt(nom_rx_pwr/rx_pwr); - end - - rx = rx .* exp(j*woffset*(1:Nsam)); - - % note variance/2 as we are using real() operator, mumble, - % reflection of -ve freq to +ve, mumble, hand wave - - noise = sqrt(variance/2)*0.5*randn(1,Nsam); - rx = real(rx) + noise; - printf("measured SNR: %3.2f dB\n", 10*log10(var(real(tx))/var(noise))+10*log10(4000) - 10*log10(3000)); - - if rx_filter - % ssb rx filter - rx = filter(b,a,sqrt(filt_gain)*rx); - end - - % adjusted by experiment to match rms power of early test signals - - Ascale = 2E5*1.1491; - - frx=fopen(filename,"wb"); fwrite(frx, Ascale*rx, "short"); fclose(frx); -endfunction