Rigol DHO914S Bode plot

[iMo]

There is a "counter" in the scope (when enabled I get it on the screen) but I do not know how to read it.. 6 digits and shows 5 stable digits with my generator.
There are many open issues with my vintage DS1062CA 2GSa/s scope - how to see the ETS actually works (I can set it in menu, up to 50GSa/sec, 10GSa/sec in my model), how to read the 5120 long buffer, how to read the counter, how to hack it to 300MHz 1302, etc..

[iMo]

A divider 1.5Meg/33ohm in a shielded box, scope averaging 8, REF amplitude aprox 11Vpp, DUT ampl aprox 200uVpp on the screen (1mV/div).
xFreq measured with scope's "counter" and used in the math, stability during the measurement within 0.2Hz.
9.5 periods on the screen, AC coupling, 1Meg input, 1x probe, BW limit ON, sinx/x ON, DS1062CA and HP3310B.
For betas 1..30.

[iMo]

A divider 1.5Meg/33ohm in a shielded box, scope averaging NONE, REF amplitude aprox 11Vpp, DUT ampl aprox 200uVpp on the screen (1mV/div), noise aprox 200uVpp.
xFreq measured with scope's "counter" and used in the math, stability during the measurement within 0.1Hz.
9.5 periods on the screen, AC coupling, 1Meg input, 1x probe, BW limit OFF, sinx/x OFF, DS1062CA and HP3310B.
For betas 1..30.
PS: both measurments clicked manually, so the number of measurements for each beta may differ (too lazy with programming python's "for loops" today )

I've tried with 1.5Meg/3.9ohm, the DUT voltage buried in the noise on the screen (my scope uses data off the display, 200points resolution), the results were -111 to -115dB, rather hairy. With 12bit less noisy scope and "ADC raw data" the  -120dB is easily doable, imho..

[gf]

sinx/x OFF

Just wonder, does the scope really honor this parameter for calculating the dispaly data at the timebase you are using here?

For betas 1..30.

I wonder if such a large beta is really useful. Sure, increasing beta increases the stop band attenuation, but it also increases the RBW and equivalent noise bandwidth of the filter, i.e. more noise can pass through the filter. If you are not dealing with strong interfering signals that must be filtered out, a too large beta is rather counterproductive. A beta of 30 also increases the minimum number of periods in the measurement window to 10.

I've tried with 1.5Meg/3.9ohm, the DUT voltage buried in the noise on the screen (my scope uses data off the display, 200points resolution), the results were -111 to -115dB, rather hairy. With 12bit less noisy scope and "ADC raw data" the  -120dB is easily doable, imho..

The aim is of course to use longer records with much more samples, not just 600.

For now, you can try to average the complex gain of (say) 100 measurements at the same frequency, before calculating magnitude and phase.
Or, if the samples you can read from the scope are 8-bit quantized, then maybe even renounce scope averaging and only do compex gain averaging in the app.

[iMo]

I asked my notebook to help me with the python:

How to measure magnitude ratio in decibels and phase in degrees of two sine waves with help of Keiser window and something like DFT at given frequency, using complex exponential function? Could you show me python code?

Code: [Select]

import numpy as np
import matplotlib.pyplot as plt

# Parameters
fs = 1000  # Sampling frequency in Hz
T = 1.0    # Signal duration in seconds
f0 = 50.0  # Frequency of the sine waves in Hz
A1, A2 = 1.0, 0.5  # Amplitudes of the sine waves
phase1, phase2 = 0.0, np.pi/4  # Phases of the sine waves in radians
beta = 8.0  # Kaiser window beta parameter

# Time array
t = np.arange(0, T, 1/fs)

# Generate the sine waves
x1 = A1 * np.sin(2 * np.pi * f0 * t + phase1)
x2 = A2 * np.sin(2 * np.pi * f0 * t + phase2)

# Kaiser window
window = np.kaiser(len(t), beta)

# Apply window to the signals
x1_windowed = x1 * window
x2_windowed = x2 * window

# Compute the DFT at the given frequency using a complex exponential
k = int(f0 * len(t) / fs)  # Index corresponding to frequency f0
exponential = np.exp(-2j * np.pi * k * np.arange(len(t)) / len(t))
X1 = np.dot(x1_windowed, exponential)  # DFT value for x1
X2 = np.dot(x2_windowed, exponential)  # DFT value for x2

# Compute magnitude ratio in decibels
magnitude_ratio_db = 20 * np.log10(np.abs(X2) / np.abs(X1))

# Compute phase difference in degrees
phase_diff_deg = np.angle(X2 / X1, deg=True)

# Print results
print(f"Magnitude Ratio (dB): {magnitude_ratio_db:.2f}")
print(f"Phase Difference (degrees): {phase_diff_deg:.2f}")

# Plot signals and window
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t, x1, label="Sine Wave 1")
plt.plot(t, x2, label="Sine Wave 2")
plt.legend()
plt.title("Original Sine Waves")

plt.subplot(2, 1, 2)
plt.plot(t, x1_windowed, label="Windowed Sine Wave 1")
plt.plot(t, x2_windowed, label="Windowed Sine Wave 2")
plt.legend()
plt.title("Windowed Sine Waves")
plt.tight_layout()
plt.show()

How to measure magnitude ratio in decibels and phase in degrees of two sine waves with help of Keiser window and something like DFT at given frequency, while not using complex exponential function and not using costly sin and cos and arctan2 functions? Could you show me python code?

Code: [Select]

import numpy as np
import matplotlib.pyplot as plt

# Parameters
fs = 1000  # Sampling frequency in Hz
T = 1.0    # Signal duration in seconds
f0 = 50.0  # Frequency of the sine waves in Hz
A1, A2 = 1.0, 0.5  # Amplitudes of the sine waves
phase1, phase2 = 0.0, np.pi/4  # Phases of the sine waves in radians
beta = 8.0  # Kaiser window beta parameter

# Time array
t = np.arange(0, T, 1/fs)

# Generate the sine waves
x1 = A1 * np.sin(2 * np.pi * f0 * t + phase1)
x2 = A2 * np.sin(2 * np.pi * f0 * t + phase2)

# Kaiser window
window = np.kaiser(len(t), beta)

# Apply window to the signals
x1_windowed = x1 * window
x2_windowed = x2 * window

# Initialize recurrence variables for sine and cosine
delta_theta = 2 * np.pi * f0 / fs
cos_curr, sin_curr = 1.0, 0.0  # cos(0) = 1, sin(0) = 0
cos_delta, sin_delta = np.cos(delta_theta), np.sin(delta_theta)

# Initialize accumulators for projections
X1_real = 0.0
X1_imag = 0.0
X2_real = 0.0
X2_imag = 0.0

# Compute projections using recurrence relations
for i in range(len(t)):
    # Accumulate projections
    X1_real += x1_windowed[i] * cos_curr
    X1_imag += x1_windowed[i] * sin_curr
    X2_real += x2_windowed[i] * cos_curr
    X2_imag += x2_windowed[i] * sin_curr
    
    # Update sine and cosine using recurrence
    cos_next = cos_curr * cos_delta - sin_curr * sin_delta
    sin_next = sin_curr * cos_delta + cos_curr * sin_delta
    cos_curr, sin_curr = cos_next, sin_next

# Normalize the projections
X1_real *= 2 / len(t)
X1_imag *= 2 / len(t)
X2_real *= 2 / len(t)
X2_imag *= 2 / len(t)

# Compute magnitudes
X1_mag = np.sqrt(X1_real**2 + X1_imag**2)
X2_mag = np.sqrt(X2_real**2 + X2_imag**2)

# Magnitude ratio in decibels
magnitude_ratio_db = 20 * np.log10(X2_mag / X1_mag)

# Approximate phase angle calculation
def approximate_angle(y, x):
    """Approximate angle using Taylor series expansion."""
    if x == 0:
        return np.pi/2 if y > 0 else -np.pi/2
    elif y == 0:
        return 0 if x > 0 else np.pi
    ratio = y / x
    angle = ratio / (1 + 0.28 * ratio**2)  # Fast approximation
    if x < 0:
        angle += np.pi if y >= 0 else -np.pi
    return angle

# Compute phases using the approximation
X1_phase = approximate_angle(X1_imag, X1_real)
X2_phase = approximate_angle(X2_imag, X2_real)

# Phase difference in degrees
phase_diff_deg = np.degrees(X2_phase - X1_phase)

# Print results
print(f"Magnitude Ratio (dB): {magnitude_ratio_db:.2f}")
print(f"Phase Difference (degrees): {phase_diff_deg:.2f}")

# Plot signals and window
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t, x1, label="Sine Wave 1")
plt.plot(t, x2, label="Sine Wave 2")
plt.legend()
plt.title("Original Sine Waves")

plt.subplot(2, 1, 2)
plt.plot(t, x1_windowed, label="Windowed Sine Wave 1")
plt.plot(t, x2_windowed, label="Windowed Sine Wave 2")
plt.legend()
plt.title("Windowed Sine Waves")
plt.tight_layout()
plt.show()

[gf]

This is something I do not understand:

Code: [Select]

# Compute the DFT at the given frequency using a complex exponential
k = int(f0 * len(t) / fs)  # Index corresponding to frequency f0
exponential = np.exp(-2j * np.pi * k * np.arange(len(t)) / len(t))

Why quantize to an integer?
For this use case we want a direct down-conversion from f0 to DC, i.e. np.exp(-2j * np.pi * f0 * t).
The above expression is almost equivalent, except for the int(...).

[iMo]

Here is the second source with previous settings..
Used as-is, just input data1 and data2 (removed offset).

Code: [Select]

# remove the offset
data1 = data1 - np.mean(data1)
data2 = data2 - np.mean(data2)

# Parameters
fs = xSAMPRATE  # Sampling frequency in Hz
T = 12 * XSCALE    # Signal duration in seconds
f0 = 4000  # Frequency of the sine waves in Hz

beta = 8.0  # Kaiser window beta parameter

[gf]

Here is the second source with previous settings..
Used as-is, just input data1 and data2 (removed offset).

Is the python loop over the samples really faster than np.exp() applied to the whole vector?
Unfortunately, python loops are not very fast.

[iMo]

You have to try to measure the loop'ss duration inside the DHO804
Do we have python/numpy which fits inside the android there??

Or simply redo it into C and do measure..
Below the Bessel I0=1 in the Kaiser it seems..

How to measure magnitude ratio in decibels and phase in degrees of two sine waves with help of Kaiser window and something like DFT at given frequency, while not using complex exponential function and not using costly sin and cos and atan2 functions? Do use atan2 approximation. Could you show me a C code?

Code: [Select]

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

// Constants
#define PI 3.14159265358979323846

// Kaiser Window function
void compute_kaiser_window(double* window, int length, double beta) {
    for (int n = 0; n < length; n++) {
        double term = (2.0 * n / (length - 1)) - 1.0;
        window[n] = exp(-beta * (1 - term * term));
    }
}

// Apply Kaiser window to a signal
void apply_window(double* signal, const double* window, int length) {
    for (int i = 0; i < length; i++) {
        signal[i] *= window[i];
    }
}

// Recursive trigonometric DFT computation
void compute_dft_at_frequency(const double* signal, int length, double freq, double sampling_rate, double* real_out, double* imag_out) {
    double omega = 2 * PI * freq / sampling_rate;
    double cos_omega = cos(omega);
    double sin_omega = sin(omega);

    double real = 0.0;
    double imag = 0.0;

    // Initialize recursive variables
    double cos_prev = 1.0;
    double sin_prev = 0.0;
    double cos_curr, sin_curr;

    for (int n = 0; n < length; n++) {
        real += signal[n] * cos_prev;
        imag -= signal[n] * sin_prev; // Negative for DFT

        // Recursive update of cos and sin
        cos_curr = cos_omega * cos_prev - sin_omega * sin_prev;
        sin_curr = sin_omega * cos_prev + cos_omega * sin_prev;

        cos_prev = cos_curr;
        sin_prev = sin_curr;
    }

    *real_out = real;
    *imag_out = imag;
}

// Approximate atan2
double approximate_atan2(double y, double x) {
    if (x == 0.0) {
        return (y > 0.0 ? PI / 2 : (y < 0.0 ? -PI / 2 : 0.0));
    }
    double abs_x = fabs(x);
    double abs_y = fabs(y);
    double ratio = y / (abs_x + sqrt(x * x + y * y));
    double angle = (PI / 4) * ratio;

    if (x > 0.0) return angle;
    if (y >= 0.0) return angle + PI;
    return angle - PI;
}

// Calculate magnitude and phase
void calculate_magnitude_and_phase(double real_x, double imag_x, double real_y, double imag_y, double* magnitude_db, double* phase_deg) {
    // Magnitudes
    double mag_x = sqrt(real_x * real_x + imag_x * imag_x);
    double mag_y = sqrt(real_y * real_y + imag_y * imag_y);

    // Magnitude ratio in dB
    *magnitude_db = 20 * log10(mag_y / mag_x);

    // Phase difference using atan2 approximation
    double phase_x = approximate_atan2(imag_x, real_x);
    double phase_y = approximate_atan2(imag_y, real_y);

    *phase_deg = (phase_y - phase_x) * 180.0 / PI;
}

// Main function
int main() {
    int length = 1024;
    double sampling_rate = 1000.0; // Hz
    double target_freq = 50.0;    // Hz
    double beta = 5.0;            // Kaiser beta

    // Signals
    double x[1024];
    double y[1024];
    double window[1024];

    // Generate test signals
    for (int n = 0; n < length; n++) {
        x[n] = sin(2 * PI * target_freq * n / sampling_rate);
        y[n] = 0.5 * sin(2 * PI * target_freq * n / sampling_rate + PI / 4);
    }

    // Compute Kaiser window
    compute_kaiser_window(window, length, beta);

    // Apply window to signals
    apply_window(x, window, length);
    apply_window(y, window, length);

    // Compute DFT at target frequency
    double real_x, imag_x, real_y, imag_y;
    compute_dft_at_frequency(x, length, target_freq, sampling_rate, &real_x, &imag_x);
    compute_dft_at_frequency(y, length, target_freq, sampling_rate, &real_y, &imag_y);

    // Calculate magnitude and phase
    double magnitude_db = 0.0;
    double phase_deg = 0.0;
    calculate_magnitude_and_phase(real_x, imag_x, real_y, imag_y, &magnitude_db, &phase_deg);

    // Print results
    printf("Magnitude Ratio (dB): %f\n", magnitude_db);
    printf("Phase Difference (degrees): %f\n", phase_deg);

    return 0;
}

Also

Code: [Select]

#include <stdio.h>
#include <math.h>

// Function to compute the modified Bessel function I0(x)
double bessel_I0(double x) {
    double sum = 1.0; // Initialize sum for the series
    double term = 1.0; // First term in the series
    double k = 1.0;    // Iteration counter
    double x2 = (x / 2.0) * (x / 2.0); // (x/2)^2 for repeated use

    while (term > 1e-10) { // Converge until the term is very small
        term *= x2 / (k * k);
        sum += term;
        k += 1.0;
    }

    return sum;
}

// Function to compute Kaiser window coefficients
void compute_kaiser_window(double* window, int length, double beta) {
    double denom = bessel_I0(beta); // Denominator I0(beta)

    for (int n = 0; n < length; n++) {
        double ratio = (2.0 * n / (length - 1)) - 1.0; // Calculate (2n/(N-1) - 1)
        double term = beta * sqrt(1.0 - ratio * ratio); // beta * sqrt(1 - ratio^2)
        window[n] = bessel_I0(term) / denom; // w[n] = I0(term) / I0(beta)
    }
}

// Test the Kaiser window computation
int main() {
    int length = 16; // Length of the window
    double beta = 5.0; // Shape parameter
    double window[16];

    // Compute the Kaiser window
    compute_kaiser_window(window, length, beta);

    // Print the coefficients
    printf("Kaiser Window Coefficients:\n");
    for (int n = 0; n < length; n++) {
        printf("w[%d] = %f\n", n, window[n]);
    }

    return 0;
}

[gf]

I do not own a DHO. In C, an explicit loop is required anyway, and a complex multiplication (-> iterative calculation) is very likely faster than cexp(). Drawback of an iterative calculation is that it accumulates errors, but I don't think this is a problem when double precision numbers are used. With single precision in conjunction with a huge number of samples, I would at least verify the accuracy.

EDIT:
Regarding apply_kaiser_window(): You want to (pre)calculate the window function only once, as long as length does not change.

[iMo]

Yep, added above the Kaiser coeff pre-computation with Bessel as well. Should be used in the previous C code (the Kaiser coeff calc function there seems to me buggy).
All sources above need to be verified before inserting into the DHO800/900, of course.. 

[iMo]

Orig GF math Python:
xMag dB is -93.5903211132
xPhaseA deg is 3.1959587217

AI no trigo functions Python:
Magnitude Ratio (dB): -93.5903211132
Phase Difference (degrees): -3.2231039924

AI no trigo functions C:
Magnitude Ratio (dB): -93.5903211132
Phase Difference (degrees): 1.2746656431

All above with:

Samples   600
Sampling rate   250000
Beta   8
Freq   4008.44
T   0.0024  (sig duration = 12div * 200us/div)

and the test .csv data file below - a real measurement, V_ref_inp first column, V_dut_out second column..

[gf]

Orig GF math Python:
xMag dB is -93.5903211132
xPhaseA deg is 3.1959587217

AI no trigo functions Python:
Magnitude Ratio (dB): -93.5903211132
Phase Difference (degrees): -3.2231039924

AI no trigo functions C:
Magnitude Ratio (dB): -93.5903211132
Phase Difference (degrees): 1.2746656431

All above with:

Samples   600
Sampling rate   250000
Beta   8
Freq   4008.44
T   0.0024  (sig duration = 12div * 200us/div)

and the test .csv data file below - a real measurement, V_ref_inp first column, V_dut_out second column..

To get a feeling for the estimation uncertainty, I ran a Monte Carlo analysis with 1000 simulated measurements using the above parameters and a simulated DUT output signal with a SNR of 0dB (this seems to be approximately the SNR in the 2nd column of the CSV file). And I got an estimation uncertainty (standard deviation) of ~0.7dB for the magnitude and ~3.5° for the phase. I would consider anything within +- 2 standard deviations as plausible result.

The samples in the CSV also look strange. Even the reference waveform is not clean (figure1.png) and the wiggles in the reference waveform don't look like random noise. And indeed, if I look at the power spectrum (figure2.png), I see an extra peak at 121kHz. I have no idea what the scope has done to the waveform. This is certainly not the waveform that came from the generator. Due to the frequency-selective detection, it should not matter for this use case, but it's nevertheless strange.

[iMo]

Here is the 5120 data lines long csv I got when I saved "csv data" off the scope onto an usb flash..
Time  CH1_ref  CH2_dut
Timestep in the file is 20ns
Interestingly it saved VOLTS into the csv file (via usb 8bit unsigned binary in a weird format only)
I have no idea how to download those data via usb, however..   

f0 was 301.373kHz from the generator
2us/div was timebase set on the display (aprox 7.2 periods)
50MSa/sec was Scope's sampling rate
So you have to derive the other stuff.. 

Below AI_python with beta=8 and 5120 samples..
While the magnitude looks nice with all three sources, the phase needs some finetuning, imho..

[gf]

The data look much more realistic now. Calculated manually in Octave, I get -88.430° from these data, with kaiser(8 ). The window function makes a small difference, but not that much. E.g. Kaiser(14) -> -88.433°, Blackman -> -88.432°, Blackman-Nuttall -> -88.432°, Hann -> -88.432°. The uncertainty due to the noise is much larger (see below).

The 89.5° you got is IMO not plausible, if it was calculated from the same data as in the CSV file. While measurement noise is of course not deterministic, the calculation itself (when applied to the same data) should be deterministic within numerical precision limits.

With iterative exp, the error of exp does of course accumulate, but after 5120 samples, the difference between direkt and iterative exp calculation is still below 10^-13, i.e. still negligible for the use case, and I still get -88.430° with Kaiser(8 ).

The attached image shows the expected variablilty of the estimated phase due to noise for repeated measurements (for 5120 points and a SNR of 40dB for data1, and 25dB for data2 -- i.e. similar level as in the CSV), from a Monte-Carlo simulation. Standard deviation is ~0.06°.

[iMo]

I've replaced the phase calc in the AI_python_no_trigo with the more traditional

Code: [Select]

xout = complex(X2_real, X2_imag)
xinp = complex(X1_real, X1_imag)
xgain = xout/xinp
# Compute phases
phase =  -180.0/np.pi * np.angle(xgain)

# Print results
print(f"Magnitude Ratio (dB): {magnitude_ratio_db: .10f}")
print(f"Phase Difference (degrees): {phase: .10f}")

It has to be -phase against GF_python, however.

I get now from the above BodeData3.csv

# Parameters
fs = 50000000  # Sampling frequency in Hz
T =  5120*20e-9    # Signal duration in seconds
f0 = 301373  # Frequency of the sine waves in Hz
beta = 8.0  # Kaiser window beta parameter

Magnitude Ratio (dB): -68.5645216867
Phase Difference (degrees): -88.4299209563

With above AI_C_no_trigo (the Kaiser coeffs calculated with the Bessel I0 code there)

Magnitude Ratio (dB): -68.5643563692
Phase Difference (degrees): -176.0556207288

The phase is off a bit, needs a replacement for the "atan2_aproximation" which seems to me too simple minded to be true..

PS: and this is AI_C code above with standard atan2:

Magnitude Ratio (dB): -68.5643563692
Phase Difference (degrees): -88.4303926343

Btw. - tried with a 3.57MHz crystal, frankly when not having a generator with 1Hz (or less) resolution you may not Bode_it properly ..
Mine is DDS based 1Hz resolution, so it works somehow but all the interesting stuff happens within 10-20Hz.. The contraption below shows +3dB "gain" in the peak, so parasitics must be considered somehow too..

[gf]

Btw. - tried with a 3.57MHz crystal, frankly when not having a generator with 1Hz (or less) resolution you may not Bode_it properly ..
Mine is DDS based 1Hz resolution, so it works somehow but all the interesting stuff happens within 10-20Hz.. The contraption below shows +3dB "gain" in the peak, so parasitics must be considered somehow too..

And with such a narrow-band DUT, you also need to allow enough dwell time after tuning the generator to the new frequency and before taking the measurement. I would guess 100ms or even more.

[iMo]

Not only that, the algorithm must be adaptive - it means before taking each Bode point it has to setup the measurement properly, I would not use internal Auto for that but something more specialized..

[iMo]

I finally found by chance the latest programming manual (2010, they added couple of commands after the latest fw upgrade) for my scope with additional scpi command for reading full 2x 5120 samples off my vintage DS1062CA.
Also I found a different math for decrypting the actual int8 values coming off the scope .

But.. The math with 5k real samples generates rather noisy results (not usable actually, almost random numbers). The 3 math sources above generate the same results with given sample set, so no issue with math itself, it seems. Data look ok as well..
IMHO it is caused by much higher sensitivity to errors of the DFT+Kaiser algo/math with larger N (like fs and f0 frequency alignment or jitter, spectral leakage, math errors, whatever..).

Tried with various betas with less success.. We have to find the way how to utilize that larger N actually.

So far the best results here are with 600 samples. Below an example - a Bode of the above RLC made in a loop while manually setting freq on the HP3310B, making AUTO and waiting couple of secs for stabilization. None averaging. Frequency reading by the scope (not by its "counter" however), so the freq readings are rather coarse.
Even so the 600points capture with beta 8 looks quite rugged, imho.

[gf]

But.. The math with 5k real samples generates rather noisy results (not usable actually). The 3 math sources above generate the same results so no issue with math itself it seems.
IMHO it is caused by much higher sensitivity to errors of the DFT+Kaiser algo/math with larger N (like fs and f0 frequency alignment or jitter, spectral leakage, math errors, whatever..).

Tried with various betas with less success.. We have to find the way how to utilize that larger N actually.

Hmm, IMO the 5120 point data in your BodeData3.csv seem to work pretty well with fs=50MSa/s and f0=301373.

Basically, the more points, the better. But the f0 estimate must be accurate enough so that the actual signal frequency still falls within the passband of the bandpass centered at f0. The more points, the narrower the bandwidth and the more accurate f0 must be. If you are using a digital signal generator, you should know what frequency you set the generator to (without measuring it). This value should already be a pretty good estimate of f0, suffering only from the deviation between the oscillator of the generator and the oscillator of the scope. For a reasonable generator and scope, I would not expect a deviation of more than say ±3 dozens of ppm. EDIT: E.g. the time base accuracy of the DHO800 is specified ±25 ppm. In order that the signal still falls within the filter bandwidht, make sure that the following constraint is fulfilled:

    f0_error <= fs / f0 / num_samples    or equivalently
    num_samples <= fs / f0 / f0_error    (if f0_error is given)

Example: With fs=50MSa/s, f0=300kHz and 5120 samples, f0_error must not exceed 50e6/300000/5210 = 0.032403, IOW the f0 estimate must not deviate more than +-3.2403% from the true value. Or, if f0_error is known to be less than 100ppm (0.0001), then records with a length up to 50e6/300000/0.0001 = 1.66M points can be used.

Can you post the data of a 5k in+out record and the corresponding fs and f0, of which you think that they do not work well?

Frequency reading by the scope (not by its "counter" however)

How accurate is the frequency measurement, and what was the sample rate?
With say 50MSa/s, f0=3MHz and 5000 points, the accuracy of the f0 estimate should be better than 0.3%.

[iMo]

I will create sets for you tomorrow. I spent the whole day today messing with the 5k issue..
The frequency shown by scope was off by aprox 700Hz at 3.57MHz for example against my DDS.
But I am using the "freq" off the scope reading in the math, so it should not matter so much, imho.. (the freq measurement time base and sample rate comes from the same oscillator inside the scope).
The issue is at kHz freqs too, the mag_dB jumping random within 20% off and phase is all over the place while making measurements at the "constant frequency of the HP generator" which is stable at say 4kHz within 0.2Hz during the measurements..

[iMo]

Here are the sets with 5120 and 600 samples measurements, DS1062CA and HP3310B, RLC DUT as above, xfreq=f0= aprox 3334 Hz in all measurements.

First raw in columns:
[SAMPRATE_scope, xSAMPRATE_meas, betaK, T_wide, ts, xfreq_w, xfreq_counter, xmagnitude_DB, xphase_rad, xphase_deg]

xfreq_w - by scope returned freq used in math
xfreq_counter - manually added from the scope's "counter" after the particular measurement
T_wide not used in math

the subsequent rows are:
data1 data2 raw1 raw2

The data1 and data2 were converted:
data1 = v1scale * (199.0-raw1) - 4.0 * v1scale
data2 = v2scale * (199.0-raw2) - 4.0 * v2scale
where
v1scale = 2.0/25
v2scale = 0.1/25
and offsets on the scope were set 0.0
and the offsets were removed in sw as well:
data1 = data1 - np.mean(data1)
data2 = data2 - np.mean(data2)

Time scale on display was 200us/div, aprox 8 periods visible on the scope screen, and ch1 2V/div, ch2 100mV/div, AC coupling, BW OFF, no averaging

Issues:

1. While looking at raw data I see already an issue in 5k raw data - the first 8 and last 9 bytes/samples returned from scope are wrong.. Was not visible in graphs..

2. While counting the number of periods in 5k samples it is aprox 34 (on the screen and in data aprox 8 with 600 samples). So I wonder whether the SAMPRATE and xSAMPRATE and "ts" is ok with 5k (when 5120/600=8.53)..

 

[gf]

2. While counting the number of periods in 5k samples it is aprox 34 (on the screen and in data aprox 8 with 600 samples). So I wonder whether the SAMPRATE and xSAMPRATE and "ts" is ok with 5k (when 5120/600=8.53)..

The raw samples have the original acquisition sample rate (500000) and obviously span a larger time interval than the screen width.
The 600 display samples are obviously downsampled to a lower sample rate and span the time inverval that corresponds to the screen width (implying a sample rate of 600 / screen_width_in_seconds).

Frequency estimate is OK. At 3.3kHz, 500kSa/s and 5000 points, up to 3% deviation are tolerated.

I checked ...5k_TEST{1,2,3}, and they give quite consistent results (using 500kSa/s):
-23.275dB/68.506°
-23.273dB/68.504°
-23.284dB/68.509°

Due to the garbled start and end, I discarded the first and last 60 samples, and used only the 5000 ones in the center of the record (no special reason to use exactly 5000, I just though is was a nice number).

[gf]

The raw samples have the original acquisition sample rate (500000) and obviously span a larger time interval than the screen width.

I just wonder how the 5k records are recorded at a very fast timebase (say 10ns/div), where the scope already interpolates.
Still at the (maximum) acquisition sample rate? Or as interpolated samples at a higher rate than the maximum acquisition sample rate of the scope?
I rather guess the former, but I would not rule out the latter.
[ The 600 display samples are certainly interpolated then. ]

[iMo]

Yep it uses sampling rate of the scope in the 5k output..
Below 90 measurements below resonance and in resonance (the RLC DUT above, res around 4030Hz) with 500k sampling rate and 5000samples (the first 60 and last 60 cut off).