Merge pull request #32 from cho45/timedomain

[proposal] Time Domain Operation

fft.h

@@ -0,0 +1,107 @@
+/* 
+ * fft.h is Based on
+ * Free FFT and convolution (C)
+ * 
+ * Copyright (c) 2019 Project Nayuki. (MIT License)
+ * https://www.nayuki.io/page/free-small-fft-in-multiple-languages
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of
+ * this software and associated documentation files (the "Software"), to deal in
+ * the Software without restriction, including without limitation the rights to
+ * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
+ * the Software, and to permit persons to whom the Software is furnished to do so,
+ * subject to the following conditions:
+ * - The above copyright notice and this permission notice shall be included in
+ *   all copies or substantial portions of the Software.
+ * - The Software is provided "as is", without warranty of any kind, express or
+ *   implied, including but not limited to the warranties of merchantability,
+ *   fitness for a particular purpose and noninfringement. In no event shall the
+ *   authors or copyright holders be liable for any claim, damages or other
+ *   liability, whether in an action of contract, tort or otherwise, arising from,
+ *   out of or in connection with the Software or the use or other dealings in the
+ *   Software.
+ */
+
+
+#include <math.h>
+#include <stdint.h>
+
+static uint8_t reverse_bits(uint8_t x, int n) {
+	uint8_t result = 0;
+	for (int i = 0; i < n; i++, x >>= 1)
+		result = (result << 1) | (x & 1U);
+	return result;
+}
+
+static const float sin_table[] = {
+	/*
+	 * float has about 7.2 digits of precision
+		for (uint8_t i = 0; i < 96; i++) {
+			printf("% .8f,%c", sin(2 * M_PI * i / n), i % 8 == 7 ? '\n' : ' ');
+		}
+	*/
+	 0.00000000,  0.04906767,  0.09801714,  0.14673047,  0.19509032,  0.24298018,  0.29028468,  0.33688985,
+	 0.38268343,  0.42755509,  0.47139674,  0.51410274,  0.55557023,  0.59569930,  0.63439328,  0.67155895,
+	 0.70710678,  0.74095113,  0.77301045,  0.80320753,  0.83146961,  0.85772861,  0.88192126,  0.90398929,
+	 0.92387953,  0.94154407,  0.95694034,  0.97003125,  0.98078528,  0.98917651,  0.99518473,  0.99879546,
+	 1.00000000,  0.99879546,  0.99518473,  0.98917651,  0.98078528,  0.97003125,  0.95694034,  0.94154407,
+	 0.92387953,  0.90398929,  0.88192126,  0.85772861,  0.83146961,  0.80320753,  0.77301045,  0.74095113,
+	 0.70710678,  0.67155895,  0.63439328,  0.59569930,  0.55557023,  0.51410274,  0.47139674,  0.42755509,
+	 0.38268343,  0.33688985,  0.29028468,  0.24298018,  0.19509032,  0.14673047,  0.09801714,  0.04906767,
+	 0.00000000, -0.04906767, -0.09801714, -0.14673047, -0.19509032, -0.24298018, -0.29028468, -0.33688985,
+	-0.38268343, -0.42755509, -0.47139674, -0.51410274, -0.55557023, -0.59569930, -0.63439328, -0.67155895,
+	-0.70710678, -0.74095113, -0.77301045, -0.80320753, -0.83146961, -0.85772861, -0.88192126, -0.90398929,
+	-0.92387953, -0.94154407, -0.95694034, -0.97003125, -0.98078528, -0.98917651, -0.99518473, -0.99879546,
+};
+
+/***
+ * dir = forward: 0, inverse: 1
+ * https://www.nayuki.io/res/free-small-fft-in-multiple-languages/fft.c
+ */
+static void fft128(float array[][2], const uint8_t dir) {
+	const uint8_t n = 128;
+	const uint8_t levels = 7; // log2(n)
+	const float* const cos_table = &sin_table[32];
+
+	const uint8_t real =   dir & 1;
+	const uint8_t imag = ~real & 1;
+
+	for (uint8_t i = 0; i < n; i++) {
+		uint8_t j = reverse_bits(i, levels);
+		if (j > i) {
+			float temp = array[i][real];
+			array[i][real] = array[j][real];
+			array[j][real] = temp;
+			temp = array[i][imag];
+			array[i][imag] = array[j][imag];
+			array[j][imag] = temp;
+		}
+	}
+
+	// Cooley-Tukey decimation-in-time radix-2 FFT
+	for (uint8_t size = 2; size <= n; size *= 2) {
+		uint8_t halfsize = size / 2;
+		uint8_t tablestep = n / size;
+		for (uint8_t i = 0; i < n; i += size) {
+			for (uint8_t j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
+				uint8_t l = j + halfsize;
+				float tpre =  array[l][real] * cos_table[k] + array[l][imag] * sin_table[k];
+				float tpim = -array[l][real] * sin_table[k] + array[l][imag] * cos_table[k] ;
+				array[l][real] = array[j][real] - tpre;
+				array[l][imag] = array[j][imag] - tpim;
+				array[j][real] += tpre;
+				array[j][imag] += tpim;
+			}
+		}
+		if (size == n)  // Prevent overflow in 'size *= 2'
+			break;
+	}
+}
+
+static inline void fft128_forward(float array[][2]) {
+	fft128(array, 0);
+}
+
+static inline void fft128_inverse(float array[][2]) {
+	fft128(array, 1);
+}

main.c

@@ -23,6 +23,7 @@
 #include "usbcfg.h"
 #include "si5351.h"
 #include "nanovna.h"
+#include "fft.h"
 
 #include <chprintf.h>
 #include <shell.h>
@@ -107,6 +108,91 @@ toggle_sweep(void)
   sweep_enabled = !sweep_enabled;
 }
 
+float bessel0(float x) {
+	const float eps = 0.0001;
+
+	float ret = 0;
+	float term = 1;
+	float m = 0;
+
+	while (term  > eps * ret) {
+		ret += term;
+		++m;
+		term *= (x*x) / (4*m*m);
+	}
+
+	return ret;
+}
+
+float kaiser_window(float k, float n, float beta) {
+	if (beta == 0.0) return 1.0;
+	float r = (2 * k) / (n - 1) - 1;
+	return bessel0(beta * sqrt(1 - r * r)) / bessel0(beta);
+}
+
+static
+void
+transform_domain(void)
+{
+  if ((domain_mode & DOMAIN_MODE) != DOMAIN_TIME) return; // nothing to do for freq domain
+  // use spi_buffer as temporary buffer
+  // and calculate ifft for time domain
+  float* tmp = (float*)spi_buffer;
+
+  uint8_t window_size, offset;
+  switch (domain_mode & TD_FUNC) {
+      case TD_FUNC_BANDPASS:
+          offset = 0;
+          window_size = 101;
+          break;
+      case TD_FUNC_LOWPASS_IMPULSE:
+      case TD_FUNC_LOWPASS_STEP:
+          offset = 101;
+          window_size = 202;
+          break;
+  }
+
+  float beta = 0.0;
+  switch (domain_mode & TD_WINDOW) {
+      case TD_WINDOW_MINIMUM:
+          beta = 0.0; // this is rectangular
+          break;
+      case TD_WINDOW_NORMAL:
+          beta = 6.0;
+          break;
+      case TD_WINDOW_MAXIMUM:
+          beta = 13;
+          break;
+  }
+
+  for (int ch = 0; ch < 2; ch++) {
+      memcpy(tmp, measured[ch], sizeof(measured[0]));
+      if (beta != 0.0) {
+          for (int i = 0; i < 101; i++) {
+              float w = kaiser_window(i+offset, window_size, beta);
+              tmp[i*2+0] *= w;
+              tmp[i*2+1] *= w;
+          }
+      }
+      for (int i = 101; i < 128; i++) {
+          tmp[i*2+0] = 0.0;
+          tmp[i*2+1] = 0.0;
+      }
+      fft128_inverse((float(*)[2])tmp);
+      memcpy(measured[ch], tmp, sizeof(measured[0]));
+      for (int i = 0; i < 101; i++) {
+          measured[ch][i][0] /= 128.0;
+          measured[ch][i][1] /= 128.0;
+      }
+      if ( (domain_mode & TD_FUNC) == TD_FUNC_LOWPASS_STEP ) {
+          for (int i = 1; i < 101; i++) {
+              measured[ch][i][0] += measured[ch][i-1][0];
+              measured[ch][i][1] += measured[ch][i-1][1];
+          }
+      }
+  }
+}
+
 static void cmd_pause(BaseSequentialStream *chp, int argc, char *argv[])
 {
     (void)chp;
@@ -515,6 +601,8 @@ properties_t current_props = {
     { 1, 30, 0 }, { 0, 40, 0 }, { 0, 60, 0 }, { 0, 80, 0 }
   },
   /* active_marker */      0,
+  /* domain_mode */        0,
+  /* velocity_factor */   70,
   /* checksum */           0
 };
 properties_t *active_props = &current_props;
@@ -599,11 +687,13 @@ void sweep(void)
     redraw_requested = FALSE;
     ui_process();
     if (redraw_requested)
-      return; // return to redraw screen asap.
+      break; // return to redraw screen asap.
       
     if (frequency_updated)
       goto rewind;
   }
+
+  transform_domain();
 }
 
 static void

nanovna.h

@@ -22,6 +22,7 @@
 /*
  * main.c
  */
+
 extern float measured[2][101][2];
 
 #define CAL_LOAD 0
@@ -49,6 +50,18 @@ extern float measured[2][101][2];
 #define ETERM_ET 3 /* error term transmission tracking */
 #define ETERM_EX 4 /* error term isolation */
 
+#define DOMAIN_MODE (1<<0)
+#define DOMAIN_FREQ (0<<0)
+#define DOMAIN_TIME (1<<0)
+#define TD_FUNC (0b11<<1)
+#define TD_FUNC_BANDPASS (0b00<<1)
+#define TD_FUNC_LOWPASS_IMPULSE (0b01<<1)
+#define TD_FUNC_LOWPASS_STEP    (0b10<<1)
+#define TD_WINDOW (0b11<<3)
+#define TD_WINDOW_NORMAL (0b00<<3)
+#define TD_WINDOW_MINIMUM (0b01<<3)
+#define TD_WINDOW_MAXIMUM (0b10<<3)
+
 void cal_collect(int type);
 void cal_done(void);
 
@@ -282,6 +295,8 @@ typedef struct {
   trace_t _trace[TRACES_MAX];
   marker_t _markers[4];
   int _active_marker;
+  uint8_t _domain_mode; /* 0bxxxxxffm : where ff: TD_FUNC m: DOMAIN_MODE */
+  uint8_t _velocity_factor; // %
 
   int32_t checksum;
 } properties_t;
@@ -305,6 +320,8 @@ extern int8_t previous_marker;
 #define trace current_props._trace
 #define markers current_props._markers
 #define active_marker current_props._active_marker
+#define domain_mode current_props._domain_mode
+#define velocity_factor current_props._velocity_factor
 
 int caldata_save(int id);
 int caldata_recall(int id);

plot.c

@@ -701,6 +701,17 @@ trace_get_info(int t, char *buf, int len)
   }
 }
 
+static float time_of_index(int idx) {
+   return 1.0 / (float)(frequencies[1] - frequencies[0]) / 128.0 * idx;
+}
+
+static float distance_of_index(int idx) {
+#define SPEED_OF_LIGHT 299792458
+   float distance = ((float)idx * (float)SPEED_OF_LIGHT) / ( (float)(frequencies[1] - frequencies[0]) * 128.0 * 2.0);
+   return distance * (velocity_factor / 100.0);
+}
+
+
 static inline void
 mark_map(int x, int y)
 {
@@ -1367,8 +1378,13 @@ cell_draw_marker_info(int m, int n, int w, int h)
   chsnprintf(buf, sizeof buf, "%d:", active_marker + 1);
   cell_drawstring_5x7(w, h, buf, xpos, ypos, 0xffff);
   xpos += 16;
-  frequency_string(buf, sizeof buf, frequencies[idx]);
+  if ((domain_mode & DOMAIN_MODE) == DOMAIN_FREQ) {
-  cell_drawstring_5x7(w, h, buf, xpos, ypos, 0xffff);
+      frequency_string(buf, sizeof buf, frequencies[idx]);
+      cell_drawstring_5x7(w, h, buf, xpos, ypos, 0xffff);
+  } else {
+      chsnprintf(buf, sizeof buf, "%d ns %.1f m", (uint16_t)(time_of_index(idx) * 1e9), distance_of_index(idx));
+      cell_drawstring_5x7(w, h, buf, xpos, ypos, 0xffff);
+  }
 
   // draw marker delta
   if (previous_marker >= 0 && active_marker != previous_marker && markers[previous_marker].enabled) {
@@ -1410,37 +1426,47 @@ void
 draw_frequencies(void)
 {
   char buf[24];
-  if (frequency1 > 0) {
+  if ((domain_mode & DOMAIN_MODE) == DOMAIN_FREQ) {
-    int start = frequency0;
+      if (frequency1 > 0) {
-    int stop = frequency1;
+        int start = frequency0;
-    strcpy(buf, "START ");
+        int stop = frequency1;
-    frequency_string(buf+6, 24-6, start);
+        strcpy(buf, "START ");
-    strcat(buf, "    ");
+        frequency_string(buf+6, 24-6, start);
-    ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000);
+        strcat(buf, "    ");
-    strcpy(buf, "STOP ");
+        ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000);
-    frequency_string(buf+5, 24-5, stop);
+        strcpy(buf, "STOP ");
-    strcat(buf, "    ");
+        frequency_string(buf+5, 24-5, stop);
-    ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000);
+        strcat(buf, "    ");
-  } else if (frequency1 < 0) {
+        ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000);
-    int fcenter = frequency0;
+      } else if (frequency1 < 0) {
-    int fspan = -frequency1;
+        int fcenter = frequency0;
-    strcpy(buf, "CENTER ");
+        int fspan = -frequency1;
-    frequency_string(buf+7, 24-7, fcenter);
+        strcpy(buf, "CENTER ");
-    strcat(buf, "    ");
+        frequency_string(buf+7, 24-7, fcenter);
-    ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000);
+        strcat(buf, "    ");
-    strcpy(buf, "SPAN ");
+        ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000);
-    frequency_string(buf+5, 24-5, fspan);
+        strcpy(buf, "SPAN ");
-    strcat(buf, "    ");
+        frequency_string(buf+5, 24-5, fspan);
-    ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000);
+        strcat(buf, "    ");
+        ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000);
+      } else {
+        int fcenter = frequency0;
+        chsnprintf(buf, 24, "CW %d.%03d %03d MHz    ",
+                   (int)(fcenter / 1000000),
+                   (int)((fcenter / 1000) % 1000),
+                   (int)(fcenter % 1000));
+        ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000);
+        chsnprintf(buf, 24, "                             ");
+        ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000);
+      }
   } else {
-    int fcenter = frequency0;
+      strcpy(buf, "START 0s        ");
-    chsnprintf(buf, 24, "CW %d.%03d %03d MHz    ",
+      ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000);
-               (int)(fcenter / 1000000),
+
-               (int)((fcenter / 1000) % 1000),
+      strcpy(buf, "STOP ");
-               (int)(fcenter % 1000));
+      chsnprintf(buf+5, 24-5, "%d ns", (uint16_t)(time_of_index(101) * 1e9));
-    ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000);
+      strcat(buf, "          ");
-    chsnprintf(buf, 24, "                             ");
+      ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000);
-    ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000);
   }
 }
 

ui.c

@@ -68,7 +68,7 @@ enum {
 };
 
 enum {
-  KM_START, KM_STOP, KM_CENTER, KM_SPAN, KM_CW, KM_SCALE, KM_REFPOS, KM_EDELAY
+  KM_START, KM_STOP, KM_CENTER, KM_SPAN, KM_CW, KM_SCALE, KM_REFPOS, KM_EDELAY, KM_VELOCITY_FACTOR
 };
 
 uint8_t ui_mode = UI_NORMAL;
@@ -665,6 +665,65 @@ menu_channel_cb(int item)
   ui_mode_normal();
 }
 
+static void
+menu_transform_window_cb(int item)
+{
+  // TODO
+  switch (item) {
+    case 0:
+      domain_mode = (domain_mode & ~TD_WINDOW) | TD_WINDOW_MINIMUM;
+      ui_mode_normal();
+      break;
+    case 1:
+      domain_mode = (domain_mode & ~TD_WINDOW) | TD_WINDOW_NORMAL;
+      ui_mode_normal();
+      break;
+    case 2:
+      domain_mode = (domain_mode & ~TD_WINDOW) | TD_WINDOW_MAXIMUM;
+      ui_mode_normal();
+      break;
+  }
+}
+
+static void
+menu_transform_cb(int item)
+{
+  int status;
+  switch (item) {
+    case 0:
+      if ((domain_mode & DOMAIN_MODE) == DOMAIN_TIME) {
+          domain_mode = (domain_mode & ~DOMAIN_MODE) | DOMAIN_FREQ;
+      } else {
+          domain_mode = (domain_mode & ~DOMAIN_MODE) | DOMAIN_TIME;
+      }
+      draw_frequencies();
+      ui_mode_normal();
+      break;
+    case 1:
+      domain_mode = (domain_mode & ~TD_FUNC) | TD_FUNC_LOWPASS_IMPULSE;
+      ui_mode_normal();
+      break;
+    case 2:
+      domain_mode = (domain_mode & ~TD_FUNC) | TD_FUNC_LOWPASS_STEP;
+      ui_mode_normal();
+      break;
+    case 3:
+      domain_mode = (domain_mode & ~TD_FUNC) | TD_FUNC_BANDPASS;
+      ui_mode_normal();
+      break;
+    case 5:
+      status = btn_wait_release();
+      if (status & EVT_BUTTON_DOWN_LONG) {
+        ui_mode_numeric(KM_VELOCITY_FACTOR);
+        ui_process_numeric();
+      } else {
+        ui_mode_keypad(KM_VELOCITY_FACTOR);
+        ui_process_keypad();
+      }
+      break;
+  }
+}
+
 static void 
 choose_active_marker(void)
 {
@@ -875,11 +934,31 @@ const menuitem_t menu_channel[] = {
   { MT_NONE, NULL, NULL } // sentinel
 };
 
+const menuitem_t menu_transform_window[] = {
+  { MT_CALLBACK, "MINIMUM", menu_transform_window_cb },
+  { MT_CALLBACK, "NORMAL", menu_transform_window_cb },
+  { MT_CALLBACK, "MAXIMUM", menu_transform_window_cb },
+  { MT_CANCEL, S_LARROW" BACK", NULL },
+  { MT_NONE, NULL, NULL } // sentinel
+};
+
+const menuitem_t menu_transform[] = {
+  { MT_CALLBACK, "\2TRANSFORM\0ON", menu_transform_cb },
+  { MT_CALLBACK, "\2LOW PASS\0IMPULSE", menu_transform_cb },
+  { MT_CALLBACK, "\2LOW PASS\0STEP", menu_transform_cb },
+  { MT_CALLBACK, "BANDPASS", menu_transform_cb },
+  { MT_SUBMENU, "WINDOW", menu_transform_window },
+  { MT_CALLBACK, "\2VELOCITY\0FACTOR", menu_transform_cb },
+  { MT_CANCEL, S_LARROW" BACK", NULL },
+  { MT_NONE, NULL, NULL } // sentinel
+};
+
 const menuitem_t menu_display[] = {
   { MT_SUBMENU, "TRACE", menu_trace },
   { MT_SUBMENU, "FORMAT", menu_format },
   { MT_SUBMENU, "SCALE", menu_scale },
   { MT_SUBMENU, "CHANNEL", menu_channel },
+  { MT_SUBMENU, "TRANSFORM", menu_transform },
   { MT_CANCEL, S_LARROW" BACK", NULL },
   { MT_NONE, NULL, NULL } // sentinel
 };
@@ -1120,11 +1199,12 @@ const keypads_t * const keypads_mode_tbl[] = {
   keypads_freq, // cw freq
   keypads_scale, // scale
   keypads_scale, // respos
-  keypads_time // electrical delay
+  keypads_time, // electrical delay
+  keypads_scale // velocity factor
 };
 
 const char * const keypad_mode_label[] = {
-  "START", "STOP", "CENTER", "SPAN", "CW FREQ", "SCALE", "REFPOS", "EDELAY"
+  "START", "STOP", "CENTER", "SPAN", "CW FREQ", "SCALE", "REFPOS", "EDELAY", "VELOCITY%"
 };
 
 void
@@ -1223,6 +1303,7 @@ menu_item_modify_attribute(const menuitem_t *menu, int item,
         || (item == 2 && (cal_status & CALSTAT_LOAD))
         || (item == 3 && (cal_status & CALSTAT_ISOLN))
         || (item == 4 && (cal_status & CALSTAT_THRU))) {
+      domain_mode = (domain_mode & ~DOMAIN_MODE) | DOMAIN_FREQ;
       *bg = 0x0000;
       *fg = 0xffff;
     }
@@ -1236,6 +1317,23 @@ menu_item_modify_attribute(const menuitem_t *menu, int item,
       *bg = 0x0000;
       *fg = 0xffff;
     }
+  } else if (menu == menu_transform) {
+      if ((item == 0 && (domain_mode & DOMAIN_MODE) == DOMAIN_TIME)
+       || (item == 1 && (domain_mode & TD_FUNC) == TD_FUNC_LOWPASS_IMPULSE)
+       || (item == 2 && (domain_mode & TD_FUNC) == TD_FUNC_LOWPASS_STEP)
+       || (item == 3 && (domain_mode & TD_FUNC) == TD_FUNC_BANDPASS)
+       ) {
+        *bg = 0x0000;
+        *fg = 0xffff;
+      }
+  } else if (menu == menu_transform_window) {
+      if ((item == 0 && (domain_mode & TD_WINDOW) == TD_WINDOW_MINIMUM)
+       || (item == 1 && (domain_mode & TD_WINDOW) == TD_WINDOW_NORMAL)
+       || (item == 2 && (domain_mode & TD_WINDOW) == TD_WINDOW_MAXIMUM)
+       ) {
+        *bg = 0x0000;
+        *fg = 0xffff;
+      }
   }
 }
 
@@ -1363,6 +1461,9 @@ fetch_numeric_target(void)
   case KM_EDELAY:
     uistat.value = get_electrical_delay();
     break;
+  case KM_VELOCITY_FACTOR:
+    uistat.value = velocity_factor;
+    break;
   }
   
   {
@@ -1402,6 +1503,9 @@ void set_numeric_value(void)
   case KM_EDELAY:
     set_electrical_delay(uistat.value);
     break;
+  case KM_VELOCITY_FACTOR:
+    velocity_factor = uistat.value;
+    break;
   }
 }
 
@@ -1576,6 +1680,9 @@ keypad_click(int key)
     case KM_EDELAY:
       set_electrical_delay(value); // pico seconds
       break;
+    case KM_VELOCITY_FACTOR:
+      velocity_factor = value;
+      break;
     }
 
     return KP_DONE;