1% estsnr.m
2% David Rowe May 2017
3%
4% estimate SNR of a sinewave in noise
5
6function snr_dB = estsnr(x, Fs=8000, Nbw = 3000)
7
8  [nr nc] = size(x);
9  if nr == 1
10    x = x';
11  end
12
13  % find peak in +ve side of spectrum, ignoring DC
14
15  L = length(x);
16  X = abs(fft(x));
17  st = floor(0.05*L);  en = floor(0.45*L);
18  [A mx_ind]= max(X(st:en));
19  mx_ind += st;
20
21  % signal energy might be spread by doppler, so sum energy
22  % in frequencies +/- 1%
23
24  s_st = floor(mx_ind*0.99); s_en = floor(mx_ind*1.01);
25  S = sum(X(s_st:s_en).^2);
26
27  % real signal, so -ve power is the same
28
29  S = 2*S;
30  SdB = 10*log10(S);
31
32  printf("Signal Power S: %3.2f dB\n", SdB);
33
34  % locate a band of noise next to it and find power in band
35
36  st = floor(mx_ind+0.05*(L/2));
37  en = st + floor(0.1*(L/2));
38
39  N = sum(X(st:en).^2);
40
41  % scale this to obtain total noise power across total bandwidth
42
43  N *= L/(en-st);
44  NdB = 10*log10(N);
45  printf("Noise Power N: %3.2f dB\n", NdB);
46
47  % scale noise to designed noise bandwidth /2 fudge factor as its a
48  % real signal, wish I had a better way to explain that!
49
50  NodB = NdB - 10*log10(Fs/2);
51  NscaleddB = NodB + 10*log10(Nbw);
52  snr_dB = SdB - NscaleddB;
53
54  figure(1); clf;
55  plot(20*log10(X(1:L/2)),'b');
56  hold on;
57  plot([s_st s_en], [NdB NdB]- 10*log10(L), 'r');
58  plot([st en], [NdB NdB]- 10*log10(L), 'r');
59  hold off;
60  top = 10*ceil(SdB/10);
61  bot = NodB - 20;
62  axis([1 L/2 bot top]);
63  grid
64  grid("minor")
65endfunction
66