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
states.code_param = code_param;
elseif diversity_en
- rate = 0.5;
-
+ rate = 0.5;
else
rate = 1;
end
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,
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);
% 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
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;
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
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;
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;');
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])
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])
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
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
end
end
- tx = ofdm_tx(states, tx_sym_lin);
+ tx = ofdm_txframe(states, tx_sym_lin);
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
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
-
#{
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)
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
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
next_state = state;
if strcmp(state,'searching')
- if Nerrs/Nbitsperframe < 0.15
+ if Nerrs == 0
next_state = 'synced';
end
end
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
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,'+');
#{
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
[ ] 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
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);
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));
projects / freetel-svn-tracking.git / commitdiff