From: drowe67 Date: Fri, 27 Apr 2018 21:19:02 +0000 (+0000) Subject: Es/No and SNR estimation, works well on AWGN of slow fading, afew dB low on high... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=0b49760a497db62b4dd10f8861733eb9c68ded87;p=freetel-svn-tracking.git Es/No and SNR estimation, works well on AWGN of slow fading, afew dB low on high SNR fast (1Hz) fading git-svn-id: https://svn.code.sf.net/p/freetel/code@3528 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/ofdm_dev.m b/codec2-dev/octave/ofdm_dev.m index 15f2f88b..68b75694 100644 --- a/codec2-dev/octave/ofdm_dev.m +++ b/codec2-dev/octave/ofdm_dev.m @@ -273,6 +273,8 @@ function [sim_out rx states] = run_sim(sim_in) snrdB = 10*log10(var(rx)/var(noise)) + 10*log10(8000) - 10*log10(3000); rx += noise; + rx *= sim_in.gain; + % some spare samples at end to avoid overflow as est windows may poke into the future a bit rx = [rx zeros(1,Nsamperframe)]; @@ -285,6 +287,7 @@ function [sim_out rx states] = run_sim(sim_in) foff_est_hz_log = []; Nerrs_log = []; Nerrs_coded_log = []; rx_bits = []; rx_np = []; rx_amp = []; + sig_var_log = []; noise_var_log = []; % reset some states for each EbNo simulation point @@ -317,13 +320,20 @@ function [sim_out rx states] = run_sim(sim_in) [arx_bits states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in); + %mean_amp = mean(mean(arx_amp)); + %mean_amp = 1/1000; + %printf("mean amp: %f\n", mean_amp); + %arx_np /= mean_amp; arx_amp /= mean_amp; + rx_bits = [rx_bits arx_bits]; rx_np = [rx_np arx_np]; rx_amp = [rx_amp arx_amp]; timing_est_log = [timing_est_log states.timing_est]; delta_t_log = [delta_t_log states.delta_t]; foff_est_hz_log = [foff_est_hz_log states.foff_est_hz]; phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log]; + sig_var_log = [sig_var_log states.sig_var]; + noise_var_log = [noise_var_log states.noise_var]; + %printf("sig_var: %f noise_var: %f\n", states.sig_var, states.noise_var); end - assert(length(rx_bits) == Nbits); % Optional de-interleave on rx QPSK symbols @@ -370,7 +380,7 @@ function [sim_out rx states] = run_sim(sim_in) st = (f-1)*Nbitsperframe/bps + 1; en = st + Nbitsperframe/bps - 1; r = rx_np(st:en); fade = rx_amp(st:en); - fade + % optional LDPC decode if ldpc_en @@ -434,7 +444,11 @@ function [sim_out rx states] = run_sim(sim_in) printf("Raw EbNodB..: % -4.1f BER: %5.4f Tbits: %5d Terrs: %5d PER: %5.4f Tpackets: %5d Tpacket_errs: %5d\n", EbNodB_raw, Terrs/Nbits, Nbits, Terrs, Tpacketerrs/Tpackets, Tpackets, Tpacketerrs); - + EsNo = mean(sig_var_log)/mean(noise_var_log); + %printf("Es/No est dB: % -4.1f\n", 10*log10(EsNo)); + sim_out.snrdB(nn) = snrdB; + sim_out.snr_estdB(nn) = 10*log10(EsNo) + 10*log10(Nc*Rs/3000); + % returns results for plotting curves if ldpc_en || diversity_en @@ -481,7 +495,7 @@ function [sim_out rx states] = run_sim(sim_in) subplot(212) stem(Nerrs_coded_log/(Nbitsperframe*rate)); title("Coded BER/frame"); - else + else title("BER/frame"); stem(Nerrs_log/Nbitsperframe); end @@ -494,7 +508,18 @@ function [sim_out rx states] = run_sim(sim_in) mx = max(Tx_dB); axis([0 Fs/2 mx-60 mx]) - + figure(7); clf; + plot(10*log10(sig_var_log),'b;Es;'); + hold on; + plot(10*log10(noise_var_log),'r;No;'); + snr_estdB = 10*log10(sig_var_log) - 10*log10(noise_var_log) + 10*log10(Nc*Rs/3000); + snr_est_smoothed_dB = filter(0.1,[1 -0.9],snr_estdB); + plot(snr_estdB,'g;SNR3k;'); + plot(snr_est_smoothed_dB,'c;SNR3k smooth;'); + + hold off; + title('Signal and Noise Power estimates'); + #{ if hf_en figure(4); clf; @@ -513,8 +538,6 @@ function [sim_out rx states] = run_sim(sim_in) plot(angle(hf_model(:,2:Nc+1))); end #} - - end end @@ -527,18 +550,19 @@ function run_single(EbNodB = 100, error_pattern_filename); sim_in.Rs = 1/Ts; sim_in.bps = 2; sim_in.Nc = 16; sim_in.Ns = 8; sim_in.Nsec = (sim_in.Ns+1)/sim_in.Rs; % one frame, make sure sim_in.interleave_frames = 1 - %sim_in.Nsec = 1; + sim_in.Nsec = 20; sim_in.EbNodB = EbNodB; - sim_in.verbose = 2; + sim_in.verbose = 1; sim_in.hf_en = 0; sim_in.foff_hz = 0; sim_in.dfoff_hz_per_sec = 0.00; sim_in.sample_clock_offset_ppm = 0; - + sim_in.gain = 1; + sim_in.timing_en = 0; sim_in.foff_est_en = 0; - sim_in.phase_est_en = 0; + sim_in.phase_est_en = 1; load HRA_112_112.txt sim_in.ldpc_code = HRA_112_112; @@ -580,6 +604,8 @@ function run_curves sim_in.verbose = 0; sim_in.foff_hz = 0; + sim_in.gain = 1; + sim_in.timing_en = 1; sim_in.foff_est_en = 1; sim_in.phase_est_en = 1; @@ -843,6 +869,58 @@ function run_curves_estimators end +% Plot SNR actual and estimated for various channels +% Note no acquistion, as this can upset results, e.g. if +% sync is lost. We average SNR over entire run, in practice +% there will be some sort of IIR averager. + +function run_curves_snr + + Nsec_awgn = 30; + Nsec_hf = 60; + + % waveform + + Ts = 0.018; sim_in.Tcp = 0.002; + sim_in.Rs = 1/Ts; sim_in.bps = 2; sim_in.Nc = 16; sim_in.Ns = 8; + + % simulation parameters + + sim_in.verbose = 1; sim_in.foff_hz = 0; sim_in.ldpc_en = 0; sim_in.gain = 1; + sim_in.phase_est_en = sim_in.timing_en = sim_in.foff_est_en = 1; + + % AWGN simulation + + sim_in.EbNodB = 0:3:12; + sim_in.hf_en = 0; sim_in.Nsec = Nsec_awgn; + + awgn = run_sim(sim_in); + + % HF simulations + + sim_in.hf_en = 1; sim_in.Nsec = Nsec_hf; sim_in.dopplerSpreadHz = 1; hf_fast = run_sim(sim_in); + sim_in.dopplerSpreadHz = 0.5; hf_mid = run_sim(sim_in); + sim_in.dopplerSpreadHz = 0.2; hf_slow = run_sim(sim_in); + + figure(1); clf; + plot(sim_in.EbNodB, awgn.snrdB,'b+-;AWGN SNR;'); + hold on; + plot(sim_in.EbNodB, awgn.snr_estdB ,'g+-;AWGN SNR est;'); + plot(sim_in.EbNodB, hf_slow.snr_estdB,'k+-;HF 0.2 Hz SNR est;'); + plot(sim_in.EbNodB, hf_mid.snr_estdB ,'r+-;HF 0.5 Hz SNR est;'); + plot(sim_in.EbNodB, hf_fast.snr_estdB,'c+-;HF 1.0 Hz SNR est;'); + hold off; + xlabel('Eb/No (dB)'); + ylabel('SNR (dB)'); + grid; + legend('boxoff'); + title('SNR Actual and Estimated'); + legend("location", "northwest"); + print('-deps', '-color', "ofdm_dev_snr.eps") + +end + + % Run an acquisition test, returning vectors of estimation errors function [delta_ct delta_foff timing_mx_log] = acquisition_test(Ntests=10, EbNodB=100, foff_hz=0, hf_en=0, fine_en=0) @@ -1084,9 +1162,10 @@ more off; init_cml('~/cml/'); -run_single +run_single(10); %run_curves %run_curves_estimators %acquisition_histograms(0, 0) %acquisition_test(Ntests=3, EbNodB=100, foff_hz=0) %sync_metrics('freq') +%run_curves_snr diff --git a/codec2-dev/octave/ofdm_lib.m b/codec2-dev/octave/ofdm_lib.m index b07d4d10..2d17f10f 100644 --- a/codec2-dev/octave/ofdm_lib.m +++ b/codec2-dev/octave/ofdm_lib.m @@ -220,6 +220,11 @@ function states = ofdm_init(bps, Rs, Tcp, Ns, Nc) states.rx_sym = zeros(1+Ns+1+1, Nc+2); + % Es/No (SNR) est states + + states.noise_var = 0; + states.sig_var = 0; + endfunction @@ -457,9 +462,12 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, aphase_est_pilot_rect += sum(rx_sym(1,cr)*pilots(cr)'); aphase_est_pilot_rect += sum(rx_sym(2+Ns+1,cr)*pilots(cr)'); - + aphase_est_pilot(c) = angle(aphase_est_pilot_rect); - aamp_est_pilot(c) = abs(aphase_est_pilot_rect/6); % amplitude is estimated over 6 rows of pilots + + % amplitude is estimated over 12 pilot symbols, so find average + + aamp_est_pilot(c) = abs(aphase_est_pilot_rect/12); end % correct phase offset using phase estimate, and demodulate @@ -506,6 +514,39 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, end end + % estimates of signal and noise power (see cohpsk.m for further explanation) + % signal power is distance from axis on complex plane + % we just measure noise power on imag axis, as it isn't affected by fading + % using all symbols in frame worked better than just pilots + + x = sum(abs(rx_np) .^ 2)/length(rx_np); + sig_var = x; + sig_rms = sqrt(x); + + sum_x = 0; + sum_xx = 0; + n = 0; + for i=1:length(rx_np) + s = rx_np(i); + if abs(real(s)) > sig_rms + % select two constellation points on real axis + sum_x += imag(s); + sum_xx += imag(s)*imag(s); + n++; + end + end + + noise_var = 0; + if n > 1 + noise_var = (n*sum_xx - sum_x*sum_x)/(n*(n-1)); + end + + % Total noise power is twice estimate of imaginary-axis noise. This + % effectively gives us the an estimate of Es/No + + states.noise_var = 2*noise_var; + states.sig_var = sig_var; + states.rx_sym = rx_sym; states.rxbuf = rxbuf; states.nin = nin; diff --git a/codec2-dev/octave/ofdm_rx.m b/codec2-dev/octave/ofdm_rx.m index 85458f93..98412551 100644 --- a/codec2-dev/octave/ofdm_rx.m +++ b/codec2-dev/octave/ofdm_rx.m @@ -30,7 +30,7 @@ function ofdm_rx(filename, error_pattern_filename) % init logs and BER stats rx_bits = []; rx_np_log = []; timing_est_log = []; delta_t_log = []; foff_est_hz_log = []; - phase_est_pilot_log = []; + phase_est_pilot_log = []; sig_var_log = []; noise_var_log = []; Terrs = Tbits = Terrs_coded = Tbits_coded = Tpackets = Tpacketerrs = frame_count = 0; Nbitspervocframe = 28; Nerrs_coded_log = Nerrs_log = []; @@ -43,7 +43,7 @@ function ofdm_rx(filename, error_pattern_filename) %states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx(prx:nin); %prx += nin; - states.verbose = 1; + states.verbose = 0; Nerrs = 0; rx_uw = zeros(1,states.Nuwbits); @@ -82,6 +82,8 @@ function ofdm_rx(filename, error_pattern_filename) delta_t_log = [delta_t_log states.delta_t]; foff_est_hz_log = [foff_est_hz_log states.foff_est_hz]; phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log]; + sig_var_log = [sig_var_log states.sig_var]; + noise_var_log = [noise_var_log states.noise_var]; % measure uncoded bit errors on modem frame @@ -104,6 +106,8 @@ function ofdm_rx(filename, error_pattern_filename) if states.sync_start Nerrs_log = []; Terrs = Tbits = frame_count = 0; + rx_np_log = []; + sig_var_log = []; noise_var_log = []; end end @@ -117,6 +121,12 @@ function ofdm_rx(filename, error_pattern_filename) Terrs -= sum(Nerrs_log(1:Ndiscard)); Tbits -= Ndiscard*Nbitsperframe; printf("BER2.: %5.4f Tbits: %5d Terrs: %5d\n", Terrs/Tbits, Tbits, Terrs); end + + %EsNo_est = mean(sig_var_log(floor(end/2):end))/mean(noise_var_log(floor(end/2):end)); + EsNo_est = mean(sig_var_log)/mean(noise_var_log); + EsNo_estdB = 10*log10(EsNo_est); + SNR_estdB = EsNo_estdB + 10*log10(Nc*Rs/3000); + printf("Es/No est dB: % -4.1f SNR3k: %3.2f %f %f\n", EsNo_estdB, SNR_estdB, mean(sig_var_log), mean(noise_var_log)); figure(1); clf; plot(rx_np_log,'+'); @@ -149,6 +159,16 @@ function ofdm_rx(filename, error_pattern_filename) title('Errors/modem frame') axis([1 length(Nerrs_log) 0 Nbitsperframe*rate/2]); + figure(6); clf; + plot(10*log10(sig_var_log),'b;Es;'); + hold on; + plot(10*log10(noise_var_log),'r;No;'); + snr_estdB = 10*log10(sig_var_log) - 10*log10(noise_var_log) + 10*log10(Nc*Rs/3000); + snr_smoothed_estdB = filter(0.1,[1 -0.9],snr_estdB); + plot(snr_smoothed_estdB,'g;SNR3k;'); + hold off; + title('Signal and Noise Power estimates'); + if nargin == 2 fep = fopen(error_pattern_filename, "wb"); fwrite(fep, error_positions, "short");