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