Initial commitHEADmaster

* codec2 cut-down version 1.2.0
* Remove codebook and generation of sources
* remove c2dec c2enc binaries
* prepare for emscripten

1 files changed, 147 insertions, 0 deletions

diff --git a/octave/esno_est.m b/octave/esno_est.m
new file mode 100644
index 0000000..590208c
--- /dev/null
+++ b/octave/esno_est.m
@@ -0,0 +1,147 @@
+% esno_est.m
+% David Rowe Mar 2017
+%
+% Functions for estimating Es/No from QPSK symbols, and supporting tests
+
+1;
+
+#{
+  ----------------------------------------------------------------------------
+  Estimate the energy and noise of received symbols.
+  
+  Signal power is distance from axis on complex
+  plane. Noise power is the distance orthogonal to each symbol, to provide an
+  estimate that is insensitive to fading that moves symbol towards he origin.
+  
+  For QAM we need to use pilots as they don't have modulation that affects 
+  estimate, for QPSK Modes we can use all rx symbols.
+  ----------------------------------------------------------------------------
+#}
+
+function EsNodB = esno_est_calc(rx_syms)
+  sig_var = sum(abs(rx_syms) .^ 2)/length(rx_syms);
+  sig_rms = sqrt(sig_var);
+   
+  sum_x = 0;
+  sum_xx = 0;
+  n = 0;
+  for i=1:length(rx_syms)
+    s = rx_syms(i);
+    % only consider symbols a reasonable distance from the origin, as these are
+    % more likely to be valid and not errors that will mess up the estimate
+    if abs(s) > sig_rms
+        % rough demodulation, determine if symbol is on real or imag axis
+        if abs(real(s)) > abs(imag(s))
+          % assume noise is orthogonal to real axis
+          sum_x  += imag(s);
+          sum_xx += imag(s)*imag(s);
+        else
+          % assume noise is orthogonal to imag axis
+          sum_x  += real(s);
+          sum_xx += real(s)*real(s);
+        end
+        n++;
+    end
+  end
+
+  % trap corner case
+  if n > 1
+    noise_var = (n*sum_xx - sum_x*sum_x)/(n*(n-1));
+  else
+    noise_var = sig_var;
+  end
+
+  % Total noise power is twice estimate of single-axis noise.
+  noise_var = 2*noise_var;
+  
+  EsNodB = 10*log10(sig_var/noise_var);  
+endfunction
+
+
+#{
+  Plot curves of Es/No estimator in action.
+  
+  Plots indicate it works OK down to Es/No=3dB,
+  where it is 1dB high.  That's acceptable as Es/No=3dB is the lower limit of 
+  our operation (ie Eb/No=0dB, 10% raw BER).
+#}
+
+function [EsNo_est rx_symbols] = esno_est_curves(EsNodB=0:20, channel="awgn", plot_en=1)
+    Nsym=1000; rand('seed',1); randn('seed',1);
+    tx_symbols = 2*(rand(1,Nsym) > 0.5) -1 + j*(2*(rand(1,Nsym) > 0.5) - 1);
+    tx_symbols *= exp(-j*pi/4)/sqrt(2);
+    
+    if strcmp(channel,"mpp")
+        % for fading we assume perfect phase recovery, so just multiply by mag
+        spread = doppler_spread(2.0, 50, length(tx_symbols));
+        tx_symbols = tx_symbols .* abs(spread);
+        % normalise power over the multipath channel run
+        S = tx_symbols*tx_symbols';
+        tx_symbols *= sqrt(Nsym/S);
+    end
+    
+    for i = 1:length(EsNodB)
+        aEsNodB = EsNodB(i);
+        EsNo = 10 .^ (aEsNodB/10);
+        N = 1/EsNo;
+        noise = sqrt(N/2)*randn(1,Nsym) +  sqrt(N/2)*j*randn(1,Nsym);  
+        S = tx_symbols*tx_symbols';
+        N = noise*noise';
+        EsNo_meas(i) = 10*log10(S/N);
+        rx_symbols = tx_symbols + noise;  
+        EsNo_est(i) = esno_est_calc(rx_symbols);
+        printf("EsNo: %5.2f EsNo_meas: %5.2f EsNo_est: %5.2f\n", aEsNodB, EsNo_meas(i), EsNo_est(i));
+    end
+    if plot_en
+        figure(1);
+        plot(EsNodB, EsNo_meas, '+-;EsNo meas;');  hold on; plot(EsNodB, EsNo_est, 'o-;EsNo est;'); hold off;
+        axis([0 max(EsNodB) 0 max(EsNodB)]); grid;
+        figure(2); plot(tx_symbols,'+');
+    end
+endfunction
+
+function esno_est_test(channel="awgn")
+    test_point_dB = 5;
+    [EsNo_est_awgn rx_syms] = esno_est_curves(test_point_dB, channel, plot_en=0);
+    if abs(EsNo_est_awgn - test_point_dB) < 1.0
+        printf("%s Pass\n",toupper(channel));
+    else
+        printf("%s Fail\n",toupper(channel));
+    end
+endfunction
+
+function esno_est_tests_octave
+    esno_est_test("awgn");
+    esno_est_test("mpp");    
+endfunction
+
+function esno_est_test_c(channel="awgn")
+    test_point_dB = 5;
+    [EsNo_est rx_syms] = esno_est_curves(test_point_dB, channel, plot_en=0);
+    rx_syms_float = zeros(1,2*length(rx_syms));
+    rx_syms_float(1:2:end) = real(rx_syms);
+    rx_syms_float(2:2:end) = imag(rx_syms);
+    f = fopen("esno_est.iqfloat","wb"); fwrite(f, rx_syms_float, "float"); fclose(f);
+
+    printf("\nRunning C version....\n");
+    path_to_unittest = "../build_linux/unittest"
+    if getenv("PATH_TO_UNITEST")
+      path_to_unittest = getenv("PATH_TO_UNITEST")
+      printf("setting path from env var to %s\n", path_to_unittest);
+    end
+    system(sprintf("%s/tesno_est %s %d > tesno_est_out.txt", path_to_unittest, "esno_est.iqfloat", length(rx_syms)));
+    load tesno_est_out.txt;
+    printf("test_point: %5.2f Octave: %5.2f C: %5.2f\n", test_point_dB, EsNo_est, tesno_est_out);
+    if abs(EsNo_est - tesno_est_out) < 0.5
+        printf("%s Pass\n",toupper(channel));
+    else
+        printf("%s Fail\n",toupper(channel));
+    end
+endfunction
+
+function esno_est_tests_c
+    esno_est_test_c("awgn");
+    esno_est_test_c("mpp");    
+endfunction
+
+