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