From: drowe67 Date: Tue, 18 Feb 2014 09:07:36 +0000 (+0000) Subject: getting good results, compared to current DQPSK FreeDV, 6dB gain using rate 3/4 and... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=7d4eb4e318a1bed18fa629e9734e18fc77b7f64a;p=freetel-svn-tracking.git getting good results, compared to current DQPSK FreeDV, 6dB gain using rate 3/4 and 10dB gain using rate 1/2 git-svn-id: https://svn.code.sf.net/p/freetel/code@1405 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/ldpc.m b/codec2-dev/octave/ldpc.m index e69d1f27..7ce4a87a 100644 --- a/codec2-dev/octave/ldpc.m +++ b/codec2-dev/octave/ldpc.m @@ -144,8 +144,8 @@ function [codeword s] = ldpc_enc(data, code_param) endfunction -function detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r, EsNo) - symbol_likelihood = Demod2D( r, code_param.S_matrix, EsNo); +function detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r, EsNo, fading) + 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 ); diff --git a/codec2-dev/octave/test_qpsk2.m b/codec2-dev/octave/test_qpsk2.m new file mode 100644 index 00000000..9052837f --- /dev/null +++ b/codec2-dev/octave/test_qpsk2.m @@ -0,0 +1,408 @@ +% test_qps2k.m +% David Rowe Feb 2014 +% +% QPSK modem simulation, version 2. Simplifed version of +% test_qpsk. 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. Uses 'genie provided' estimates for timing estimation, +% frame sync. + +1; + +% main test function + +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; + phase_offset = sim_in.phase_offset; + phase_est = sim_in.phase_est; + w_offset = sim_in.w_offset; + plot_scatter = sim_in.plot_scatter; + Rs = sim_in.Rs; + hf_sim = sim_in.hf_sim; + nhfdelay = sim_in.hf_delay_ms*Rs/1000; + hf_phase_only = sim_in.hf_phase_only; + hf_mag_only = sim_in.hf_mag_only; + Nc = sim_in.Nc; + + bps = 2; + Nsymb = framesize/bps; + prev_sym_tx = qpsk_mod([0 0]); + prev_sym_rx = qpsk_mod([0 0]); + + 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; + + ldpc_code = sim_in.ldpc_code; + + % Init LDPC -------------------------------------------------------------------- + + 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 = sim_in.ldpc_code_rate; + 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 + + % Init HF channel model from stored sample files of spreading signal ---------------------------------- + + % 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; M = Fs/Rs; + fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb"); + spread1k = fread(fspread, "int16")/10000; + fclose(fspread); + fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb"); + spread1k_2ms = fread(fspread, "int16")/10000; + fclose(fspread); + + % down convert to complex baseband + spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))'); + spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))'); + + % remove -2000 Hz image + b = fir1(50, 5/Fs); + spread = filter(b,1,spreadbb); + spread_2ms = filter(b,1,spreadbb_2ms); + + % discard first 1000 samples as these were near 0, probably as + % PathSim states were ramping up + + spread = spread(1000:length(spread)); + spread_2ms = spread_2ms(1000:length(spread_2ms)); + + % decimate down to Rs + + spread = spread(1:M:length(spread)); + spread_2ms = spread_2ms(1:M:length(spread_2ms)); + + % 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)); + + % Start Simulation ---------------------------------------------------------------- + + for ne = 1:length(Esvec) + EsNodB = Esvec(ne); + EsNo = 10^(EsNodB/10); + + variance = 1/EsNo; + if verbose > 1 + printf("EsNo (dB): %f EsNo: %f variance: %f\n", EsNodB, EsNo, variance); + end + + Terrs = 0; Tbits = 0; Terrsldpc = 0; Tbitsldpc = 0; Ferrsldpc = 0; + + tx_symb_log = []; + rx_symb_log = []; + noise_log = []; + + % init HF channel + + hf_n = 1; + hf_angle_log = []; + hf_fading = ones(1,Nsymb); % default input for ldpc dec + hf_model = ones(Ntrials*Nsymb/Nc, Nc); % defaults for plotting surface + + for nn = 1: Ntrials + + tx_bits = round( rand( 1, framesize*rate ) ); + + % modulate -------------------------------------------- + + if ldpc_code + [tx_bits, s] = ldpc_enc(tx_bits, code_param); + else + s = zeros(1, Nsymb); + for i=1:Nsymb + tx_symb = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); + if strcmp(modulation,'dqpsk') + tx_symb *= prev_sym_tx; + prev_sym_tx = tx_symb; + end + s(i) = tx_symb; + end + end + s_ch = s; + + % HF channel simulation ------------------------------------ + + if hf_sim + + % separation between carriers. Note this is + % effectively under samples at Rs, I dont think this + % matters. Equivalent to doing freq shift at Fs, then + % decimating to Rs. + + wsep = 2*pi*(1+0.5); % e.g. 75Hz spacing at Rs=50Hz, alpha=0.5 filters + + if Nsymb/Nc != floor(Nsymb/Nc) + printf("Error: Nsymb/Nc must be an integrer\n") + return; + end + + % arrange symbols in Nsymb/Nc by Nc matrix + + for i=1:Nc:Nsymb + + % Determine HF channel at each carrier for this symbol + + for k=1:Nc + hf_model(hf_n, k) = hf_gain*(spread(hf_n) + exp(-j*k*wsep*nhfdelay)*spread_2ms(hf_n)); + hf_fading(i+k-1) = abs(hf_model(hf_n, k)); + if hf_mag_only + s_ch(i+k-1) *= abs(hf_model(hf_n, k)); + else + s_ch(i+k-1) *= hf_model(hf_n, k); + end + end + hf_n++; + end + end + + tx_symb_log = [tx_symb_log s_ch]; + + % AWGN noise and phase/freq offset channel simulation + % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im + + noise = sqrt(variance*0.5)*(randn(1,Nsymb) + j*randn(1,Nsymb)); + noise_log = [noise_log noise]; + s_ch = s_ch.*exp(j*phase_offset) + noise; + phase_offset += w_offset; + + % de-modulate + + rx_bits = zeros(1, framesize); + for i=1:Nsymb + rx_symb = s_ch(i); + if strcmp(modulation,'dqpsk') + tmp = rx_symb; + rx_symb *= conj(prev_sym_rx/abs(prev_sym_rx)); + prev_sym_rx = tmp; + end + rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb); + rx_symb_log = [rx_symb_log rx_symb]; + end + + % Measure BER + + error_positions = xor(rx_bits, tx_bits); + Nerrs = sum(error_positions); + Terrs += Nerrs; + Tbits += length(tx_bits); + + % Optionally LDPC decode + + if ldpc_code + detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, s_ch, min(100,EsNo), hf_fading); + error_positions = xor( detected_data(1:framesize*rate), tx_bits(1:framesize*rate) ); + Nerrs = sum(error_positions); + if Nerrs + Ferrsldpc++; + end + Terrsldpc += Nerrs; + Tbitsldpc += framesize*rate; + end + end + + TERvec(ne) = Terrs; + BERvec(ne) = Terrs/Tbits; + 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", EsNodB, Terrs, + Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2))); + if ldpc_code + printf(" LDPC: Terrs: %d BER: %f Ferrs: %d FER: %f", + Terrsldpc, Terrsldpc/Tbitsldpc, Ferrsldpc, Ferrsldpc/(Ntrials-1)); + 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_symb_log), var(noise_log), + var(tx_symb_log), var(noise_log), var(tx_symb_log)/var(noise_log)); + end + end + + Ebvec = Esvec - 10*log10(bps); + 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 .* exp(j*pi/4); + plot(real(scat), imag(scat),'+'); + + figure(3); + clf; + + y = 1:Rs*2; + x = 1:Nc; + EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance); + mesh(x,y,EsNodBSurface); + grid + %axis([1 Nc 1 Rs*2 -10 10]) + + figure(4); + clf; + mesh(x,y,unwrap(angle(hf_model(y,:)))); + + figure(5); + clf; + subplot(211) + plot(real(spread)); + hold on; + plot(imag(spread),'g'); + hold off; + subplot(212) + plot(real(spread_2ms)); + hold on; + plot(imag(spread_2ms),'g'); + hold off; + + figure(6) + tmp = []; + for i = 1:hf_n-1 + tmp = [tmp abs(hf_model(i,:))]; + end + hist(tmp); + end + +endfunction + +% Gray coded QPSK modulation function + +function symbol = qpsk_mod(two_bits) + two_bits_decimal = sum(two_bits .* [2 1]); + switch(two_bits_decimal) + case (0) symbol = 1; + case (1) symbol = j; + case (2) symbol = -j; + case (3) symbol = -1; + endswitch +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]; +endfunction + +% Start simulations --------------------------------------- + +more off; +sim_in.verbose = 2; +sim_in.plot_scatter = 1; + +sim_in.Esvec = 5; +sim_in.Ntrials = 30; +sim_in.framesize = 576; +sim_in.Rs = 100; +sim_in.Nc = 8; + +sim_in.phase_offset = 0; +sim_in.w_offset = 0; + +sim_in.hf_delay_ms = 2; +sim_in.hf_sim = 1; +sim_in.hf_phase_only = 0; +sim_in.hf_mag_only = 1; + +sim_in.phase_est = 0; +sim_in.Np = 6; +sim_in.Ns = 5; + +sim_in.ldpc_code_rate = 3/4; +sim_in.ldpc_code = 1; + +sim_qpsk_hf = ber_test(sim_in, 'qpsk'); + +sim_in.hf_sim = 0; +sim_in.plot_scatter = 0; +sim_in.Esvec = 2:15; +sim_in.ldpc_code = 0; +Ebvec = sim_in.Esvec - 10*log10(2); +BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); +sim_qpsk = ber_test(sim_in, 'qpsk'); +sim_dqpsk = ber_test(sim_in, 'dqpsk'); +sim_in.hf_sim = 1; +sim_in.Esvec = 2:15; +sim_qpsk_hf = ber_test(sim_in, 'qpsk'); +sim_dqpsk_hf = ber_test(sim_in, 'dqpsk'); +sim_in.ldpc_code = 1; +sim_qpsk_hf_ldpc1 = ber_test(sim_in, 'qpsk'); +sim_in.ldpc_code_rate = 1/2; +sim_qpsk_hf_ldpc2 = ber_test(sim_in, 'qpsk'); +sim_in.ldpc_code_rate = 3/4; +sim_in.hf_sim = 0; +sim_qpsk_awgn_ldpc = 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 AWGN;') +semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'r;QPSK HF;') +semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;') +semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'m;DQPSK HF;') +semilogy(sim_qpsk_hf_ldpc1.Ebvec, sim_qpsk_hf_ldpc1.BERldpcvec,'k;QPSK HF LDPC 3/4;') +semilogy(sim_qpsk_hf_ldpc2.Ebvec, sim_qpsk_hf_ldpc2.BERldpcvec,'b;QPSK HF LDPC 1/2;') +semilogy(sim_qpsk_awgn_ldpc.Ebvec, sim_qpsk_awgn_ldpc.BERldpcvec,'k;QPSK AWGN LDPC 3/4;') + +semilogy([min(Ebvec) max(Ebvec)], [3E-2 3E-2], 'r;Codec 2 target BER;') +hold off; +xlabel('Eb/N0') +ylabel('BER') +grid("minor") +axis([min(Ebvec) max(Ebvec) 1E-2 1])