diff --git a/fft.h b/fft.h new file mode 100644 index 0000000..1a46515 --- /dev/null +++ b/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 +#include + +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); +} diff --git a/main.c b/main.c index 4cff0df..3653ca0 100644 --- a/main.c +++ b/main.c @@ -23,6 +23,7 @@ #include "usbcfg.h" #include "si5351.h" #include "nanovna.h" +#include "fft.h" #include #include @@ -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 = ¤t_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 diff --git a/nanovna.h b/nanovna.h index 655e78e..e60924d 100644 --- a/nanovna.h +++ b/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); diff --git a/plot.c b/plot.c index 2367372..83ed0fb 100644 --- a/plot.c +++ b/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]); - cell_drawstring_5x7(w, h, buf, xpos, ypos, 0xffff); + if ((domain_mode & DOMAIN_MODE) == DOMAIN_FREQ) { + 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) { - int start = frequency0; - int stop = frequency1; - strcpy(buf, "START "); - frequency_string(buf+6, 24-6, start); - strcat(buf, " "); - ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000); - strcpy(buf, "STOP "); - frequency_string(buf+5, 24-5, stop); - strcat(buf, " "); - ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000); - } else if (frequency1 < 0) { - int fcenter = frequency0; - int fspan = -frequency1; - strcpy(buf, "CENTER "); - frequency_string(buf+7, 24-7, fcenter); - strcat(buf, " "); - ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000); - strcpy(buf, "SPAN "); - frequency_string(buf+5, 24-5, fspan); - strcat(buf, " "); - ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000); + if ((domain_mode & DOMAIN_MODE) == DOMAIN_FREQ) { + if (frequency1 > 0) { + int start = frequency0; + int stop = frequency1; + strcpy(buf, "START "); + frequency_string(buf+6, 24-6, start); + strcat(buf, " "); + ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000); + strcpy(buf, "STOP "); + frequency_string(buf+5, 24-5, stop); + strcat(buf, " "); + ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000); + } else if (frequency1 < 0) { + int fcenter = frequency0; + int fspan = -frequency1; + strcpy(buf, "CENTER "); + frequency_string(buf+7, 24-7, fcenter); + strcat(buf, " "); + ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000); + strcpy(buf, "SPAN "); + frequency_string(buf+5, 24-5, fspan); + 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; - 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); + strcpy(buf, "START 0s "); + ili9341_drawstring_5x7(buf, OFFSETX, 233, 0xffff, 0x0000); + + strcpy(buf, "STOP "); + chsnprintf(buf+5, 24-5, "%d ns", (uint16_t)(time_of_index(101) * 1e9)); + strcat(buf, " "); + ili9341_drawstring_5x7(buf, 205, 233, 0xffff, 0x0000); } } diff --git a/ui.c b/ui.c index 7eb8576..bc4acd0 100644 --- a/ui.c +++ b/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;