From ab9d8866c41f972b578a9656b3323395d85d662d Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sat, 13 May 2017 21:59:39 +0000 Subject: [PATCH] one frame ldpc enc/dec working for ofdm_rx/ofdm_tx git-svn-id: https://svn.code.sf.net/p/freetel/code@3134 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/ldpc.m | 48 +++++++++++++++++++++++++++---- codec2-dev/octave/ofdm_dev.m | 31 ++++++++++++-------- codec2-dev/octave/ofdm_lib.m | 47 ++---------------------------- codec2-dev/octave/ofdm_rx.m | 34 +++++++++++++++++----- codec2-dev/octave/ofdm_tx.m | 55 +++++++++++++++++++++++------------- 5 files changed, 128 insertions(+), 87 deletions(-) diff --git a/codec2-dev/octave/ldpc.m b/codec2-dev/octave/ldpc.m index 1151adfd..3cbf7fe6 100644 --- a/codec2-dev/octave/ldpc.m +++ b/codec2-dev/octave/ldpc.m @@ -17,21 +17,59 @@ 1; - -function code_param = ldpc_init(rate, framesize, modulation, mod_order, mapping) +% NOTE: You will need to set the CML path in the call to init_cml() below +% for you CML install. See lpdc.m for instructions on how to install +% CML library + +function init_cml(path_to_cml) + currentdir = pwd; + + if exist(path_to_cml, 'dir') == 7 + cd(path_to_cml) + CmlStartup + cd(currentdir); + else + printf("\n---------------------------------------------------\n"); + printf("Can't start CML in path: %s\n", path_to_cml); + printf("See CML path instructions at top of this script\n"); + printf("-----------------------------------------------------\n\n"); + assert(0); + end +end + +% init using built in WiMax code + +function code_param = ldpc_init_wimax(rate, framesize, modulation, mod_order, mapping) [code_param.H_rows, code_param.H_cols, code_param.P_matrix] = InitializeWiMaxLDPC( rate, framesize, 0 ); code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix ); code_param.S_matrix = CreateConstellation( modulation, mod_order, mapping ); code_param.bits_per_symbol = log2(mod_order); code_param.code_bits_per_frame = framesize; - code_param.symbols_per_frame = framesize/2; + code_param.symbols_per_frame = framesize/code_param.bits_per_symbol; endfunction +% init using user supplied code + +function [code_param framesize rate] = ldpc_init_user(HRA, modulation, mod_order, mapping) + [Nr Nc] = size(HRA); + rate = (Nc-Nr)/Nc; + framesize = Nc; + [H_rows, H_cols] = Mat2Hrows(HRA); + code_param.H_rows = H_rows; + code_param.H_cols = H_cols; + code_param.P_matrix = []; + code_param.data_bits_per_frame = length(code_param.H_cols) - length(code_param.P_matrix); + code_param.S_matrix = CreateConstellation( modulation, mod_order, mapping ); + code_param.bits_per_symbol = log2(mod_order); + code_param.code_bits_per_frame = framesize; + code_param.symbols_per_frame = framesize/code_param.bits_per_symbol; +endfunction + function [codeword s] = ldpc_enc(data, code_param) - codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix ); - s = Modulate( codeword, code_param.S_matrix ); + codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix ); + s = Modulate( codeword, code_param.S_matrix ); endfunction diff --git a/codec2-dev/octave/ofdm_dev.m b/codec2-dev/octave/ofdm_dev.m index 2bace4f4..bad8e5f7 100644 --- a/codec2-dev/octave/ofdm_dev.m +++ b/codec2-dev/octave/ofdm_dev.m @@ -125,8 +125,7 @@ function [sim_out rx states] = run_sim(sim_in) states.code_param = code_param; elseif diversity_en - rate = 0.5; - + rate = 0.5; else rate = 1; end @@ -216,7 +215,7 @@ function [sim_out rx states] = run_sim(sim_in) tx = []; for f=1:Nframes st = (f-1)*Nbitsperframe/bps+1; en = st + Nbitsperframe/bps - 1; - tx = [tx ofdm_tx(states, tx_symbols(st:en))]; + tx = [tx ofdm_txframe(states, tx_symbols(st:en))]; end % add extra row of pilots at end, to allow one frame simulations, @@ -247,12 +246,14 @@ function [sim_out rx states] = run_sim(sim_in) end end - rx = tx; if hf_en + rx = tx; + + if hf_en 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 (e.g. 1/M) + % normalise rx power to same as tx nom_rx_pwr = 2/(Ns*(M*M)) + Nc/(M*M); rx_pwr = var(rx); @@ -426,7 +427,7 @@ function [sim_out rx states] = run_sim(sim_in) % returns results for plotting curves - if ldpc_en + if ldpc_en || diversity_en sim_out.ber(nn) = Terrs_coded/(Nbits*rate); sim_out.per(nn) = Tpacketerrs_coded/Tpackets_coded; else @@ -568,7 +569,7 @@ function run_curves sim_in.ldpc_en = 0; sim_in.hf_en = 0; - sim_in.Nsec = 10; + sim_in.Nsec = 20; sim_in.EbNodB = 0:8; awgn_EbNodB = sim_in.EbNodB; @@ -580,12 +581,13 @@ function run_curves sim_in.hf_en = 1; sim_in.ldpc_en = 0; sim_in.Nsec = 60; - sim_in.EbNodB = 4:1:14; + sim_in.EbNodB = 4:2:14; EbNoLin = 10.^(sim_in.EbNodB/10); hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); hf = run_sim(sim_in); + sim_in.diversity_en = 1; hf_diversity = run_sim(sim_in); sim_in.diversity_en = 0; sim_in.ldpc_en = 1; hf_ldpc = run_sim(sim_in); % try a few interleavers @@ -594,7 +596,7 @@ function run_curves sim_in.interleave_frames = 8; hf_ldpc_8 = run_sim(sim_in); sim_in.interleave_frames = 16; hf_ldpc_16 = run_sim(sim_in); sim_in.interleave_frames = 32; hf_ldpc_32 = run_sim(sim_in); - + % Rate Fs modem BER curves of various coding schemes figure(1); clf; @@ -603,6 +605,7 @@ function run_curves semilogy(sim_in.EbNodB, hf_theory,'b+-;HF theory;'); semilogy(awgn_EbNodB, awgn.ber,'r+-;AWGN;'); semilogy(sim_in.EbNodB, hf.ber,'r+-;HF;'); + semilogy(sim_in.EbNodB, hf_diversity.ber,'ro-;HF diversity;'); semilogy(awgn_EbNodB, awgn_ldpc.ber,'c+-;AWGN LDPC (224,112);'); semilogy(sim_in.EbNodB, hf_ldpc.ber,'c+-;HF LDPC (224,112);'); semilogy(sim_in.EbNodB, hf_ldpc_1.ber,'m+-;HF LDPC (224,112) interleave 1;'); @@ -630,11 +633,12 @@ function run_curves semilogy(sim_in.EbNodB, hf_per_theory,'b+-;HF theory;'); semilogy(awgn_EbNodB, awgn.per,'r+-;AWGN;'); semilogy(sim_in.EbNodB, hf.per,'r+-;HF sim;'); + semilogy(sim_in.EbNodB, hf_diversity.per,'ro-;HF diversity;'); semilogy(awgn_EbNodB, awgn_ldpc.per,'c+-;AWGN LDPC (224,112);'); semilogy(sim_in.EbNodB, hf_ldpc.per,'c+-;HF LDPC (224,112);'); semilogy(sim_in.EbNodB, hf_ldpc_1.per,'m+-;HF LDPC (224,112) interleave 1;'); semilogy(sim_in.EbNodB, hf_ldpc_8.per,'g+-;HF LDPC (224,112) interleave 8;'); - semilogy(sim_in.EbNodB, hf_ldpc_16.per,'k+-;HF LDPC (224,112) interleave 16;'); + semilogy(sim_in.EbNodB, hf_ldpc_16.per,'ko-;HF LDPC (224,112) interleave 16;'); semilogy(sim_in.EbNodB, hf_ldpc_32.per,'k+-;HF LDPC (224,112) interleave 32;'); hold off; axis([0 14 1E-2 1]) @@ -693,17 +697,20 @@ function run_curves snr_awgn = snr_awgn_theory + overhead_dB; snr_hf = sim_in.EbNodB + 10*log10(700/3000) + overhead_dB; - % 2/6 symbols are pilots, 1dB implementation loss + % est 700C: 2/6 symbols are pilots, 1dB implementation loss snr_awgn_700c = awgn_EbNodB + 10*log10(700/3000) + 10*log10(6/4) + 1; + snr_hf_700c = sim_in.EbNodB + 10*log10(700/3000) + 10*log10(6/4) + 1; figure(5); clf; semilogy(snr_awgn_theory, awgn_theory,'b+-;AWGN theory;'); hold on; semilogy(snr_awgn_700c, awgn_theory,'g+-;AWGN 700C;'); + semilogy(snr_hf_700c, hf_diversity.ber,'go-;AWGN 700C;'); semilogy(snr_hf_theory, hf_theory,'b+-;HF theory;'); semilogy(snr_awgn, awgn_ldpc.ber,'c+-;AWGN LDPC (224,112);'); semilogy(snr_hf, hf_ldpc.ber,'c+-;HF LDPC (224,112);'); + semilogy(snr_hf, hf_diversity.ber,'bo-;HF diversity;'); semilogy(snr_hf, hf_ldpc_16.ber,'k+-;HF LDPC (224,112) interleave 16;'); hold off; axis([-5 8 1E-3 2E-1]) @@ -712,7 +719,7 @@ function run_curves grid; grid minor on; legend('boxoff'); legend("location", "southwest"); - title('Rate Fs modem BER versus SNR inclduing pilot/CP overhead'); + title('Rate Fs modem BER versus SNR including pilot/CP overhead'); print('-deps', '-color', "ofdm_dev_ber_snr.eps") end diff --git a/codec2-dev/octave/ofdm_lib.m b/codec2-dev/octave/ofdm_lib.m index 1eb7ded2..605d6bfe 100644 --- a/codec2-dev/octave/ofdm_lib.m +++ b/codec2-dev/octave/ofdm_lib.m @@ -174,33 +174,13 @@ function states = ofdm_init(bps, Rs, Tcp, Ns, Nc) endfunction -% NOTE: You will need to set the CML path in the call to init_cml() below -% for you CML install. See lpdc.m for instructions on how to install -% CML library - -function init_cml(path_to_cml) - currentdir = pwd; - - if exist(path_to_cml, 'dir') == 7 - cd(path_to_cml) - CmlStartup - cd(currentdir); - else - printf("\n---------------------------------------------------\n"); - printf("Can't start CML in path: %s\n", path_to_cml); - printf("See CML path instructions at top of this script\n"); - printf("-----------------------------------------------------\n\n"); - assert(0); - end -end - % -------------------------------------- % 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 @@ -214,7 +194,7 @@ function tx = ofdm_mod(states, tx_bits) end end - tx = ofdm_tx(states, tx_sym_lin); + tx = ofdm_txframe(states, tx_sym_lin); endfunction % ----------------------------------------- @@ -234,9 +214,8 @@ endfunction #} -function tx = ofdm_tx(states, tx_sym_lin) +function tx = ofdm_txframe(states, tx_sym_lin) ofdm_load_const; - assert(length(tx_sym_lin) == Nbitsperframe/bps); % place symbols in multi-carrier frame with pilots and boundary carriers @@ -435,23 +414,3 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, states.foff_est_hz = foff_est_hz; endfunction - -function detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r, EsNo, fading) - if nargin == 6 - fading = ones(1, length(r)); - end - - symbol_likelihood = Demod2D( r, code_param.S_matrix, EsNo, fading); - - % initialize the extrinsic decoder input - - input_somap_c = zeros(1, code_param.code_bits_per_frame ); - bit_likelihood = Somap( symbol_likelihood, demod_type, input_somap_c ); - - input_decoder_c = bit_likelihood(1:code_param.code_bits_per_frame); - - x_hat= MpDecode( -input_decoder_c, code_param.H_rows, code_param.H_cols, ... - max_iterations, decoder_type, 1, 1); - detected_data = x_hat(max_iterations,:); -endfunction - diff --git a/codec2-dev/octave/ofdm_rx.m b/codec2-dev/octave/ofdm_rx.m index b2401b13..ba3b1425 100644 --- a/codec2-dev/octave/ofdm_rx.m +++ b/codec2-dev/octave/ofdm_rx.m @@ -5,6 +5,7 @@ #{ 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 @@ -12,16 +13,33 @@ function ofdm_rx(filename) ofdm_lib; + ldpc; + 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); + % fixed test frame of tx bits rand('seed', 100); - tx_bits = rand(1,Nbitsperframe) > 0.5; + tx_bits = round(rand(1,code_param.data_bits_per_frame)); + tx_codeword = LdpcEncode(tx_bits, code_param.H_rows, code_param.P_matrix); % init logs and BER stats @@ -61,7 +79,9 @@ function ofdm_rx(filename) end prx += states.nin; - [rx_bits states aphase_est_pilot_log arx_np] = ofdm_demod(states, rxbuf_in); + [rx_bits_raw states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in); + rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, arx_np, min(EsNo,30), arx_amp); + rx_bits = rx_codeword(1:code_param.data_bits_per_frame); % measure errors and iterate start machine @@ -70,7 +90,7 @@ function ofdm_rx(filename) next_state = state; if strcmp(state,'searching') - if Nerrs/Nbitsperframe < 0.15 + if Nerrs == 0 next_state = 'synced'; end end @@ -83,7 +103,7 @@ function ofdm_rx(filename) 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("ct_est: %4d foff_est: %3.1f\n", ct_est, foff_est); + 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 @@ -109,11 +129,11 @@ function ofdm_rx(filename) Terrs += Nerrs; Nerrs_log(f) = Nerrs; - Tbits += Nbitsperframe; + Tbits += code_param.data_bits_per_frame; end end - printf("BER: %5.4f Tbits: %d Terrs: %d\n", Terrs/Tbits, Tbits, Terrs); + printf("Coded BER: %5.4f Tbits: %d Terrs: %d\n", Terrs/Tbits, Tbits, Terrs); figure(1); clf; plot(rx_np,'+'); diff --git a/codec2-dev/octave/ofdm_tx.m b/codec2-dev/octave/ofdm_tx.m index 4afba9a2..a414fc4a 100644 --- a/codec2-dev/octave/ofdm_tx.m +++ b/codec2-dev/octave/ofdm_tx.m @@ -6,7 +6,6 @@ #{ TODO: - [ ] Optional LDPC code [ ] measure and report raw and coded BER [ ] maybe 10s worth of frames, sync up to any one automatically + or start with maybe 10 frames @@ -14,14 +13,31 @@ [ ] 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, EbNodB=100, channel='awgn', freq_offset_Hz=0) ofdm_lib; + ldpc; + + % 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 % todo: maybe extend this to 4 or 8 frames, one is a bit short @@ -29,35 +45,32 @@ function ofdm_tx(filename, Nsec, EbNodB=100, channel='awgn', freq_offset_Hz=0) Nframes = floor((Nrows-1)/Ns); rand('seed', 100); - tx_bits = rand(1,Nbitsperframe) > 0.5; + tx_bits = round(rand(1,code_param.data_bits_per_frame)); + tx_codeword = LdpcEncode(tx_bits, code_param.H_rows, code_param.P_matrix); tx = []; for f=1:Nframes - tx = [tx ofdm_mod(states, tx_bits)]; + tx = [tx ofdm_mod(states, tx_codeword)]; end Nsam = length(tx); % channel simulation - EbNo = bps * (10 .^ (EbNodB/10)); - variance = 1/(M*EbNo/2); + EsNo = rate * bps * (10 .^ (EbNodB/10)); + variance = 1/(M*EsNo/2); woffset = 2*pi*freq_offset_Hz/Fs; - % SNR calculation is probably a bit off or confused, need to compare - % to cohpsk_frame_design.ods and make sure I've got it right. PLus - % 3 is used as we are dealing with real channels because mumble - % mumble not understood very well mumble magic number. - - SNRdB = EbNodB + 10*log10(bps*Rs*Nc/Fs) + 3; - printf("EbNo: %3.1f dB SNR(3k): %3.1f dB foff: %3.1fHz\n", EbNodB, SNRdB, freq_offset_Hz); + SNRdB = EbNodB + 10*log10(700/3000); + printf("EbNo: %3.1f dB SNR(3k) est: %3.1f dB foff: %3.1fHz\n", 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.0; path_delay_ms = 1; + 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); @@ -66,22 +79,26 @@ function ofdm_tx(filename, Nsec, EbNodB=100, channel='awgn', freq_offset_Hz=0) randn('seed',1); - spread1 = doppler_spread(dopplerSpreadHz, Fs, (Nsec*(M+Ncp)/M+0.2)*Fs); - spread2 = doppler_spread(dopplerSpreadHz, Fs, (Nsec*(M+Ncp)/M+0.2)*Fs); + 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"); - - hf_gain = 1.0/sqrt(var(spread1)+var(spread2)); end rx = tx; if strcmp(channel, 'hf') - rx = hf_gain * tx(1:Nsam) .* spread1(1:Nsam); - rx += hf_gain * [zeros(1,path_delay_samples) tx(1:Nsam-path_delay_samples)] .* spread2(1:Nsam); + 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)); -- 2.25.1