From 4d64ef6e48dc69f5f475c1289fac4b26d7e87bb4 Mon Sep 17 00:00:00 2001
From: Christian Zietz <czietz@gmx.net>
Date: Sun, 25 Oct 2020 13:31:08 +0100
Subject: [PATCH] Compensate IFFT window / zero-padding loss in TD
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Depending on the chosen window and mode, the magnitude of the
impulse response in time-domain previously was to low. This can
be explained by looking at the signal processing. For example,
in bandpass mode with normal window, it applies a 101 point Kaiser
window (shape factor 6) and zero-pads to do a 256 point IFFT.
Therefore, the loss is 20*log10(256/sum(kaiser(101,6))) ≈ 14.2 dB.

This change compensates the signal processing losses in bandpass
and lowpass impulse mode depending on the window type, which makes
the time-domain results similar to other VNAs.
---
 main.c | 46 ++++++++++++++++++++++++++++++----------------
 1 file changed, 30 insertions(+), 16 deletions(-)

diff --git a/main.c b/main.c
index 5dea463..d3d6a9d 100644
--- a/main.c
+++ b/main.c
@@ -198,40 +198,54 @@ transform_domain(void)
   // and calculate ifft for time domain
   float* tmp = (float*)spi_buffer;
 
+  // correct IFFT window and zero-padding loss
+  // assuming a 2*POINT_COUNT window
+  float wincorr = 1.0f;
+  float beta = 0.0;
+  switch (domain_mode & TD_WINDOW) {
+    case TD_WINDOW_MINIMUM:
+      beta = 0.0; // this is rectangular
+      // loss by zero-padding 202 to 256 points
+      wincorr = (float)FFT_SIZE / (float)(2*POINTS_COUNT);
+      break;
+    case TD_WINDOW_NORMAL:
+      beta = 6.0;
+      // additional window loss: 1/mean(kaiser(202,6)) = 2.01
+      wincorr = (float)FFT_SIZE / (float)(2*POINTS_COUNT) * 2.01f;
+      break;
+    case TD_WINDOW_MAXIMUM:
+      beta = 13;
+      // additional window loss: 1/mean(kaiser(202,13)) = 2.92
+      wincorr = (float)FFT_SIZE / (float)(2*POINTS_COUNT) * 2.92f;
+      break;
+  }
+
   uint8_t window_size = POINTS_COUNT, offset = 0;
   uint8_t is_lowpass = FALSE;
   switch (domain_mode & TD_FUNC) {
     case TD_FUNC_BANDPASS:
       offset = 0;
       window_size = POINTS_COUNT;
+      // window size is half the size as assumed above => twice the IFFT loss
+      wincorr *= 2.0f;
       break;
-    case TD_FUNC_LOWPASS_IMPULSE:
     case TD_FUNC_LOWPASS_STEP:
+      // no IFFT losses need to be considered to calculate the step response
+      wincorr = 1.0f;
+      // fall-through
+    case TD_FUNC_LOWPASS_IMPULSE:
       is_lowpass = TRUE;
       offset = POINTS_COUNT;
       window_size = POINTS_COUNT * 2;
       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]));
     for (int i = 0; i < POINTS_COUNT; i++) {
       float w = kaiser_window(i + offset, window_size, beta);
-      tmp[i * 2 + 0] *= w;
-      tmp[i * 2 + 1] *= w;
+      tmp[i * 2 + 0] *= w * wincorr;
+      tmp[i * 2 + 1] *= w * wincorr;
     }
     for (int i = POINTS_COUNT; i < FFT_SIZE; i++) {
       tmp[i * 2 + 0] = 0.0;