--- /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
+ newldpc = sim_in.newldpc;\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
+ phase_est_method = sim_in.phase_est_method;\r
+ if phase_est_method == 2\r
+ Np = sim_in.Np;\r
+ Ns = sim_in.Ns;\r
+ if Np/2 == floor(Np/2)\r
+ printf("Np must be odd\n");\r
+ return;\r
+ end\r
+ Nps = (Np-1)*Ns+1; \r
+ else\r
+ Nps = sim_in.Np; \r
+ end\r
+ r_delay_line = zeros(Nc, Nps);\r
+ s_delay_line = zeros(Nc, Nps);\r
+ ph_est_log = [];\r
+\r
+ phase_noise_amp = sim_in.phase_noise_amp;\r
+\r
+ ldpc_code = sim_in.ldpc_code;\r
+\r
+ tx_bits_buf = zeros(1,2*framesize);\r
+ rx_bits_buf = zeros(1,2*framesize);\r
+ rx_symb_buf = zeros(1,2*Nsymb);\r
+ hf_fading_buf = zeros(1,2*Nsymb);\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
+ mod_strip_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
+ tx_bits_buf(1:framesize) = tx_bits_buf(framesize+1:2*framesize);\r
+ tx_bits_buf(framesize+1:2*framesize) = tx_bits;\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
+ phase_noise = phase_noise_amp*(2.0*rand(1,Nsymb)-1.0);\r
+ \r
+ % organise into carriers to apply frequency and phase offset\r
+\r
+ for i=1:Nc:Nsymb\r
+ for k=1:Nc\r
+ s_ch(i+k-1) = s_ch(i+k-1)*exp(j*(phase_offset+phase_noise(i+k-1))) + noise(i+k-1);\r
+ end \r
+ phase_offset += w_offset;\r
+ end\r
+ \r
+ % phase estimation\r
+ \r
+ ph_est = zeros(Nc,1);\r
+\r
+ if phase_est\r
+\r
+ % organise into carriers\r
+\r
+ for i=1:Nc:Nsymb\r
+\r
+ for k=1:Nc\r
+ centre = floor(Nps/2)+1;\r
+\r
+ % delay line for phase est window\r
+\r
+ r_delay_line(k,1:Nps-1) = r_delay_line(k,2:Nps);\r
+ r_delay_line(k,Nps) = s_ch(i+k-1);\r
+\r
+ % delay in tx data to compensate data for phase est window\r
+\r
+ s_delay_line(k,1:Nps-1) = s_delay_line(k,2:Nps);\r
+ s_delay_line(k,Nps) = s(i+k-1);\r
+ \r
+ if phase_est_method == 1\r
+ % QPSK modulation strip and phase est\r
+\r
+ mod_strip_pol = angle(r_delay_line(k,:)) * 4;\r
+ mod_strip_rect = exp(j*mod_strip_pol);\r
+\r
+ ph_est_pol = atan2(sum(imag(mod_strip_rect)),sum(real(mod_strip_rect)))/4;\r
+ ph_est(k) = exp(j*ph_est_pol);\r
+\r
+ s_ch(i+k-1) = r_delay_line(k,centre).*exp(-j*ph_est_pol);\r
+ % s_ch(i+k-1) = r_delay_line(k,centre);\r
+ end\r
+\r
+ if phase_est_method == 3\r
+ % QPSK modulation strip and phase est with original symbol mags\r
+\r
+ mod_strip_pol = angle(r_delay_line(k,:)) * 4;\r
+ mod_strip_rect = abs(r_delay_line(k,:)) .* exp(j*mod_strip_pol);\r
+\r
+ ph_est_pol = atan2(sum(imag(mod_strip_rect)),sum(real(mod_strip_rect)))/4;\r
+ ph_est(k) = exp(j*ph_est_pol);\r
+\r
+ s_ch(i+k-1) = r_delay_line(k,centre).*exp(-j*ph_est_pol);\r
+ % s_ch(i+k-1) = r_delay_line(k,centre);\r
+ end\r
+\r
+ if phase_est_method == 2\r
+\r
+ % estimate phase from surrounding known pilot symbols and correct\r
+\r
+ corr = 0;\r
+ for m=1:Ns:Nps\r
+ if (m != centre)\r
+ corr += s_delay_line(k,m) * r_delay_line(k,m)';\r
+ end\r
+ end\r
+ ph_est(k) = conj(corr/(1E-6+abs(corr)));\r
+ s_ch(i+k-1) = r_delay_line(k,centre).*exp(j*angle(corr));\r
+ %s_ch(i+k-1) = r_delay_line(k,centre);\r
+ end\r
+\r
+ end\r
+ \r
+ ph_est_log = [ph_est_log ph_est];\r
+ end \r
+ %printf("corr: %f angle: %f\n", corr, angle(corr));\r
+ end\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
+if newldpc\r
+ rx_bits_buf(1:framesize) = rx_bits_buf(framesize+1:2*framesize);\r
+ rx_bits_buf(framesize+1:2*framesize) = rx_bits;\r
+ rx_symb_buf(1:Nsymb) = rx_symb_buf(Nsymb+1:2*Nsymb);\r
+ rx_symb_buf(Nsymb+1:2*Nsymb) = s_ch;\r
+ hf_fading_buf(1:Nsymb) = hf_fading_buf(Nsymb+1:2*Nsymb);\r
+ hf_fading_buf(Nsymb+1:2*Nsymb) = hf_fading;\r
+\r
+ % determine location of start and end of frame depending on processing delays\r
+\r
+ if phase_est\r
+ st_rx_bits = 1+(floor(Nps/2)+1-1)*Nc*2;\r
+ st_rx_symb = 1+(floor(Nps/2)+1-1)*Nc;\r
+ else\r
+ st_rx_bits = 1;\r
+ st_rx_symb = 1;\r
+ end\r
+ en_rx_bits = st_rx_bits+framesize-1;\r
+ en_rx_symb = st_rx_symb+Nsymb-1;\r
+\r
+ if nn > 1\r
+ % Measure BER\r
+\r
+ %printf("nn: %d centre: %d\n", nn, floor(Nps/2)+1);\r
+ %tx_bits_buf(1:20)\r
+ %rx_bits_buf(st_rx_bits:st_rx_bits+20-1)\r
+ error_positions = xor(rx_bits_buf(st_rx_bits:en_rx_bits), tx_bits_buf(1:framesize));\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, ...\r
+ rx_symb_buf(st_rx_symb:en_rx_symb), min(100,EsNo), hf_fading_buf(1:Nsymb));\r
+ error_positions = xor( detected_data(1:framesize*rate), tx_bits_buf(1:framesize*rate) );\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
+else \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
+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(Nps*Nc:length(scat))), imag(scat(Nps*Nc:length(scat))),'+');\r
+ title('Scatter plot');\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
+ title('HF Channel Es/No');\r
+\r
+ figure(4);\r
+ clf;\r
+ %mesh(x,y,unwrap(angle(hf_model(y,:))));\r
+ subplot(211)\r
+ plot(y,abs(hf_model(y,1)))\r
+ title('HF Channel Carrier 1 Mag');\r
+ subplot(212)\r
+ plot(y,angle(hf_model(y,1)))\r
+ title('HF Channel Carrier 1 Phase');\r
+\r
+ if phase_est\r
+ scat = ph_est_log(1,floor(Nps/2):Rs*2+floor(Nps/2)-1);\r
+ hold on;\r
+ plot(angle(scat),'r');\r
+ hold off;\r
+\r
+ figure(5)\r
+ clf;\r
+ scat = ph_est_log(1,y);\r
+ plot(real(scat), imag(scat),'+');\r
+ title('Carrier 1 Phase Est');\r
+ axis([-1 1 -1 1])\r
+ end\r
+if 0 \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
+ 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
+function sim_in = standard_init\r
+ sim_in.verbose = 1;\r
+ sim_in.plot_scatter = 0;\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
+ sim_in.phase_noise_amp = 0;\r
+\r
+ sim_in.hf_delay_ms = 2;\r
+ sim_in.hf_sim = 0;\r
+ sim_in.hf_phase_only = 0;\r
+ sim_in.hf_mag_only = 1;\r
+\r
+ sim_in.phase_est = 0;\r
+ sim_in.phase_est_method = 1;\r
+ sim_in.Np = 5;\r
+ sim_in.Ns = 5;\r
+\r
+ sim_in.ldpc_code_rate = 1/2;\r
+ sim_in.ldpc_code = 1;\r
+endfunction\r
+\r
+function ideal\r
+\r
+ sim_in = standard_init();\r
+\r
+ sim_in.verbose = 1;\r
+ sim_in.plot_scatter = 1;\r
+\r
+ sim_in.Esvec = 5; \r
+ sim_in.hf_sim = 1;\r
+ sim_in.Ntrials = 100;\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
+\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_in.ldpc_code_rate = 3/4;\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
+ hold off;\r
+ xlabel('Eb/N0')\r
+ ylabel('BER')\r
+ grid("minor")\r
+ axis([min(Ebvec) max(Ebvec) 1E-3 1])\r
+endfunction\r
+\r
+function phase_noise\r
+ sim_in = standard_init();\r
+\r
+ sim_in.verbose = 1;\r
+ sim_in.plot_scatter = 1;\r
+\r
+ sim_in.Esvec = 100; \r
+ sim_in.Ntrials = 100;\r
+\r
+ sim_in.ldpc_code_rate = 1/2;\r
+ sim_in.ldpc_code = 1;\r
+\r
+ sim_in.phase_noise_amp = pi/16;\r
+ tmp = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.plot_scatter = 0;\r
+ sim_in.Esvec = 2:8; \r
+ sim_qpsk_hf = ber_test(sim_in, 'qpsk');\r
+\r
+ Ebvec = sim_in.Esvec - 10*log10(2);\r
+ BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10)));\r
+\r
+ sim_in.phase_noise_amp = 0;\r
+ sim_qpsk = ber_test(sim_in, 'qpsk');\r
+ sim_in.phase_noise_amp = pi/8;\r
+ sim_qpsk_pn8 = ber_test(sim_in, 'qpsk');\r
+ sim_in.phase_noise_amp = pi/16;\r
+ sim_qpsk_pn16 = ber_test(sim_in, 'qpsk');\r
+ sim_in.phase_noise_amp = pi/32;\r
+ sim_qpsk_pn32 = ber_test(sim_in, 'qpsk');\r
+\r
+ figure(1); \r
+ clf;\r
+ semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec,'g;QPSK phase noise 0;')\r
+ hold on;\r
+ semilogy(sim_qpsk_pn8.Ebvec, sim_qpsk_pn8.BERvec,'c;QPSK phase noise +/- pi/8;')\r
+ semilogy(sim_qpsk_pn16.Ebvec, sim_qpsk_pn16.BERvec,'b;QPSK phase noise +/- pi/16;')\r
+ semilogy(sim_qpsk_pn32.Ebvec, sim_qpsk_pn32.BERvec,'k;QPSK phase noise +/- pi/32;')\r
+\r
+ semilogy(sim_qpsk.Ebvec, sim_qpsk.BERldpcvec,'g;QPSK phase noise 0 ldpc;')\r
+ semilogy(sim_qpsk_pn8.Ebvec, sim_qpsk_pn8.BERldpcvec,'c;QPSK phase noise +/- pi/8 ldpc;')\r
+ semilogy(sim_qpsk_pn16.Ebvec, sim_qpsk_pn16.BERldpcvec,'b;QPSK phase noise +/- pi/16 ldpc;')\r
+ semilogy(sim_qpsk_pn32.Ebvec, sim_qpsk_pn32.BERldpcvec,'k;QPSK phase noise +/- pi/32 ldpc;')\r
+\r
+ hold off;\r
+ xlabel('Eb/N0')\r
+ ylabel('BER')\r
+ grid("minor")\r
+ axis([min(Ebvec) max(Ebvec) 1E-2 1])\r
+endfunction\r
+\r
+function phase_est_hf\r
+ sim_in = standard_init();\r
+\r
+ sim_in.Rs = 200;\r
+ sim_in.Nc = 4;\r
+\r
+ sim_in.verbose = 1;\r
+ sim_in.plot_scatter = 0;\r
+\r
+ sim_in.Esvec = 2:8; \r
+ sim_in.Ntrials = 30;\r
+\r
+ sim_in.newldpc = 1;\r
+ sim_in.ldpc_code_rate = 1/2;\r
+ sim_in.ldpc_code = 1;\r
+\r
+ sim_in.phase_est = 0;\r
+ sim_in.phase_est_method = 2;\r
+ sim_in.Np = 3;\r
+ sim_in.phase_offset = 0;\r
+ sim_in.w_offset = 0;\r
+\r
+ sim_in.hf_sim = 1;\r
+ sim_in.hf_mag_only = 1;\r
+\r
+ Ebvec = sim_in.Esvec - 10*log10(2);\r
+\r
+ baseline = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.hf_mag_only = 1;\r
+ sim_in.phase_est_method = 2;\r
+ sim_in.phase_est = 1;\r
+ sim_in.Np = 7;\r
+ pilot_7 = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.hf_mag_only = 0;\r
+ pilot_7b = ber_test(sim_in, 'qpsk');\r
+\r
+if 0\r
+ sim_in.phase_est = 0;\r
+ dqpsk = ber_test(sim_in, 'dqpsk');\r
+\r
+ sim_in.phase_est = 1;\r
+ sim_in.phase_est_method = 3;\r
+ sim_in.Np = 41;\r
+ dqpsk_strip_41 = ber_test(sim_in, 'dqpsk');\r
+end\r
+ \r
+ figure(1); \r
+ clf;\r
+ semilogy(baseline.Ebvec, baseline.BERvec,'r;QPSK CCIR poor;')\r
+ hold on;\r
+ semilogy(baseline.Ebvec, baseline.BERldpcvec,'r;QPSK CCIR poor ldpc;')\r
+ semilogy(pilot_7.Ebvec, pilot_7.BERvec,'b;QPSK CCIR poor ldpc pilot 7;')\r
+ semilogy(pilot_7.Ebvec, pilot_7.BERldpcvec,'b;QPSK CCIR poor ldpc pilot 7;')\r
+ semilogy(pilot_7b.Ebvec, pilot_7b.BERvec,'g;QPSK CCIR poor ldpc pilot 7b;')\r
+ semilogy(pilot_7b.Ebvec, pilot_7b.BERldpcvec,'g;QPSK CCIR poor ldpc pilot 7b;')\r
+if 0\r
+ semilogy(dqpsk.Ebvec, dqpsk.BERvec,'c;DQPSK CCIR poor ldpc;')\r
+ semilogy(dqpsk.Ebvec, dqpsk.BERldpcvec,'c;DQPSK CCIR poor ldpc;')\r
+ semilogy(dqpsk_strip_41.Ebvec, dqpsk_strip_41.BERvec,'m;DQPSK CCIR poor ldpc strip 41;')\r
+ semilogy(dqpsk_strip_41.Ebvec, dqpsk_strip_41.BERldpcvec,'m;DQPSK CCIR poor ldpc strip 41;')\r
+end\r
+\r
+ hold off;\r
+ xlabel('Eb/N0')\r
+ ylabel('BER')\r
+ grid("minor")\r
+ axis([min(Ebvec) max(Ebvec) 1E-2 1])\r
+\r
+endfunction\r
+\r
+function phase_est_awgn\r
+ sim_in = standard_init();\r
+\r
+ sim_in.Rs = 100;\r
+ sim_in.Nc = 8;\r
+\r
+ sim_in.verbose = 1;\r
+ sim_in.plot_scatter = 0;\r
+\r
+ sim_in.Esvec = 0:0.5:3; \r
+ sim_in.Ntrials = 30;\r
+\r
+ sim_in.newldpc = 1;\r
+ sim_in.ldpc_code_rate = 1/2;\r
+ sim_in.ldpc_code = 1;\r
+\r
+ sim_in.phase_est = 0;\r
+ sim_in.phase_est_method = 1;\r
+ sim_in.Np = 3;\r
+ sim_in.phase_offset = 0;\r
+ sim_in.w_offset = 0;\r
+\r
+ sim_in.hf_sim = 0;\r
+ sim_in.hf_mag_only = 1;\r
+\r
+ ideal = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.phase_est = 1;\r
+ sim_in.Np = 21;\r
+ sim_in.phase_est_method = 3;\r
+ strip_21_mag = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.Np = 41;\r
+ strip_41_mag = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.phase_est_method = 1;\r
+ sim_in.Np = 21;\r
+ strip_21 = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.Np = 41;\r
+ strip_41 = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.Np = 7;\r
+ sim_in.phase_est_method = 2;\r
+ pilot_7 = ber_test(sim_in, 'qpsk');\r
+\r
+ Ebvec = sim_in.Esvec - 10*log10(2);\r
+\r
+ figure(1); \r
+ clf;\r
+ semilogy(ideal.Ebvec, ideal.BERvec,'r;QPSK;')\r
+ hold on;\r
+ semilogy(ideal.Ebvec, ideal.BERldpcvec,'r;QPSK LDPC;')\r
+ semilogy(strip_21.Ebvec, strip_21.BERvec,'g;QPSK strip 21;')\r
+ semilogy(strip_21.Ebvec, strip_21.BERldpcvec,'g;QPSK LDPC strip 21;')\r
+ semilogy(strip_41.Ebvec, strip_41.BERvec,'b;QPSK strip 41;')\r
+ semilogy(strip_41.Ebvec, strip_41.BERldpcvec,'b;QPSK LDPC strip 41;')\r
+ semilogy(strip_21_mag.Ebvec, strip_21_mag.BERvec,'m;QPSK strip 21 mag;')\r
+ semilogy(strip_21_mag.Ebvec, strip_21_mag.BERldpcvec,'m;QPSK LDPC strip 21 mag;')\r
+ semilogy(strip_41_mag.Ebvec, strip_41_mag.BERvec,'c;QPSK strip 41 mag;')\r
+ semilogy(strip_41_mag.Ebvec, strip_41_mag.BERldpcvec,'c;QPSK LDPC strip 41 mag;')\r
+ semilogy(pilot_7.Ebvec, pilot_7.BERvec,'k;QPSK pilot 7;')\r
+ semilogy(pilot_7.Ebvec, pilot_7.BERldpcvec,'k;QPSK LDPC pilot 7;')\r
+\r
+ hold off;\r
+ xlabel('Eb/N0')\r
+ ylabel('BER')\r
+ grid("minor")\r
+ axis([min(Ebvec) max(Ebvec) 1E-2 1])\r
+endfunction\r
+\r
+% Start simulations ---------------------------------------\r
+\r
+more off;\r
+\r
+%ideal();\r
+phase_est_hf();\r
+%phase_est_awgn();\r
projects / freetel-svn-tracking.git / commitdiff