rtl_fm: esbensen's discriminant

src/rtl_fm.c

@@ -215,7 +215,7 @@ void usage(void)
 		"\t[-F fir_size (default: off)]\n"
 		"\t    enables low-leakage downsample filter\n"
 		"\t    size can be 0 or 9.  0 has bad roll off\n"
-		"\t[-A std/fast/lut choose atan math (default: std)]\n"
+		"\t[-A std/fast/lut/ale choose atan math (default: std)]\n"
 		//"\t[-C clip_path (default: off)\n"
 		//"\t (create time stamped raw clips, requires squelch)\n"
 		//"\t (path must have '\%s' and will expand to date_time_freq)\n"
@@ -515,6 +515,24 @@ int polar_disc_lut(int ar, int aj, int br, int bj)
 	return 0;
 }
 
+int esbensen(int ar, int aj, int br, int bj)
+/*
+  input signal: s(t) = a*exp(-i*w*t+p)
+  a = amplitude, w = angular freq, p = phase difference
+  solve w
+  s' = -i(w)*a*exp(-i*w*t+p)
+  s'*conj(s) = -i*w*a*a
+  s'*conj(s) / |s|^2 = -i*w
+*/
+{
+	int cj, dr, dj;
+	int scaled_pi = 2608; /* 1<<14 / (2*pi) */
+	dr = (br - ar) * 2;
+	dj = (bj - aj) * 2;
+	cj = bj*dr - br*dj; /* imag(ds*conj(s)) */
+	return (scaled_pi * cj / (ar*ar+aj*aj+1));
+}
+
 void fm_demod(struct demod_state *fm)
 {
 	int i, pcm;
@@ -536,6 +554,10 @@ void fm_demod(struct demod_state *fm)
 			pcm = polar_disc_lut(lp[i], lp[i+1],
 				lp[i-2], lp[i-1]);
 			break;
+		case 3:
+			pcm = esbensen(lp[i], lp[i+1],
+				lp[i-2], lp[i-1]);
+			break;
 		}
 		fm->result[i/2] = (int16_t)pcm;
 	}
@@ -1122,6 +1144,8 @@ int main(int argc, char **argv)
 			if (strcmp("lut",  optarg) == 0) {
 				atan_lut_init();
 				demod.custom_atan = 2;}
+			if (strcmp("ale", optarg) == 0) {
+				demod.custom_atan = 3;}
 			break;
 		case 'M':
 			if (strcmp("fm",  optarg) == 0) {