freq_pol(c) = 2*pi*carrier_freq/Fs;
freq(c) = exp(j*freq_pol(c));
end
-for c=Nc/2+1:Nc
+for c=floor(Nc/2)+1:Nc
%carrier_freq = (-Nc/2 + c)*Fsep + Fcentre;
carrier_freq = (-Nc/2 + c)*Fsep;
freq_pol(c) = 2*pi*carrier_freq/Fs;
for nn = 1: Ntrials \r
st = (nn-1)*code_param.symbols_per_frame + 1;\r
en = (nn)*code_param.symbols_per_frame;\r
- detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r(st:en), EsNo);\r
+ detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r(st:en), EsNo, ones(1,code_param.symbols_per_frame));\r
st = (nn-1)*code_param.data_bits_per_frame + 1;\r
en = (nn)*code_param.data_bits_per_frame;\r
error_positions = xor( detected_data(1:code_param.data_bits_per_frame), data(st:en) );\r
--- /dev/null
+% test_ml.m
+% David Rowe Oct 2014
+%
+
+% Simulation to test FDM QPSK with maximum likelihood decoding on
+% fading channels.
+
+1;
+
+% main test function
+
+function sim_out = ber_test(sim_in, modulation)
+ Fs = 8000;
+
+ verbose = sim_in.verbose;
+ Ntrials = sim_in.Ntrials;
+ Esvec = sim_in.Esvec;
+ 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_mag_only = sim_in.hf_mag_only;
+ framesize = sim_in.framesize;
+ ml = sim_in.ml;
+
+ bps = 2;
+ Nc = Nsymb = framesize/bps; % total number of symbols
+
+ prev_sym_tx = qpsk_mod([0 0])*ones(1,Nc);
+ prev_sym_rx = qpsk_mod([0 0])*ones(1,Nc);
+ r = qpsk_mod([0 0])*ones(Nc,4);
+
+ % 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;
+
+ tx_symb_log = [];
+ rx_symb_log = [];
+ errors_log = [];
+ Nerrs_log = [];
+
+ % init HF channel
+
+ hf_n = 1;
+
+ % simulation starts here-----------------------------------
+
+ for nn = 1: Ntrials
+
+ tx_bits = round(rand(1,framesize));
+
+ % modulate --------------------------------------------
+
+ for c=1:Nc
+ tx_symb(c) = qpsk_mod(tx_bits((2*(c-1)+1):(2*c)));
+ if strcmp(modulation,'dqpsk')
+ tx_symb(c) *= prev_sym_tx(c);
+ prev_sym_tx(c) = tx_symb(c);
+ end
+ end
+ s_ch = tx_symb;
+
+ % HF channel simulation ------------------------------------
+
+ if hf_sim
+
+ % separation between carriers. Note this 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
+
+ for c=1:Nc
+ ahf_model = hf_gain*(spread(hf_n) + exp(-j*c*wsep*nhfdelay)*spread_2ms(hf_n));
+ if hf_mag_only
+ s_ch(c) *= abs(ahf_model);
+ else
+ s_ch(c) *= ahf_model;
+ end
+ hf_model(hf_n, c) = ahf_model;
+ end
+ hf_n++;
+ 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,Nc) + j*randn(1,Nc));
+
+ s_ch = s_ch + noise;
+
+ % de-modulate
+
+ rx_bits = zeros(1, framesize);
+ for c=1:Nc
+
+ r(c,1:3) = r(c, 2:4);
+ r(c,4) = s_ch(c);
+
+ rx_symb(c) = s_ch(c);
+ if strcmp(modulation,'dqpsk')
+ tmp = rx_symb(c);
+ rx_symb(c) *= conj(prev_sym_rx(c)/abs(prev_sym_rx(c)));
+ prev_sym_rx(c) = tmp;
+ end
+
+ % r(c,:)
+
+ if ml == 0
+ rx_bits((2*(c-1)+1):(2*c)) = qpsk_demod(rx_symb(c));
+ else
+ tx = [1 j -j -1];
+ max_eta = 0; max_symb = tx(1);
+ for k=1:4
+ for k_1=1:4
+ for k_2=1:4
+ %eta = abs(r(c,1) + r(c,3)*tx(k)'*tx(k_1)' + r(c,2)*tx(k_1)')^2;
+ eta = abs(r(c,1) + r(c,4)*tx(k)'*tx(k_1)'*tx(k_2)' + r(c,3)*tx(k_1)'*tx(k_2)' + r(c,2)*tx(k_2)')^2;
+ %printf(" %d %d %f \n", k_1, k, eta);
+ if eta > max_eta
+ max_eta = eta;
+ max_symb = tx(k);
+ end
+ end
+ end
+ end
+ rx_bits((2*(c-1)+1):(2*c)) = qpsk_demod(max_symb);
+ end
+
+ rx_symb_log = [rx_symb_log rx_symb(c)];
+ end
+
+ % Measure BER
+
+ error_positions = xor(rx_bits, tx_bits);
+ Nerrs = sum(error_positions);
+ Terrs += Nerrs;
+ Tbits += length(tx_bits);
+ errors_log = [errors_log error_positions];
+ Nerrs_log = [Nerrs_log Nerrs];
+ end
+
+ TERvec(ne) = Terrs;
+ BERvec(ne) = Terrs/Tbits;
+
+ if verbose
+ av_tx_pwr = (tx_symb_log * tx_symb_log')/length(tx_symb_log)
+
+ printf("EsNo (dB): %3.1f Terrs: %d BER %4.3f QPSK BER theory %4.3f av_tx_pwr: %3.2f", EsNodB, Terrs,
+ Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), av_tx_pwr);
+ printf("\n");
+ end
+ end
+
+ Ebvec = Esvec - 10*log10(bps);
+ sim_out.BERvec = BERvec;
+ sim_out.Ebvec = Ebvec;
+ sim_out.TERvec = TERvec;
+ sim_out.errors_log = errors_log;
+
+ if plot_scatter
+ figure(2);
+ clf;
+ scat = rx_symb_log .* exp(j*pi/4);
+ plot(real(scat), imag(scat),'+');
+ title('Scatter plot');
+
+ if hf_sim
+ figure(3);
+ clf;
+
+ y = 1:(hf_n-1);
+ x = 1:Nc;
+ EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance);
+ EsNodBSurface(find(EsNodBSurface < -5)) = -5;
+ mesh(x,y,EsNodBSurface);
+ grid
+ axis([1 Nc 1 Rs*5 -5 15])
+ title('HF Channel Es/No');
+
+ if verbose
+ av_hf_pwr = sum(abs(hf_model(y)).^2)/(hf_n-1);
+ printf("average HF power: %3.2f over %d symbols\n", av_hf_pwr, hf_n-1);
+ end
+ end
+
+ figure(4)
+ clf
+ stem(Nerrs_log)
+ 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
+
+
+function sim_in = standard_init
+ sim_in.verbose = 1;
+ sim_in.plot_scatter = 0;
+
+ sim_in.Esvec = 5;
+ sim_in.Ntrials = 30;
+ sim_in.Rs = 50;
+ sim_in.framesize = 8;
+ sim_in.ml = 0;
+
+ sim_in.hf_delay_ms = 2;
+ sim_in.hf_sim = 0;
+ sim_in.hf_mag_only = 0;
+endfunction
+
+
+function test_curves
+
+ sim_in = standard_init();
+
+ sim_in.Ntrials = 1000;
+
+ sim_in.hf_sim = 0;
+ sim_in.plot_scatter = 0;
+ sim_in.Esvec = 5:15;
+ Ebvec = sim_in.Esvec - 10*log10(2);
+ BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10)));
+ sim_in.ml = 0;
+ sim_dqpsk = ber_test(sim_in, 'dqpsk');
+ sim_in.ml = 1;
+ sim_dqpsk_ml = ber_test(sim_in, 'dqpsk');
+ sim_in.hf_sim = 1;
+ sim_in.ml = 0;
+ sim_dqpsk_hf = ber_test(sim_in, 'dqpsk');
+ sim_in.ml = 1;
+ sim_dqpsk_ml_hf = ber_test(sim_in, 'dqpsk');
+
+ figure(1);
+ clf;
+ semilogy(Ebvec, BER_theory,'r;QPSK theory;')
+ hold on;
+ semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;')
+ semilogy(sim_dqpsk_ml.Ebvec, sim_dqpsk_ml.BERvec,'k;DQPSK ML AWGN;')
+ semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'c;DQPSK HF;')
+ semilogy(sim_dqpsk_ml_hf.Ebvec, sim_dqpsk_ml_hf.BERvec,'k;DQPSK ML HF;')
+ hold off;
+
+ xlabel('Eb/N0')
+ ylabel('BER')
+ grid("minor")
+ axis([min(Ebvec) max(Ebvec) 1E-3 1])
+endfunction
+
+
+function test_single
+ sim_in = standard_init();
+
+ sim_in.verbose = 1;
+ sim_in.plot_scatter = 1;
+ sim_in.Ntrials = 500;
+
+ sim_in.hf_mag_only = 0;
+ sim_in.hf_sim = 1;
+ sim_in.ml = 0;
+ sim_in.Esvec = 10;
+
+ sim_qpsk_hf = ber_test(sim_in, 'dqpsk');
+endfunction
+
+
+% Start simulations ---------------------------------------
+
+more off;
+
+test_curves();
+%test_single();
\r
sim_in.hf_sim = 0;\r
sim_in.plot_scatter = 0;\r
- sim_in.Esvec = 2:15; \r
+ sim_in.Esvec = 2:10; \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_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
+ sim_in.ldpc_code = 1;\r
+ sim_qpsk_hf_ldpc = ber_test(sim_in, 'qpsk');\r
+ sim_in.hf_mag_only = 0;\r
+ sim_dqpsk_hf_ldpc = ber_test(sim_in, 'dqpsk');\r
\r
figure(1); \r
clf;\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
+ semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'g;DQPSK AWGN;')\r
+ semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'r;DQPSK HF;')\r
+ semilogy(sim_qpsk_hf_ldpc.Ebvec, sim_qpsk_hf_ldpc.BERldpcvec,'b;QPSK HF LDPC 1/2;')\r
+ semilogy(sim_dqpsk_hf_ldpc.Ebvec, sim_dqpsk_hf_ldpc.BERldpcvec,'b;DQPSK HF LDPC 1/2;')\r
\r
hold off;\r
xlabel('Eb/N0')\r
\r
more off;\r
\r
-test_phase_est();\r
+ideal();\r
+++ /dev/null
-% test_ucinter.m
-% David Rowe August 2014
-%
-
-% FDM QPSK modem simulation to test uncododed interleaving ideas on HF
-% channels without building a full blown modem.
-
-% [X] baseline QPSK simulation AWGN
-% [ ] Spreading function
-% [ ] visualise different carriers with and without spreading
-% [ ] prove perf same as AWGN when one carrier is knocked out
-% + AWGN and "faded channel" with same average SNR
-% + use contrived example
-% + then try less contrived but still well behaived maths channel
-
-% Ideas:
-% + decode quickly but then slow down while playing, using full interleave
-% + exp type window so we can decode using current symbols alone in gd SNR, or extend window over greater time
-% + like root nyq filtering, combining multiple symbols, weighted
-% + spreading could make it worse, e.g. short term average might be very low
-% + SSB spreads out over a long way. We could so this too! Like spread spectrum. Don't have to be related to
-% carrier width. Minimise amount of stuff that gets nailed. But didn't we try this with "1600 wide"? OK, so
- key issue was BER. That, is what we need to impove, using "high areas".
-
-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;
- 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]);
-
- tx_bits_buf = zeros(1,2*framesize);
- rx_bits_buf = zeros(1,2*framesize);
- rx_symb_buf = zeros(1,2*Nsymb);
-
- % 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;
-
- 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) );
-
- % modulate --------------------------------------------
-
- 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
- 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];
-
- % organise into carriers to apply frequency and phase offset
-
- for i=1:Nc:Nsymb
- for k=1:Nc
- s_ch(i+k-1) = s_ch(i+k-1)*exp(j*phase_offset) + noise(i+k-1);
- end
- phase_offset += w_offset;
- end
-
- % 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);
- end
-
- TERvec(ne) = Terrs;
- BERvec(ne) = Terrs/Tbits;
-
- if verbose
- printf("EsNo (dB): %f Terrs: %d BER %f BER theory %f", EsNodB, Terrs,
- Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)));
- 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 plot_scatter
- figure(2);
- clf;
- scat = rx_symb_log .* exp(j*pi/4);
- plot(real(scat), imag(scat),'+');
- title('Scatter plot');
-
- 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 -5 15])
- title('HF Channel Es/No');
-
- figure(4);
- clf;
- %mesh(x,y,unwrap(angle(hf_model(y,:))));
- subplot(211)
- plot(y,abs(hf_model(y,1)))
- title('HF Channel Carrier 1 Mag');
- subplot(212)
- plot(y,angle(hf_model(y,1)))
- title('HF Channel Carrier 1 Phase');
-
- 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
-
-function sim_in = standard_init
- sim_in.verbose = 1;
- sim_in.plot_scatter = 0;
-
- 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.phase_noise_amp = 0;
-
- sim_in.hf_delay_ms = 2;
- sim_in.hf_sim = 0;
- sim_in.hf_phase_only = 0;
- sim_in.hf_mag_only = 1;
-endfunction
-
-function test_ideal
-
- sim_in = standard_init();
-
- sim_in.verbose = 1;
- sim_in.plot_scatter = 1;
-
- sim_in.Esvec = 5;
- sim_in.hf_sim = 1;
- sim_in.Ntrials = 30;
-
- sim_qpsk_hf = ber_test(sim_in, 'qpsk');
-
- sim_in.hf_sim = 0;
- sim_in.plot_scatter = 0;
- sim_in.Esvec = 2:5;
- 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');
-
- figure(1);
- clf;
- semilogy(Ebvec, BER_theory,'r;QPSK theory;')
- hold on;
- semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec,'g;QPSK AWGN;')
- semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;')
- hold off;
-
- xlabel('Eb/N0')
- ylabel('BER')
- grid("minor")
- axis([min(Ebvec) max(Ebvec) 1E-3 1])
-endfunction
-
-
-% Start simulations ---------------------------------------
-
-more off;
-
-test_ideal();
projects / freetel-svn-tracking.git / commitdiff