From 9cdd8f4554ec96a0b8e839efe1ffafe67d29fea7 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sun, 6 Apr 2014 09:29:30 +0000 Subject: [PATCH] in the process of adding HF model to test_dqpsk2 git-svn-id: https://svn.code.sf.net/p/freetel/code@1489 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/test_dqpsk.m | 26 ++++-- codec2-dev/octave/test_dqpsk2.m | 161 +++++++++++++++++++++----------- 2 files changed, 125 insertions(+), 62 deletions(-) diff --git a/codec2-dev/octave/test_dqpsk.m b/codec2-dev/octave/test_dqpsk.m index 3d33bc37..76beaf9c 100644 --- a/codec2-dev/octave/test_dqpsk.m +++ b/codec2-dev/octave/test_dqpsk.m @@ -199,25 +199,24 @@ function sim_out = ber_test(sim_in) end error_positions = xor(rx_bits, tx_bits); - sim_out.errors_log = [sim_out.errors_log error_positions]; Nerrs = sum(error_positions); sim_out.Nerrs = [sim_out.Nerrs Nerrs]; Terrs += Nerrs; Tbits += length(tx_bits); - + + sim_out.errors_log = [sim_out.errors_log error_positions]; 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("EsNo (dB): %f Terrs: %d BER %f ", EsNodB, Terrs, Terrs/Tbits); 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), + printf("Terrs: %d BER %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs, + Terrs/Tbits, var(tx_symb_log), var(noise_log), var(tx_symb_log), var(noise_log), var(tx_symb_log)/var(noise_log)); end end @@ -344,15 +343,16 @@ pdB = pdB(4:4:length(pdB)); % Use linear mapping function in dB domain to map to symbol power power_map_x = [ 0 20 24 40 50 ]; -power_map_y = [-3 -3 0 6 6]; +power_map_y = [-6 -6 0 6 6]; mapped_pdB = interp1(power_map_x, power_map_y, pdB); -sim_in.symbol_amp = 10 .^ (mapped_pdB/20); -%sim_in.symbol_amp = ones(1,length(pdB)); +%sim_in.symbol_amp = 10 .^ (mapped_pdB/20); +sim_in.symbol_amp = ones(1,length(pdB)); sim_in.plot_scatter = 1; sim_in.verbose = 2; sim_in.hf_sim = 1; sim_in.Esvec = 10; +sim_in.Ntrials = 400; dqpsk_pwr_hf = ber_test(sim_in); @@ -382,7 +382,13 @@ plot((1:sim_in.Ntrials)*M, 20*log10(sim_in.symbol_amp(1:sim_in.Ntrials)),'b;Es ( hold on; plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.hf_model_pwr),'g;Fading (dB);'); plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.snr_log),'r;Es/No (dB);'); + +ber = dqpsk_pwr_hf.Nerrs/sim_in.framesize; +ber_clip = ber; +ber_clip(find(ber > 0.2)) = 0.2; +plot((1:sim_in.Ntrials)*M, -20+100*ber_clip,'k;BER (0-20%);'); hold off; +axis([1 sim_in.Ntrials*M -20 20]) fep=fopen("dqpsk_errors_pwr.bin","wb"); fwrite(fep, dqpsk_pwr_hf.errors_log, "short"); fclose(fep); -fmute=fopen("dqpsk_mute.bin","wb"); fwrite(fmute, mute, "short"); fclose(fmute); +fber=fopen("ber.bin","wb"); fwrite(fber, ber, "float"); fclose(fber); diff --git a/codec2-dev/octave/test_dqpsk2.m b/codec2-dev/octave/test_dqpsk2.m index a1806b6b..67f699f7 100644 --- a/codec2-dev/octave/test_dqpsk2.m +++ b/codec2-dev/octave/test_dqpsk2.m @@ -21,6 +21,7 @@ function sim_out = ber_test(sim_in) plot_scatter = sim_in.plot_scatter; Rs = sim_in.Rs; hf_sim = sim_in.hf_sim; + Nhfdelay = floor(sim_in.hf_delay_ms*Fs/1000); Nc = sim_in.Nc; symbol_amp = sim_in.symbol_amp; @@ -41,13 +42,47 @@ function sim_out = ber_test(sim_in) hrn = gen_rn_coeffs(alpha, T, Rs, Nfiltsym, M); Nfilter = length(hrn); + % 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; + 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. Transpose for convenience + + spread = transpose(spread(1000:length(spread))); + spread_2ms = transpose(spread_2ms(1000: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; + variance = Fs/(Rs*EsNo); if verbose > 1 printf("EsNo (dB): %f EsNo: %f variance: %f\n", EsNodB, EsNo, variance); end @@ -60,14 +95,26 @@ function sim_out = ber_test(sim_in) sim_out.errors_log = []; sim_out.tx_baseband_log = []; sim_out.rx_filt_log = []; + symbol_amp_index = 1; + + % init filter memories and LOs + + tx_filter_memory = zeros(Nc, Nfilter); + rx_filter_memory = zeros(Nc, Nfilter); + s_delay_line_filt = zeros(Nc, Nfiltsym); + phase_tx = ones(1,Nc); + phase_rx = ones(1,Nc); + Fcentre = 1500; Fsep = (1+alpha)*Rs; + freq = Fcentre + Fsep*((-Nc/2+0.5):(Nc/2-0.5)); + freq = exp(j*freq*2*pi/Fs); % 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 - symbol_amp_index = 1; + sc = 1; hf_n = 1; + hf_sim_delay_line = zeros(1,M+Nhfdelay); + freq_sample = Fcentre + ((Fsep*(-Nc/2)):50:(Fsep*(Nc/2))); + freq_sample = (2*pi/Fs)*freq_sample; + hf_model = ones(Ntrials*Nsymb/Nc, length(freq_sample)); % defaults for plotting surface % bunch of outputs we log for graphing @@ -76,17 +123,7 @@ function sim_out = ber_test(sim_in) sim_out.snr_log = []; sim_out.hf_model_pwr = []; sim_out.tx_fdm_log = []; - - % init filter memories and LOs - - tx_filter_memory = zeros(Nc, Nfilter); - rx_filter_memory = zeros(Nc, Nfilter); - s_delay_line_filt = zeros(Nc, Nfiltsym); - phase_tx = ones(1,Nc); - phase_rx = ones(1,Nc); - Fcentre = 1500; Fsep = (1+alpha)*Rs; - freq = (2*pi/Fs)*(Fcentre/2 + Fsep*(1:Nc)); - freq = exp(j*freq); + C_log = []; for nn = 1: Ntrials @@ -107,11 +144,13 @@ function sim_out = ber_test(sim_in) symbol_amp_index++; s_ch = s; + % Now we start processing frame Nc symbols at a time to model parallel carriers + for i=1:Nc:Nsymb - % Delay tx symbols to match delay due to filters qpsk + % Delay tx symbols to match delay due to filters. qpsk % (rather than dqpsk) symbols used for convenience as - % it's easy to shift symbols than pair of bits + % it's easy to shift symbols than pairs of bits s_delay_line_filt(:,1:Nfiltsym-1) = s_delay_line_filt(:,2:Nfiltsym); s_delay_line_filt(:,Nfiltsym) = s_qpsk(i:i+Nc-1); @@ -137,7 +176,7 @@ function sim_out = ber_test(sim_in) sim_out.tx_baseband_log = [sim_out.tx_baseband_log tx_baseband]; - % upcovert + % upconvert tx_fdm = zeros(1,M); @@ -152,7 +191,33 @@ function sim_out = ber_test(sim_in) % HF channel - rx_fdm = tx_fdm; + if hf_sim + hf_sim_delay_line(1:Nhfdelay) = hf_sim_delay_line(M+1:M+Nhfdelay); + hf_sim_delay_line(Nhfdelay+1:M+Nhfdelay) = tx_fdm; + + tx_fdm = tx_fdm.*spread(sc:sc+M-1) + hf_sim_delay_line(1:M).*spread_2ms(sc:sc+M-1); + tx_fdm *= hf_gain; + + % sample HF channel spectrum in middle of this symbol for plotting + + hf_model(hf_n,:) = hf_gain*(spread(sc+M/2) + exp(-j*freq_sample*Nhfdelay)*spread_2ms(sc+M/2)); + + sc += M; + hf_n++; + end + + % 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,M) + j*randn(1,M)); + noise_log = [noise_log noise]; + + % apply frequency and phase offset and noise + + for k=1:M + rx_fdm(k) = tx_fdm(k)*exp(j*phase_offset) + noise(k); + phase_offset += w_offset; + end % downconvert @@ -174,32 +239,20 @@ function sim_out = ber_test(sim_in) s_ch(i:i+Nc-1) = rx_filt; end - % HF channel simulation ------------------------------------ - + % est HF model power for entire code frame (which could be several symbols) + if hf_sim + frame_hf_model = reshape(hf_model(hf_n-Nsymb/Nc:hf_n-1,:),1,(Nsymb/Nc)*length(freq_sample)); + sim_out.hf_model_pwr = [sim_out.hf_model_pwr mean(abs(frame_hf_model).^2)]; + else + sim_out.hf_model_pwr = [sim_out.hf_model_pwr 1]; end % "genie" SNR estimate - snr = (s_ch*s_ch')/(Nsymb*variance); + snr = (tx_fdm*tx_fdm')/(M*variance); sim_out.snr_log = [sim_out.snr_log snr]; - sim_out.hf_model_pwr = [sim_out.hf_model_pwr mean(hf_fading.^2)]; - - % 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); @@ -232,14 +285,13 @@ function sim_out = ber_test(sim_in) 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("EsNo (dB): %f Terrs: %d BER %f ", EsNodB, Terrs, Terrs/Tbits); 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)); + printf("Terrs: %d BER %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs, + Terrs/Tbits, var(sim_out.tx_fdm_log), var(noise_log), + var(sim_out.tx_fdm_log)/(Nc*Rs), var(noise_log)/Fs, (var(sim_out.tx_fdm_log)/(Nc*Rs))/(var(noise_log)/Fs)); end end @@ -259,8 +311,9 @@ function sim_out = ber_test(sim_in) figure(3); clf; y = 1:Rs*2; - x = 1:Nc; - EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance); + [rr,cc] = size(hf_model); + x = 1:cc; + EsNodBSurface = 20*log10(abs(hf_model(y,:))) + EsNodB; mesh(x,y,EsNodBSurface); grid title('HF Channel Es/No'); @@ -357,12 +410,13 @@ power_map_x = [ 0 20 24 40 50 ]; power_map_y = [-6 -6 0 6 6]; mapped_pdB = interp1(power_map_x, power_map_y, pdB); -%sim_in.symbol_amp = 10 .^ (mapped_pdB/20); -sim_in.symbol_amp = ones(1,length(pdB)); +sim_in.symbol_amp = 10 .^ (mapped_pdB/20); +%sim_in.symbol_amp = ones(1,length(pdB)); sim_in.plot_scatter = 1; sim_in.verbose = 2; +sim_in.hf_delay_ms = 2; sim_in.hf_sim = 1; -sim_in.Esvec = 100; +sim_in.Esvec = 10; sim_in.Ntrials = 100; dqpsk_pwr_hf = ber_test(sim_in); @@ -378,14 +432,17 @@ hold off; figure(5) clf; - s = load_raw("../raw/ve9qrp.raw"); M=320; M_on_2 = M/2; % processing delay between input speech and centre of analysis window +subplot(211) plot(M_on_2:(M_on_2-1+sim_in.Ntrials*M),s(1:sim_in.Ntrials*M)) hold on; plot((1:sim_in.Ntrials)*M, 5000*sim_in.symbol_amp(1:sim_in.Ntrials),'r'); hold off; axis([1 sim_in.Ntrials*M -3E4 3E4]); +subplot(212) +plot(real(dqpsk_pwr_hf.tx_fdm_log)); + figure(6) clf; -- 2.25.1