end\r
\r
error_positions = xor(rx_bits, tx_bits);\r
- sim_out.errors_log = [sim_out.errors_log error_positions];\r
Nerrs = sum(error_positions);\r
sim_out.Nerrs = [sim_out.Nerrs Nerrs];\r
Terrs += Nerrs;\r
Tbits += length(tx_bits);\r
-\r
+ \r
+ sim_out.errors_log = [sim_out.errors_log error_positions];\r
end\r
\r
TERvec(ne) = Terrs;\r
BERvec(ne) = Terrs/Tbits;\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
+ printf("EsNo (dB): %f Terrs: %d BER %f ", EsNodB, Terrs, Terrs/Tbits);\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
+ printf("Terrs: %d BER %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs,\r
+ Terrs/Tbits, 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
% Use linear mapping function in dB domain to map to symbol power\r
\r
power_map_x = [ 0 20 24 40 50 ];\r
-power_map_y = [-3 -3 0 6 6];\r
+power_map_y = [-6 -6 0 6 6];\r
mapped_pdB = interp1(power_map_x, power_map_y, pdB);\r
\r
-sim_in.symbol_amp = 10 .^ (mapped_pdB/20);\r
-%sim_in.symbol_amp = ones(1,length(pdB));\r
+%sim_in.symbol_amp = 10 .^ (mapped_pdB/20);\r
+sim_in.symbol_amp = ones(1,length(pdB));\r
sim_in.plot_scatter = 1;\r
sim_in.verbose = 2;\r
sim_in.hf_sim = 1;\r
sim_in.Esvec = 10;\r
+sim_in.Ntrials = 400;\r
\r
dqpsk_pwr_hf = ber_test(sim_in);\r
\r
hold on;\r
plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.hf_model_pwr),'g;Fading (dB);');\r
plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.snr_log),'r;Es/No (dB);');\r
+\r
+ber = dqpsk_pwr_hf.Nerrs/sim_in.framesize;\r
+ber_clip = ber;\r
+ber_clip(find(ber > 0.2)) = 0.2;\r
+plot((1:sim_in.Ntrials)*M, -20+100*ber_clip,'k;BER (0-20%);');\r
hold off;\r
+axis([1 sim_in.Ntrials*M -20 20])\r
\r
fep=fopen("dqpsk_errors_pwr.bin","wb"); fwrite(fep, dqpsk_pwr_hf.errors_log, "short"); fclose(fep);\r
-fmute=fopen("dqpsk_mute.bin","wb"); fwrite(fmute, mute, "short"); fclose(fmute);\r
+fber=fopen("ber.bin","wb"); fwrite(fber, ber, "float"); fclose(fber);\r
plot_scatter = sim_in.plot_scatter;\r
Rs = sim_in.Rs;\r
hf_sim = sim_in.hf_sim;\r
+ Nhfdelay = floor(sim_in.hf_delay_ms*Fs/1000);\r
Nc = sim_in.Nc;\r
symbol_amp = sim_in.symbol_amp;\r
\r
hrn = gen_rn_coeffs(alpha, T, Rs, Nfiltsym, M);\r
Nfilter = length(hrn);\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;\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. Transpose for convenience\r
+\r
+ spread = transpose(spread(1000:length(spread)));\r
+ spread_2ms = transpose(spread_2ms(1000: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
+ variance = Fs/(Rs*EsNo);\r
if verbose > 1\r
printf("EsNo (dB): %f EsNo: %f variance: %f\n", EsNodB, EsNo, variance);\r
end\r
sim_out.errors_log = [];\r
sim_out.tx_baseband_log = [];\r
sim_out.rx_filt_log = [];\r
+ symbol_amp_index = 1;\r
+ \r
+ % init filter memories and LOs\r
+\r
+ tx_filter_memory = zeros(Nc, Nfilter);\r
+ rx_filter_memory = zeros(Nc, Nfilter);\r
+ s_delay_line_filt = zeros(Nc, Nfiltsym);\r
+ phase_tx = ones(1,Nc);\r
+ phase_rx = ones(1,Nc);\r
+ Fcentre = 1500; Fsep = (1+alpha)*Rs;\r
+ freq = Fcentre + Fsep*((-Nc/2+0.5):(Nc/2-0.5));\r
+ freq = exp(j*freq*2*pi/Fs);\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
- symbol_amp_index = 1;\r
+ sc = 1; hf_n = 1;\r
+ hf_sim_delay_line = zeros(1,M+Nhfdelay);\r
+ freq_sample = Fcentre + ((Fsep*(-Nc/2)):50:(Fsep*(Nc/2)));\r
+ freq_sample = (2*pi/Fs)*freq_sample;\r
+ hf_model = ones(Ntrials*Nsymb/Nc, length(freq_sample)); % defaults for plotting surface\r
\r
% bunch of outputs we log for graphing\r
\r
sim_out.snr_log = [];\r
sim_out.hf_model_pwr = [];\r
sim_out.tx_fdm_log = [];\r
-\r
- % init filter memories and LOs\r
-\r
- tx_filter_memory = zeros(Nc, Nfilter);\r
- rx_filter_memory = zeros(Nc, Nfilter);\r
- s_delay_line_filt = zeros(Nc, Nfiltsym);\r
- phase_tx = ones(1,Nc);\r
- phase_rx = ones(1,Nc);\r
- Fcentre = 1500; Fsep = (1+alpha)*Rs;\r
- freq = (2*pi/Fs)*(Fcentre/2 + Fsep*(1:Nc));\r
- freq = exp(j*freq);\r
+ C_log = [];\r
\r
for nn = 1: Ntrials\r
\r
symbol_amp_index++;\r
s_ch = s;\r
\r
+ % Now we start processing frame Nc symbols at a time to model parallel carriers\r
+\r
for i=1:Nc:Nsymb\r
\r
- % Delay tx symbols to match delay due to filters qpsk\r
+ % Delay tx symbols to match delay due to filters. qpsk\r
% (rather than dqpsk) symbols used for convenience as\r
- % it's easy to shift symbols than pair of bits\r
+ % it's easy to shift symbols than pairs of bits\r
\r
s_delay_line_filt(:,1:Nfiltsym-1) = s_delay_line_filt(:,2:Nfiltsym);\r
s_delay_line_filt(:,Nfiltsym) = s_qpsk(i:i+Nc-1);\r
\r
sim_out.tx_baseband_log = [sim_out.tx_baseband_log tx_baseband];\r
\r
- % upcovert\r
+ % upconvert\r
\r
tx_fdm = zeros(1,M);\r
\r
\r
% HF channel\r
\r
- rx_fdm = tx_fdm;\r
+ if hf_sim\r
+ hf_sim_delay_line(1:Nhfdelay) = hf_sim_delay_line(M+1:M+Nhfdelay);\r
+ hf_sim_delay_line(Nhfdelay+1:M+Nhfdelay) = tx_fdm;\r
+\r
+ tx_fdm = tx_fdm.*spread(sc:sc+M-1) + hf_sim_delay_line(1:M).*spread_2ms(sc:sc+M-1);\r
+ tx_fdm *= hf_gain;\r
+\r
+ % sample HF channel spectrum in middle of this symbol for plotting\r
+\r
+ hf_model(hf_n,:) = hf_gain*(spread(sc+M/2) + exp(-j*freq_sample*Nhfdelay)*spread_2ms(sc+M/2));\r
+\r
+ sc += M;\r
+ hf_n++;\r
+ end\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,M) + j*randn(1,M));\r
+ noise_log = [noise_log noise];\r
+\r
+ % apply frequency and phase offset and noise\r
+\r
+ for k=1:M\r
+ rx_fdm(k) = tx_fdm(k)*exp(j*phase_offset) + noise(k);\r
+ phase_offset += w_offset;\r
+ end\r
\r
% downconvert\r
\r
s_ch(i:i+Nc-1) = rx_filt;\r
end\r
\r
- % HF channel simulation ------------------------------------\r
- \r
+ % est HF model power for entire code frame (which could be several symbols)\r
+\r
if hf_sim\r
+ frame_hf_model = reshape(hf_model(hf_n-Nsymb/Nc:hf_n-1,:),1,(Nsymb/Nc)*length(freq_sample)); \r
+ sim_out.hf_model_pwr = [sim_out.hf_model_pwr mean(abs(frame_hf_model).^2)];\r
+ else \r
+ sim_out.hf_model_pwr = [sim_out.hf_model_pwr 1];\r
end\r
\r
% "genie" SNR estimate \r
\r
- snr = (s_ch*s_ch')/(Nsymb*variance);\r
+ snr = (tx_fdm*tx_fdm')/(M*variance);\r
sim_out.snr_log = [sim_out.snr_log snr];\r
- sim_out.hf_model_pwr = [sim_out.hf_model_pwr mean(hf_fading.^2)];\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
- \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) + noise(i+k-1);\r
- end \r
- phase_offset += w_offset;\r
- end\r
- \r
+ \r
% de-modulate\r
\r
rx_bits = zeros(1, framesize);\r
BERvec(ne) = Terrs/Tbits;\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
+ printf("EsNo (dB): %f Terrs: %d BER %f ", EsNodB, Terrs, Terrs/Tbits);\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
+ printf("Terrs: %d BER %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs,\r
+ Terrs/Tbits, var(sim_out.tx_fdm_log), var(noise_log),\r
+ var(sim_out.tx_fdm_log)/(Nc*Rs), var(noise_log)/Fs, (var(sim_out.tx_fdm_log)/(Nc*Rs))/(var(noise_log)/Fs));\r
end\r
end\r
\r
figure(3);\r
clf; \r
y = 1:Rs*2;\r
- x = 1:Nc;\r
- EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance);\r
+ [rr,cc] = size(hf_model);\r
+ x = 1:cc;\r
+ EsNodBSurface = 20*log10(abs(hf_model(y,:))) + EsNodB;\r
mesh(x,y,EsNodBSurface);\r
grid\r
title('HF Channel Es/No');\r
power_map_y = [-6 -6 0 6 6];\r
mapped_pdB = interp1(power_map_x, power_map_y, pdB);\r
\r
-%sim_in.symbol_amp = 10 .^ (mapped_pdB/20);\r
-sim_in.symbol_amp = ones(1,length(pdB));\r
+sim_in.symbol_amp = 10 .^ (mapped_pdB/20);\r
+%sim_in.symbol_amp = ones(1,length(pdB));\r
sim_in.plot_scatter = 1;\r
sim_in.verbose = 2;\r
+sim_in.hf_delay_ms = 2;\r
sim_in.hf_sim = 1;\r
-sim_in.Esvec = 100;\r
+sim_in.Esvec = 10;\r
sim_in.Ntrials = 100;\r
\r
dqpsk_pwr_hf = ber_test(sim_in);\r
\r
figure(5)\r
clf;\r
-\r
s = load_raw("../raw/ve9qrp.raw");\r
M=320; M_on_2 = M/2; % processing delay between input speech and centre of analysis window\r
+subplot(211)\r
plot(M_on_2:(M_on_2-1+sim_in.Ntrials*M),s(1:sim_in.Ntrials*M))\r
hold on;\r
plot((1:sim_in.Ntrials)*M, 5000*sim_in.symbol_amp(1:sim_in.Ntrials),'r');\r
hold off;\r
axis([1 sim_in.Ntrials*M -3E4 3E4]);\r
+subplot(212)\r
+plot(real(dqpsk_pwr_hf.tx_fdm_log));\r
+\r
\r
figure(6)\r
clf;\r
projects / freetel-svn-tracking.git / commitdiff