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
+ phase_est_method = sim_in.phase_est_method;\r
+ if phase_est_method == 1\r
+ Nps = sim_in.Np; \r
+ else\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
+ 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
+\r
% Init LDPC --------------------------------------------------------------------\r
\r
if ldpc_code\r
tx_symb_log = [];\r
rx_symb_log = [];\r
noise_log = [];\r
+ mod_strip_log = [];\r
\r
% init HF channel\r
\r
for nn = 1: Ntrials\r
\r
tx_bits = round( rand( 1, framesize*rate ) );\r
-\r
+ \r
% modulate --------------------------------------------\r
\r
if ldpc_code\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
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
+ 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
+ %tx_bits(2*(i+k-1-1)+1:2*(i+k-1)) = qpsk_demod(s_delay_line(k,centre));\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
+ else\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
+ 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
rx_symb_log = [rx_symb_log rx_symb];\r
end\r
\r
- % Measure BER\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
+\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
- error_positions = xor(rx_bits, tx_bits);\r
- Nerrs = sum(error_positions);\r
- Terrs += Nerrs;\r
- Tbits += length(tx_bits);\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
+ % 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
+ if ldpc_code\r
+ detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symb_buf(st_rx_symb:en_rx_symb), min(100,EsNo), hf_fading);\r
+ %for m=1:20\r
+ % printf("%f ", qpsk_demod(rx_symb_buf(m)));\r
+ %end\r
+ %detected_data(1:19)\r
+ error_positions = xor( detected_data(1:framesize*rate), tx_bits_buf(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
figure(2);\r
clf;\r
scat = rx_symb_log .* exp(j*pi/4);\r
- plot(real(scat), imag(scat),'+');\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
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
-\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
tmp = [tmp abs(hf_model(i,:))];\r
end\r
hist(tmp);\r
+end\r
end\r
\r
endfunction\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 = 30;\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_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 = 30;\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 test_phase_est\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 = 1;\r
+\r
+ sim_in.Esvec = 5; \r
+ sim_in.Ntrials = 30;\r
+\r
+ sim_in.ldpc_code_rate = 1/2;\r
+ sim_in.ldpc_code = 1;\r
+\r
+ sim_in.phase_est = 1;\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 = 0;\r
+\r
+ tmp = ber_test(sim_in, 'qpsk');\r
+\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
+\r
+test_phase_est();\r
projects / freetel-svn-tracking.git / commitdiff