% QPSK modem simulation, initially based on code by Bill Cowley\r
% Generates curves BER versus E/No curves for different modems.\r
% Design to test coherent demodulation ideas on HF channels without\r
-% building a full blown modem. Lacks timing estimation, frame sync.\r
-\r
-% TODO\r
-% [ ] put spreading sample files sonewhere useful\r
+% building a full blown modem. Uses 'genie provided' estimates for\r
+% timing estimation, frame sync.\r
\r
1;\r
\r
function sim_out = ber_test(sim_in, modulation)\r
Fs = 8000;\r
\r
+ verbose = sim_in.verbose;\r
framesize = sim_in.framesize;\r
Ntrials = sim_in.Ntrials;\r
Esvec = sim_in.Esvec;\r
Rs = sim_in.Rs;\r
hf_sim = sim_in.hf_sim;\r
Nhfdelay = floor(sim_in.hf_delay_ms*1000/Fs);\r
+ hf_phase_only = sim_in.hf_phase_only;\r
+ hf_mag_only = sim_in.hf_mag_only;\r
\r
bps = 2;\r
Nsymb = framesize/bps;\r
prev_sym_rx = qpsk_mod([0 0]);\r
rx_symb_log = [];\r
\r
- Np = 5;\r
- r_delay_line = zeros(1,Np);\r
- s_delay_line = zeros(1,Np);\r
+ Np = sim_in.Np; % number of pilot symbols to use in phase est\r
+ Ns = sim_in.Ns; % spacing of pilot symbols, so (Nps-1) data symbols between every pilot\r
+ Nps = Np*Ns;\r
+ r_delay_line = zeros(1,Nps+1);\r
+ s_delay_line = zeros(1,Nps+1);\r
spread_main_phi = 0;\r
spread_delay_phi = 0;\r
spread_main_phi_log = [];\r
\r
- % convert "spreading" samples from 1kHz carrier at Fs to complex baseband\r
+ ldpc_code = sim_in.ldpc_code;\r
+\r
+ if ldpc_code\r
+ % Start CML library\r
+\r
+ currentdir = pwd;\r
+ addpath '/home/david/tmp/cml/mat' % assume the source files stored here\r
+ cd /home/david/tmp/cml\r
+ CmlStartup % note that this is not in the cml path!\r
+ cd(currentdir)\r
+ \r
+ % Our LDPC library\r
+\r
+ ldpc;\r
\r
+ rate = 1/2; \r
+ mod_order = 4; \r
+ modulation = 'QPSK';\r
+ mapping = 'gray';\r
+\r
+ demod_type = 0;\r
+ decoder_type = 0;\r
+ max_iterations = 100;\r
+\r
+ code_param = ldpc_init(rate, framesize, modulation, mod_order, mapping);\r
+ code_param.code_bits_per_frame = framesize;\r
+ code_param.symbols_per_frame = framesize/bps;\r
+ else\r
+ rate = 1;\r
+ end\r
+\r
+ % convert "spreading" samples from 1kHz carrier at Fs to complex\r
+ % baseband, generated by passing a 1kHz sine wave through PathSim\r
+ % with the ccir-poor model, enabling one path at a time.\r
+ \r
Fc = 1000;\r
fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb");\r
spread1k = fread(fspread, "int16")/10000;\r
spread = filter(b,1,spreadbb);\r
spread_2ms = filter(b,1,spreadbb_2ms);\r
\r
- % Determine "gain" of HF channel model, so we can normalise later\r
+ % discard first 1000 samples as these were near 0, probably as\r
+ % PathSim states were ramping up\r
\r
- hf_gain = 1.0/sqrt(var(spread)+var(spread_2ms));\r
+ spread = spread(1000:length(spread));\r
+ spread_2ms = spread_2ms(1000:length(spread_2ms));\r
\r
- sc = 1;\r
+ % Determine "gain" of HF channel model, so we can normalise\r
+ % carrier power during HF channel sim to calibrate SNR. I imagine\r
+ % different implementations of ccir-poor would do this in\r
+ % different ways, leading to different BER results. Oh Well!\r
+\r
+ hf_gain = 1.0/sqrt(var(spread)+var(spread_2ms));\r
\r
% design root nyquist (root raised cosine) filter and init tx and rx filter states\r
\r
alpha = 0.5; T=1/Fs; Nfiltsym=7; M=Fs/Rs;\r
if floor(Fs/Rs) != Fs/Rs\r
printf("oversampling ratio must be an integer\n");\r
- exit;\r
+ return;\r
end\r
hrn = gen_rn_coeffs(alpha, T, Rs, Nfiltsym, M);\r
Nfilter = length(hrn);\r
% No = CFs/(Rs(Es/No))\r
\r
variance = Fs/(Rs*EsNo);\r
- Terrs = 0; Tbits = 0; Ferrs = 0;\r
- printf("EsNo (dB): %f EsNo: %f variance: %f\n", Es, EsNo, variance);\r
+ Terrs = 0; Tbits = 0; Terrsldpc = 0; Tbitsldpc = 0; Ferrsldpc = 0;\r
+ if verbose > 1\r
+ printf("EsNo (dB): %f EsNo: %f variance: %f\n", Es, EsNo, variance);\r
+ end\r
+ \r
+ % init HF channel\r
+ sc = 1;\r
\r
tx_filt_log = [];\r
rx_filt_log = [];\r
rx_baseband_log = [];\r
tx_baseband_log = [];\r
noise_log = [];\r
- c_est_log = [];\r
- c_est = 0;\r
-\r
+ hf_angle_log = [];\r
tx_phase = rx_phase = 0;\r
+ tx_data_buffer = zeros(1,2*framesize);\r
+ s_data_buffer = zeros(1,2*Nsymb);\r
\r
for nn = 1: Ntrials\r
\r
- tx_bits = round( rand( 1, framesize ) );\r
+ %tx_bits = round( rand( 1, framesize*rate ) );\r
+ tx_bits = [1 0 zeros(1,framesize*rate-2)];\r
\r
% modulate\r
\r
- s = zeros(1, Nsymb);\r
- for i=1:Nsymb\r
- tx_symb = qpsk_mod(tx_bits(2*(i-1)+1:2*i));\r
- %printf("shift: %f prev_sym: %f ", tx_symb, prev_sym_tx);\r
- if strcmp(modulation,'dqpsk')\r
- tx_symb *= prev_sym_tx;\r
- %printf("tx_symb: %f\n", tx_symb);\r
- prev_sym_tx = tx_symb;\r
- end \r
- s(i) = tx_symb;\r
+ if ldpc_code\r
+ [tx_bits, s] = ldpc_enc(tx_bits, code_param);\r
+ t2 = tx_bits;\r
+ s2 = s;\r
+ else\r
+ s = zeros(1, Nsymb);\r
+ for i=1:Nsymb\r
+ tx_symb = qpsk_mod(tx_bits(2*(i-1)+1:2*i));\r
+ %printf("shift: %f prev_sym: %f ", tx_symb, prev_sym_tx);\r
+ if strcmp(modulation,'dqpsk')\r
+ tx_symb *= prev_sym_tx;\r
+ %printf("tx_symb: %f\n", tx_symb);\r
+ prev_sym_tx = tx_symb;\r
+ end \r
+ s(i) = tx_symb;\r
+ end\r
end\r
-\r
s_ch = s;\r
\r
% root nyquist filter symbols\r
% HF channel simulation\r
\r
if hf_sim\r
+\r
hf_sim_delay_line(1:Nhfdelay) = hf_sim_delay_line(M+1:M+Nhfdelay);\r
hf_sim_delay_line(Nhfdelay+1:M+Nhfdelay) = tx_filt;\r
\r
- if ((sc+M) > length(spread)) || ((sc+M) > length(spread_2ms))\r
- sc =1 ;\r
+ % disable as a wrap around will cause a nasty phase jump. Best to generate\r
+ % longer files\r
+ %if ((sc+M) > length(spread)) || ((sc+M) > length(spread_2ms))\r
+ % sc =1 ;\r
+ %end\r
+ if hf_phase_only\r
+ comb = conj(spread(sc:sc+M-1))' + conj(spread_2ms(sc:sc+M-1))';\r
+ tx_filt = tx_filt.*exp(j*angle(comb));\r
+ hf_angle_log = [hf_angle_log angle(comb)];\r
+ else\r
+ if hf_mag_only\r
+ comb = conj(spread(sc:sc+M-1))' + conj(spread_2ms(sc:sc+M-1))';\r
+ tx_filt = tx_filt.*abs(comb); \r
+ else\r
+ % regular HF channel sim\r
+ tx_filt = tx_filt.*conj(spread(sc:sc+M-1))' + hf_sim_delay_line(1:M).*conj(spread_2ms(sc:sc+M-1))';\r
+ end\r
end\r
- tx_filt = tx_filt.*spread(sc:sc+M-1)' + hf_sim_delay_line(1:M).*spread_2ms(sc:sc+M-1)';\r
sc += M;\r
-\r
+ \r
% normalise so average HF power C=1\r
\r
- tx_filt *= hf_gain;\r
+ if hf_phase_only == 0 % C already 1 if we are just tweaking phase\r
+ tx_filt *= hf_gain;\r
+ end\r
end\r
tx_filt_log = [tx_filt_log tx_filt];\r
\r
s_ch(k) = rx_filt; \r
end\r
\r
- % coherent demod phase estimation and correction\r
- % need sliding window\r
- % use known (pilot) symbols\r
- % start with all symbols pilots, then gradually decimate, e.g. 1 in 5 pilots\r
+ % coherent demod phase estimation and correction using pilot symbols\r
+ % cheating a bit here, we use fact that all tx-ed symbols are known\r
\r
if phase_est\r
for i=1:Nsymb\r
\r
% delay line for phase est window\r
\r
- r_delay_line(1:Np-1) = r_delay_line(2:Np);\r
- r_delay_line(Np) = s_ch(i);\r
+ r_delay_line(1:Nps-1) = r_delay_line(2:Nps);\r
+ r_delay_line(Nps) = s_ch(i);\r
\r
% delay in tx data to compensate data for phase est window\r
\r
- s_delay_line(1:Np-1) = s_delay_line(2:Np);\r
- s_delay_line(Np) = s(i);\r
- tx_bits(2*(i-1)+1:2*i) = qpsk_demod(s_delay_line(floor(Np/2)+1));\r
+ s_delay_line(1:Nps-1) = s_delay_line(2:Nps);\r
+ s_delay_line(Nps) = s(i);\r
+ tx_bits(2*(i-1)+1:2*i) = qpsk_demod(s_delay_line(floor(Nps/2)+1));\r
\r
- % estimate phase from known symbols and correct\r
+ % estimate phase from surrounding known pilot symbols and correct\r
\r
- corr = s_delay_line * r_delay_line';\r
- s_ch(i) = r_delay_line(floor(Np/2)+1).*exp(j*angle(corr));\r
+ corr = 0; centre = floor(Nps/2)+1;\r
+ for k=1:Ns:(Nps+1)\r
+ if (k != centre)\r
+ corr += s_delay_line(k) * r_delay_line(k)';\r
+ end\r
+ end\r
+ s_ch(i) = r_delay_line(centre).*exp(j*angle(corr));\r
end \r
- %printf("corr: %f angle: %f\n", corr, angle(corr));\r
+ %printf("corr: %f angle: %f\n", corr, angle(corr));\r
end\r
\r
% de-modulate\r
\r
% Measure BER\r
\r
- % discard bits from first 2*Nfiltsym symbols as tx and rx filter memories not full\r
+ % discard bits from first 2*Nfiltsym+Nps+1 symbols as tx\r
+ % and rx filter and phase est memories not full\r
\r
+ skip = bps*(2*Nfiltsym+1+Nps+1);\r
if nn == 1\r
- tx_bits = tx_bits(2*bps*Nfiltsym+1:length(tx_bits));\r
- rx_bits = rx_bits(2*bps*Nfiltsym+1:length(rx_bits));\r
+ tx_bits_tmp = tx_bits(skip:length(tx_bits));\r
+ rx_bits_tmp = rx_bits(skip:length(rx_bits));\r
end\r
\r
- error_positions = xor( rx_bits, tx_bits );\r
+ error_positions = xor( rx_bits_tmp, tx_bits_tmp );\r
Nerrs = sum(error_positions);\r
Terrs += Nerrs;\r
- Tbits = Tbits + length(tx_bits);\r
+ Tbits += length(tx_bits_tmp);\r
+\r
+ % Optionally LDPC decode\r
+ \r
+ if ldpc_code\r
+ % filter memories etc screw up frame alignment so we need to buffer a frame\r
+\r
+ tx_data_buffer(1:framesize) = tx_data_buffer(framesize+1:2*framesize);\r
+ s_data_buffer(1:Nsymb) = s_data_buffer(Nsymb+1:2*Nsymb);\r
+ tx_data_buffer(framesize+1:2*framesize) = tx_bits;\r
+ s_data_buffer(Nsymb+1:2*Nsymb) = s_ch;\r
+\r
+ st_tx = (Nfiltsym-1)*bps+1;\r
+ st_s = (Nfiltsym-1);\r
+\r
+ detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, ...\r
+ s_data_buffer(st_s+1:st_s+Nsymb), min(100,EsNo));\r
+\r
+ % ignore first frame as filter, phase est memories filling up\r
+ if nn != 1\r
+ error_positions = xor( detected_data(1:framesize*rate), ...\r
+ tx_data_buffer(st_tx:st_tx+framesize*rate-1) );\r
+ Nerrs = sum(error_positions);\r
+ if Nerrs\r
+ Ferrsldpc++;\r
+ end\r
+ Terrsldpc += Nerrs;\r
+ Tbitsldpc += framesize*rate;\r
+ end\r
+ end\r
\r
end\r
\r
TERvec(ne) = Terrs;\r
- FERvec(ne) = Ferrs;\r
BERvec(ne) = Terrs/Tbits;\r
- %printf(" Terrs: %d BER %f BER theory %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs,\r
- % Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), var(tx_filt_log), var(noise_log),\r
- % var(tx_filt_log)/Rs, var(noise_log)/Fs, (var(tx_filt_log)/Rs)/(var(noise_log)/Fs));\r
+ if ldpc_code\r
+ TERldpcvec(ne) = Terrsldpc;\r
+ FERldpcvec(ne) = Ferrsldpc;\r
+ BERldpcvec(ne) = Terrsldpc/Tbitsldpc;\r
+ end\r
+\r
+ if verbose \r
+ printf("EsNo (dB): %f Terrs: %d BER %f BER theory %f", Es, Terrs,\r
+ Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)));\r
+ if ldpc_code\r
+ printf(" LDPC: Terrs: %d BER: %f", Terrsldpc, Terrsldpc/Tbitsldpc);\r
+ end\r
+ printf("\n");\r
+ end\r
+ if verbose > 1\r
+ printf("Terrs: %d BER %f BER theory %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs,\r
+ Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), var(tx_filt_log), var(noise_log),\r
+ var(tx_filt_log)/Rs, var(noise_log)/Fs, (var(tx_filt_log)/Rs)/(var(noise_log)/Fs));\r
+ end\r
end\r
\r
Ebvec = Esvec - 10*log10(bps);\r
- sim_out.BERvec = BERvec;\r
- sim_out.BER_theoryvec = 0.5*erfc(sqrt(10.^(Ebvec/10)));\r
- sim_out.Ebvec = Ebvec;\r
- sim_out.FERvec = FERvec;\r
- sim_out.TERvec = TERvec;\r
+ sim_out.BERvec = BERvec;\r
+ sim_out.Ebvec = Ebvec;\r
+ sim_out.TERvec = TERvec;\r
+ if ldpc_code\r
+ sim_out.BERldpcvec = BERldpcvec;\r
+ sim_out.TERldpcvec = TERldpcvec;\r
+ sim_out.FERldpcvec = FERldpcvec;\r
+ end\r
\r
if plot_scatter\r
figure(2);\r
clf;\r
scat = rx_symb_log(2*Nfiltsym:length(rx_symb_log)) .* exp(j*pi/4);\r
plot(real(scat), imag(scat),'+');\r
- figure(3);\r
- plot(c_est_log);\r
+\r
+ if hf_sim\r
+ figure(3);\r
+\r
+ if hf_phase_only\r
+ plot(hf_angle_log);\r
+ axis([1 10000 min(hf_angle_log) max(hf_angle_log)])\r
+ else\r
+ plot(abs(spread));\r
+ hold on;\r
+ plot(abs(spread_2ms),'g');\r
+ hold off;\r
+ axis([1 10000 0 1.4])\r
+ end\r
+ end\r
end\r
endfunction\r
\r
% Gray coded QPSK demodulation function\r
\r
function two_bits = qpsk_demod(symbol)\r
+ if isscalar(symbol) == 0\r
+ printf("only works with scalars\n");\r
+ return;\r
+ end\r
bit0 = real(symbol*exp(j*pi/4)) < 0;\r
bit1 = imag(symbol*exp(j*pi/4)) < 0;\r
two_bits = [bit1 bit0];\r
\r
% Start simulations ---------------------------------------\r
\r
-sim_in.Esvec = 1:12; \r
-sim_in.Ntrials = 100;\r
-sim_in.framesize = 100;\r
-sim_in.Rs = 100;\r
+more off;\r
+sim_in.verbose = 1;\r
+\r
+sim_in.Esvec = 1:5; \r
+sim_in.Ntrials = 3;\r
+sim_in.framesize = 576;\r
+sim_in.Rs = 400;\r
sim_in.phase_offset = 0;\r
sim_in.phase_est = 0;\r
sim_in.w_offset = 0;\r
sim_in.plot_scatter = 0;\r
+sim_in.hf_delay_ms = 2;\r
sim_in.hf_sim = 0;\r
+sim_in.Np = 6;\r
+sim_in.Ns = 5;\r
+sim_in.hf_phase_only = 0;\r
+sim_in.hf_mag_only = 0;\r
+sim_in.ldpc_code = 1;\r
\r
-sim_qpsk = ber_test(sim_in, 'qpsk');\r
-sim_dqpsk = ber_test(sim_in, 'dqpsk');\r
+Ebvec = sim_in.Esvec - 10*log10(2);\r
+BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10)));\r
+\r
+sim_qpsk = ber_test(sim_in, 'qpsk');\r
+\r
+figure(1); \r
+clf;\r
+semilogy(Ebvec, BER_theory,'r;QPSK theory;')\r
+hold on;\r
+semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec,'g;QPSK;')\r
+hold off;\r
+xlabel('Eb/N0')\r
+ylabel('BER')\r
+grid("minor")\r
+\r
+\r
+if 0\r
+sim_in.hf_mag_only = 1;\r
+sim_qpsk_mag = ber_test(sim_in, 'qpsk');\r
+\r
+sim_in.hf_mag_only = 0;\r
+sim_in.hf_phase_only = 1;\r
sim_in.phase_est = 1;\r
-sim_qpsk_coh = ber_test(sim_in, 'qpsk');\r
-sim_in.hf_delay_ms = 2;\r
-sim_in.hf_sim = 1;\r
-sim_qpsk_hf = ber_test(sim_in, 'qpsk');\r
+sim_qpsk_phase = ber_test(sim_in, 'qpsk');\r
+\r
+sim_in.hf_phase_only = 0;\r
+sim_qpsk_coh_6_5 = ber_test(sim_in, 'qpsk');\r
+\r
sim_in.phase_est = 0;\r
-sim_dqpsk_hf = ber_test(sim_in, 'dqpsk');\r
+sim_dqpsk = ber_test(sim_in, 'dqpsk');\r
\r
figure(1); \r
clf;\r
-semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec,'b;qpsk;')\r
+semilogy(Ebvec, BER_theory,'r;QPSK theory;')\r
hold on;\r
-semilogy(sim_qpsk.Ebvec, sim_qpsk.BER_theoryvec,'r;qpsk theory;')\r
-semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'g;dqpsk;')\r
-semilogy(sim_qpsk_coh.Ebvec, sim_qpsk_coh.BERvec,'k;qpsk coh;')\r
-semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'c;qpsk hf;')\r
-semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'m;dqpsk hf;')\r
+semilogy(sim_qpsk_mag.Ebvec, sim_qpsk_mag.BERvec,'g;QPSK CCIR poor mag;')\r
+semilogy(sim_qpsk_phase.Ebvec, sim_qpsk_phase.BERvec,'k;QPSK CCIR poor phase;')\r
+semilogy(sim_qpsk_coh_6_5.Ebvec, sim_qpsk_coh_6_5.BERvec,'c;QPSK CCIR poor Np=6 Ns=5;')\r
+semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'b;DQPSK CCIR poor;')\r
+%semilogy(sim_qpsk_coh_5_24.Ebvec, sim_qpsk_coh_5_24.BERvec,'k;QPSK Ns=5 Np=24;')\r
+%semilogy(sim_qpsk_coh_2_12.Ebvec, sim_qpsk_coh_2_12.BERvec,'c;QPSK Ns=2 Np=12;')\r
hold off;\r
xlabel('Eb/N0')\r
ylabel('BER')\r
grid("minor")\r
-\r
+axis([min(Ebvec)-1 max(Ebvec)+1 1E-2 1])\r
+end\r
projects / freetel-svn-tracking.git / commitdiff