From 526f6dd41518fe86b6972d9ad88114ca7f447036 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Sat, 27 May 2017 21:54:31 +0000
Subject: [PATCH] SNR estimation for EME

git-svn-id: https://svn.code.sf.net/p/freetel/code@3143 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/octave/estsnr.m        | 65 +++++++++++++++++++++++++++++++
 codec2-dev/octave/plot_specgram.m | 30 ++++++++++++++
 2 files changed, 95 insertions(+)
 create mode 100644 codec2-dev/octave/estsnr.m
 create mode 100644 codec2-dev/octave/plot_specgram.m

diff --git a/codec2-dev/octave/estsnr.m b/codec2-dev/octave/estsnr.m
new file mode 100644
index 00000000..5a00bd80
--- /dev/null
+++ b/codec2-dev/octave/estsnr.m
@@ -0,0 +1,65 @@
+% estsnr.m
+% David Rowe May 2017
+%
+% estimate SNR of a sinewave in noise
+
+function snr_dB = estsnr(x, Fs=8000, Nbw = 3000)
+
+  [nr nc] = size(x);
+  if nr == 1
+    x = x';
+  end
+
+  % find peak in +ve side of spectrum, ignoring DC
+
+  L = length(x);
+  X = abs(fft(x));
+  st = floor(0.05*L);  en = floor(0.45*L);
+  [A mx_ind]= max(X(st:en));
+  mx_ind += st;
+
+  % signal energy might be spread by doppler, so sum energy
+  % in frequencies +/- 1%
+
+  s_st = floor(mx_ind*0.99); s_en = floor(mx_ind*1.01); 
+  S = sum(X(s_st:s_en).^2);
+
+  % real signal, so -ve power is the same
+
+  S = 2*S;
+  SdB = 10*log10(S);
+
+  printf("Signal Power S: %3.2f dB\n", SdB);
+
+  % locate a band of noise next to it and find power in band
+
+  st = floor(mx_ind+0.05*(L/2));
+  en = st + floor(0.1*(L/2));
+  
+  N = sum(X(st:en).^2);
+
+  % scale this to obtain total noise power across total bandwidth
+
+  N *= L/(en-st);
+  NdB = 10*log10(N);
+  printf("Noise Power N: %3.2f dB\n", NdB);
+
+  % scale noise to designed noise bandwidth /2 fudge factor as its a
+  % real signal, wish I had a better way to explain that!
+
+  NodB = NdB - 10*log10(Fs/2);
+  NscaleddB = NodB + 10*log10(Nbw);
+  snr_dB = SdB - NscaleddB;
+
+  figure(1); clf;
+  plot(20*log10(X(1:L/2)),'b');
+  hold on;
+  plot([s_st s_en], [NdB NdB]- 10*log10(L), 'r');
+  plot([st en], [NdB NdB]- 10*log10(L), 'r');
+  hold off;
+  top = 10*ceil(SdB/10);
+  bot = NodB - 20;
+  axis([1 L/2 bot top]);
+  grid
+  grid("minor")
+endfunction
diff --git a/codec2-dev/octave/plot_specgram.m b/codec2-dev/octave/plot_specgram.m
new file mode 100644
index 00000000..480253eb
--- /dev/null
+++ b/codec2-dev/octave/plot_specgram.m
@@ -0,0 +1,30 @@
+% plot_specgram.m
+% David Rowe May 2017
+%
+% As the name suggests.....
+
+
+function S = plot_specgram(x)
+
+  % set step size so we end up with an aray of pixels that is square.  Otherwise
+  % imagesc will only use a ribbon on the figure.
+
+  Fs = 8000;
+  l = length(x);
+
+  Nfft = 512; Nfft2 = Nfft/2; window = 512; nstep = window/2;
+  nsamples = floor(l/nstep) - 1;
+  S = zeros(nsamples,Nfft);
+  h = hanning(Nfft);
+  for i=1:nsamples
+    st = (i-1)*nstep+1; en = st+window-1;
+    %printf("i: %d st: %d en: %d l: %d\n", i, st, en, l);
+    S(i,:) = 20*log10(abs(fft(x(st:en).*h')));
+  end
+  mx = ceil(max(max(S))/10)*10;
+  S = max(S,mx-30);
+  mesh((1:Nfft2)*4000/Nfft2, (1:nsamples)*nstep, S(:,1:Nfft2));
+  view(90,-90);
+  %xlabel('Frequency (Hz)');
+  %ylabel('Sample');
+en
-- 
2.25.1