Bode plot for scope owners that don't have the option

[shabaz]

I just took a very simple approach, since it will cover many of my use-cases without much effort. Although I have the functionality built-in too, not all 'scopes support it, and besides, sometimes it's just preferred to have the custom chart rather than a 'scope screenshot.

My simple approach was to use SCPI from Python code to control the sine-wave spot frequencies, and it's not too granular, in fact I have to list each discrete frequency in the code (not a big deal, especially since I mostly am not looking for any more detail than a few tens of points total). I made no effort to optimize the volts-per-div for dynamic range; I expect to go into the code and tweak it if needed (but usually I don't bother).

I don't use a GUI; instead, I use the Python command prompt to issue a command to run the scan, and then another command to plot the output. I can call additional functions to set the parameters for the capture, if I'm not happy with the defaults, or if I don't edit the source code.

The chart shows an audio preamp circuit I was experimenting with, with the bass turned down low in that capture.

The chart looks a bit odd without the horizontal grid lines; that was supposed to be temporary because I was experimenting with trying to align the left and right side axes so that I could have common grid lines, but was unsuccessful in figuring out how to do that in the time I had; I need to think on that more when I get around to it.

Anyway, obviously such a quick hack won't suit all, but sometimes that's good enough.

[TheoB]

..When you push plot all datasets are reduced to one dataset. That works good for two datasets (like 100Hz and 10kHz), but for sets with more overlap I have some issues with accuracy. Work in progress.
Anyhow, I placed a first version on GitHub.

The first naive approach could be - simply merge all bode points (a Bpoint is {freq, magDB, phase}) from all bode data sets into a single array/record and sort all the Bpoints based on the ascending freq, in case there are Bpoints with the same freq (it may happen) do an average of them.

When averaging duplicate Bode points of the same frequency, you should average the complex gain, not the magnitude/phase separately. And to ensure that noisy points due to a small stimulus amplitude don't ruin the result, I suggest calculating a weighted average (weighted by the stimulus amplitude at the particular frequency).

Another option is keeping the point with the largest stimulus amplitude if there are duplicates with the same frequency. However, this approach cannot improve the SNR, while averaging may enhance the SNR if there are multiple points with a large stimulus amplitude.

Another question is how to determine whether a particular FFT frequency point is considered a "Bode point" or not. If the stimulus is known to be a sum of sine waves (e.g., fundamental+harmonics), you could search for the peaks in the spectrum. For a stimulus with a continuous spectrum, or for a noise stimulus, this does not work, and the only way coming into my mind is to exclude FFT points with a stimulus power below a certain threshold, and treat all other FFT points as "Bode points". In this case, I would also prefer to average duplicate Bode points of the same frequency and not just keep the points with the largest stimulus amplitude, since there is a chance to encounter multiple acquisitions with a high stimulus amplitude, where averaging may enhance the SNR.

I averaged the complex FFT of a number of acquisitions. That can also be automated for say N=16 but I tried it first manually. That works as expected.
It's already in my code but I changed it to picking the strongest of N of CH1. The same mask (which sample to use from N) is also applied to CH2. This works best to combine acquisitions from different stimuli frequencies. Using both 100 Hz and 10kHz improves the SNR by 40dB from 10kHz onwards. Averaging with N=16 gives 12dB SNR improvement. I will add a radio button to select between these two options. Running 32 FFT over 1Msamples is time consuming though. Using two input frequencies needs only 4 FFT calculations.
The points to keep are all points that are less than 60dB below the max point. That should roughly be three decades higher than the fundamental. It's this function that determines the quality of the waves.
If all acquisitions use the same stimulus then averaging before masking is the best. If the stimuli is not the same averaging should not be used. The 100th harmonic of 100Hz might fall on the fundamental at 10kHz. Averaging the two will not improve the SNR as the amplitude levels are vastly different.

I apply a Hanning window over the data at the cost of some SNR. That's more needed for plotting single channel FFT. The same code is reused for the Bode plot. I now realize this is perhaps not needed. The frequency bleeding will be the same for both channels for which I plot the difference.

[gf]

I averaged the complex FFT of a number of acquisitions.

You mean averaging fft(CH1) and fft(CH2) separately?

FFT is a linear operation (multiplication of the time domain signal vector of length N with a huge complex NxN matrix), and average is basically just a sum and division by a constant, so the distributive law applies. Therefore average(fft(CH1)) is equivalent to fft(average(CH1)). So what you get is the fft of the average of the CH1 traces. [ Btw, this also offers a cheaper alternative to calculate this average: Average the time domain traces and calculate the FFT only once. ]

Note that this kind of averaging will only give you the desired result if the traces to be averaged are in phase for each frequency contained in the signal. In the end, you still rely on the trigger to do this alignment. But this does not work in general. Of course, it is supposed to work if the same periodic waveform is captured multiple times. But it will not do what you want it to do if you average two traces of, say, a 25% duty cycle square wave, where one trace captures a 10 Hz tone and the second trace captures a 100 Hz tone, both triggered on the rising edge. Then the 10th harmonic of the 10 Hz tone in the first trace is not in phase with the fundamental of the 100 Hz tone in the 2nd trace (-pi/4 vs. -pi/2). Therefore, they do not qualify for this kind of averaging.

What I actually had in mind instead is average(fft(CH2)/fft(CH1)), i.e. averaging the complex gains of multiple acquisitions. The expected value for the phase angles of this quotient is the phase of the DUT's transfer function, which is invariant with respect to the acquisition phase. Therefore, this average does not depend on the trigger, but relies only on the coherence between CH1 and CH2, which is achieved by sampling CH1 and CH2 simultaneously.

I apply a Hanning window over the data at the cost of some SNR. That's more needed for plotting single channel FFT. The same code is reused for the Bode plot. I now realize this is perhaps not needed. The frequency bleeding will be the same for both channels for which I plot the difference.

There are two main considerations:

1) Scalloping loss: You are correct, it cancels out in the complex gain fft(CH2) / fft(CH1), but of course it also reduces the SNR because a lower signal level is detected by the FFT bin. OTOH, for displaying the spectrum of a periodic signal you usually want to avoid scalloping loss by using a flattop window in order to show the amplitudes of the harmonics correctly, even for arbitrary signal frequencies.

2) Bleeding between harmonics does not cancel out, and therefore also matters for the complex gain. Especially for a narrowband DUT, you likely want a window function with a high selectivity and high stopband attenuation (maybe even 100+ dB), like e.g. Blackman-Harris, Blackman-Nuttall, or Kaiser with an appropriate beta.

EDIT:

Plotting fft(CH1) and fft(CH2) is, of course, only for your diagnosis.
It's not part of the actual Bode plot.

EDIT:

For better understanding, here is pseudo code for what I had in mind:

Code: [Select]

// initialize
for f in 0...FFT_SIZE-1 {
    gain[f] = 0 + 0i; // complex
    weight_sum[f] = 0;
}
for each acquisition {
    Vin = fft(CH1);
    Vout = fft(CH2);
    for f in 0...FFT_SIZE-1 {
        weight = abs(Vin[f]); // magnitude of reference signal at frequency f
        // skip frequencies where the drive level is too low
        if (weight > threshold) {
            gain[f] += Vout[f] / Vin[f] * weight;
            weight_sum[f] += weight;
        }
    }
}
for f in 0...FFT_SIZE-1 {
    if (weight_sum[f] > 0)
        gain[f] /= weight_sum[f];
}
discard/ignore all points gain[f] where weight_sum[f] == 0
plot magnitude and phase of the remaining points.
[code]


[TheoB]

Yes, good point averaging can just as well be done by the scope in the time domain 
I try to work on arrays and not loop over all frequency points. The numpy library is optimized for working on large arrays.

Code: [Select]

        if mode=='fftmean':
            fft_val1=np.mean(fft1,axis=0)
            fft_val2=np.mean(fft2,axis=0)
fft1 and fft2 are two dimensional arrays of FFT results for CH1 and CH2. But I use the max of each now, not the mean.

Code: [Select]

        if mode=='fftmax':
            fft_val1=np.max(fft1,axis=0)
            fft_val2=np.max(fft2,axis=0)
next  I select samples which are not smaller that 60dB from the peak value.

Code: [Select]

      xbode = Vrms1 > (max(Vrms1[2:]) / 1000)
This xbode is a selector (array of booleans) that select only the samples I want. Samples near DC I also filter out.

Code: [Select]

        # Use range input from user
        xbode[xf < float(self.bodefmin.get())] = False
        xbode[xf > float(self.bodefmax.get())] = False

        bodedb = 20 * np.log10(Vrms2[xbode]) - 20 * np.log10(Vrms1[xbode])
        PhaseDiff=np.angle(fft_val1[:data["N"] // 2])-np.angle(fft_val2[:data["N"] // 2])
        bodephase = np.unwrap(PhaseDiff[xbode], period=360)
Your pseudo code looks similar except that you still average if more than one acquisition contains more energy than the threshold. I don't see the practical use for that. Averaging is done by the scope.

[gf]

Code: [Select]

        if mode=='fftmax':
            fft_val1=np.max(fft1,axis=0)
            fft_val2=np.max(fft2,axis=0)


Be careful with np.max(). For complex numbers, it seems not to do what you want - it does not select the largest magnitude:

Code: [Select]

>>> np.max([0+2j,1+0j])
(1+0j)

Sort order for complex in python seems to be not by magnitude, but lexical, using the real component as primary key and the imaginary component as secondary key.

[tggzzz]

Code: [Select]

        if mode=='fftmax':
            fft_val1=np.max(fft1,axis=0)
            fft_val2=np.max(fft2,axis=0)


Be careful with np.max(). For complex numbers, it seems not to do what you want - it does not select the largest magnitude:

Code: [Select]

>>> np.max([0+2j,1+0j])
(1+0j)

Sort order for complex in python seems to be not by magnitude, but lexical, using the real component as primary key and the imaginary component as secondary key.

Oh, wow.

I'd like to know the definitive root cause of that "surprising" result, since it might affect far more than just that one function.

[iMo]

My dsp math gets a bit rusty with the decades so let me ask following:
1. I do bode1 with 1000Hz fundamental (any starting phase, I do for example Kaiser_beta_8)
2. my data set is bode1[freq, real, imag]
3. I do bode2 with 1200Hz fundamental (any starting phase, I do for example Kaiser_beta_8)
4. my data set is bode2[freq, real, imag]
5. I do bode3 with 2000Hz fundamental (any starting phase, I do for example Kaiser_beta_8)
6. my data set is bode3[freq, real, imag]
7. I do sort bode1 and bode2 and bode3 ascending based on freq
8. I can see the freq=6000 is common for bode1_6th and bode2_5th and bode3_3rd harmonic
9. How do I make the average of the three points?

[RoGeorge]

There are drawbacks when using a signal with a rich spectrum to deduce the Bode plot (thinking here using noise as an input signal, or a frequency comb, or a Dirac, etc.).  By drawbacks I mean limiting things like these:
- low number of bits (8..12bit) in an oscilloscope ADC
- low dynamic range of an oscilloscope analog stage (can't use a big input pulse(s), might saturate the input of an oscilloscope, and the effect of the short impulse(s) might be too low in amplitude for the ADC)
- intermodulation distortions (because of many frequencies are present at once, and interact with each other)

Looks like would be easier to measure if we have only one frequency at a time.  However, this will take a lot of measurements, and the measuring loop has to wait for the signal to stabilise at each frequency step, so repeatedly measuring the amplitude at different frequency might take too long.


Here's an idea:
- we apply a continuous frequency sweep at the input, at a constant amplitude, and a known start/stop frequency range
- sample the entire duration of the frequency sweep in single mode trigger (my 10 years old DSO, Rigol DS1054Z can take up to 24 million consecutive ADC samples, then transfer them later to a PC)
- on the PC, apply the Wavelet Transform and plot the spectrum over time

The problem is that I'm not very confident with my lame math skills, so I'm not very sure if the Wavelet Transform idea will work in practice.  My understanding when I've looked into the Wavelet Transform for the first time (https://www.eevblog.com/forum/chat/fun-for-nerds/msg3554257/#msg3554257), and very first application I thought of after reading the PDFs, was to use a continuous frequency sweep, sample, then apply Wavelet Transform, and plot the output, which will be the spectrum over time.  (There are already Python modules that can calculate the WT.)

Since we know the frequency of the input signal at every moment during the frequency sweep, we can apply WT on the output of the DUT, and the projection of the WT on the amplitude plane would be the Bode plot. 

This looks like some sort of "Captain Obvious" type of application for the Wavelet Transform.
Is this a known technique, to produce a Bode plot by using a frequency sweep then WT? 

[TheoB]

Code: [Select]

        if mode=='fftmax':
            fft_val1=np.max(fft1,axis=0)
            fft_val2=np.max(fft2,axis=0)


Be careful with np.max(). For complex numbers, it seems not to do what you want - it does not select the largest magnitude:

Code: [Select]

>>> np.max([0+2j,1+0j])
(1+0j)

Sort order for complex in python seems to be not by magnitude, but lexical, using the real component as primary key and the imaginary component as secondary key.

Very good point. This is a bug due to my wrong assumption that max considers magnitude. I was already investigating some wrong results I sometimes got with multiple captures. Will fix.

[TheoB]

My dsp math gets a bit rusty with the decades so let me ask following:
1. I do bode1 with 1000Hz fundamental (any starting phase, I do for example Kaiser_beta_8)
2. my data set is bode1[freq, real, imag]
3. I do bode2 with 1200Hz fundamental (any starting phase, I do for example Kaiser_beta_8)
4. my data set is bode2[freq, real, imag]
5. I do bode3 with 2000Hz fundamental (any starting phase, I do for example Kaiser_beta_8)
6. my data set is bode3[freq, real, imag]
7. I do sort bode1 and bode2 and bode3 ascending based on freq
8. I can see the freq=6000 is common for bode1_6th and bode2_5th and bode3_3rd harmonic
9. How do I make the average of the three points?

In this example you have three results for the same frequency (bin).  If the amplitude would be the same, averaging would be easy, giving the same weight to all three observations. Here the amplitude ratio is 2:1.5:1. For simplicity assume the amplitude ratio is two using two samples. This means the noise in the weak sample is twice as high as the noise in the other sample. Both samples represent the same CH2/CH1 ratio. Goal is to lower the noise by averaging. Using equal weight to both samples would decrease the SNR by 1dB where we normally expect an increase of 3dB. A weighted average would still increase the SNR but by only 1dB. Averaging  is only efficient if the signal amplitude is about the same.

[TheoB]

There are drawbacks when using a signal with a rich spectrum to deduce the Bode plot (thinking here using noise as an input signal, or a frequency comb, or a Dirac, etc.).  By drawbacks I mean limiting things like these:
- low number of bits (8..12bit) in an oscilloscope ADC
- low dynamic range of an oscilloscope analog stage (can't use a big input pulse(s), might saturate the input of an oscilloscope, and the effect of the short impulse(s) might be too low in amplitude for the ADC)
- intermodulation distortions (because of many frequencies are present at once, and interact with each other)

Looks like would be easier to measure if we have only one frequency at a time.  However, this will take a lot of measurements, and the measuring loop has to wait for the signal to stabilise at each frequency step, so repeatedly measuring the amplitude at different frequency might take too long.


Here's an idea:
- we apply a continuous frequency sweep at the input, at a constant amplitude, and a known start/stop frequency range
- sample the entire duration of the frequency sweep in single mode trigger (my 10 years old DSO, Rigol DS1054Z can take up to 24 million consecutive ADC samples, then transfer them later to a PC)
- on the PC, apply the Wavelet Transform and plot the spectrum over time

The problem is that I'm not very confident with my lame math skills, so I'm not very sure if the Wavelet Transform idea will work in practice.  My understanding when I've looked into the Wavelet Transform for the first time (https://www.eevblog.com/forum/chat/fun-for-nerds/msg3554257/#msg3554257), and very first application I thought of after reading the PDFs, was to use a continuous frequency sweep, sample, then apply Wavelet Transform, and plot the output, which will be the spectrum over time.  (There are already Python modules that can calculate the WT.)

Since we know the frequency of the input signal at every moment during the frequency sweep, we can apply WT on the output of the DUT, and the projection of the WT on the amplitude plane would be the Bode plot. 

This looks like some sort of "Captain Obvious" type of application for the Wavelet Transform.
Is this a known technique, to produce a Bode plot by using a frequency sweep then WT? 

It's fully understood that is not the best Bode plot implementation. Where it wins is in price performance ratio. I find 100dB dynamic range reasonable when measuring bandpass filters. Limited by the scope's dynamic range. That can only be increased by changing the scope frontend to give more gain when a weak out of band signal needs to be measured. And then you must use a sine signal.
The WaveLet idea gives time information when a chirp signal is applied. That's not needed for this application. Just capture an entire sweep in the time domain and consider it to be stationary (it repeats forever) and perform FFT. I don't expect much better performance from a chirp compared to a constant frequency signal. The generator will only very shortly output a tone for each individual frequency bin. It needs to sweep linearly as the FFT bins are all the same size. The sweep time and capture time need to be carefully aligned. If you limit the sweep range to that of a crystal resonator, you do have a big advantage as all energy is concentrated around the frequency of interest. An FM modulated signal should work as well.
Why not try it out 

[Smokey]

Bode plot from swept sine (with matlab code):
https://www.eevblog.com/forum/projects/bode-plot-from-swept-sine-input-and-scope-captured-data/
https://pages.hmc.edu/mspencer/e151/AutoBodeGuide.pdf

Another paper describing the same thing:
https://www.thinksrs.com/downloads/PDFs/ApplicationNotes/SR1_SweptSine.pdf

[gf]

Yes, good point averaging can just as well be done by the scope in the time domain 

... if the scope returns the averaged trace at a resolution higher than the ADC resolution, and if the scope still supports full memory depth for the averaged trace.

[TheoB]

Spend my weekend to fix two bugs as some things were actually
I updated dho-remote to v1.1

• The best sample from multiple captures is correctly chosen.
• Phase display is displayed in degrees and not in radians.

Added an input field to specify the amplitude of the harmonics to be used relative to 0dB full scale. That is needed if CH2 has a much smaller signal like when measuring a bandpass filter or crystal.

Yes, good point averaging can just as well be done by the scope in the time domain 

... if the scope returns the averaged trace at a resolution higher than the ADC resolution, and if the scope still supports full memory depth for the averaged trace.

The waveform data of the Rigol is 16 bits and that is used for averaged data. Although it improves the SNR, it also limits the memory depth to 1M  . For higher memory depth the averaging could be done in software.

I also experimented with a sine wave as stimulus. That gives one point per capture. A complete Bode diagram can be build in this way as well but doing so is a lot of work.

[gf]

I also experimented with a sine wave as stimulus. That gives one point per capture. A complete Bode diagram can be build in this way as well but doing so is a lot of work.

Well, if the generator is controlled from the program, then it's not much work. But it takes quite a while...

Since a chirp was also proposed in this thread, I tried to simulate a measurement of a 4th order Butterworth notch filter with -3dB notch width of 3.995...4.005 Mhz (similar to yours, but with higher Q) with a single chirp impulse sweeping from 3.97...4.03 MHz. Measurement window is 2000000 samples @625 MSa/s, chirp impulse duration is 1600µs (1000000 samples @625 MSa/s). No averaging. I also included random noise of 0.1% RMS of full scale (SNR -51 dB_FS) for each of the 2 channels, and 12-bit quantization in the simulation. The resulting plots (see attachment) look stunningly clean Expected dynamic range is ~85 dB. The processing gain (and therefore the expected DR) does not depend on the number of samples, but on the sample_rate/sweep_span ratio, i.e. for high DR we want a high sample rate, and a narrow span.

Just capture an entire sweep in the time domain and consider it to be stationary (it repeats forever) and perform FFT.
...
The sweep time and capture time need to be carefully aligned.

Just treat the chirp impulse as what it is, a single impulse with finite support. For a time-limited discrete-time signal (which is essentially zero outside the support region), the DFT calculates discrete-frequency samples of the signal's continuous DTFT spectrum. That's exactly what we want. Don't use a window function, a time-limited impulse is already (self-)windowed.

There is no need for exact alignment of the measurement window. Just make sure that both the support of the chirp impulse and also the support of the DUT's response to the chirp impulse fall completely within the measurement window. Note that the DUT may delay and/or stretch the stimulus impulse, so the measurement window must be wider than the support of the chirp impulse to capture the full DUT response until it has decayed/settled essentially to zero after the end of the chirp impulse. Note that the settling time of a DUT with narrow bandwidth can be quite long. It doesn't matter if you include some additional leading and trailing "silence" (signal=0) in the measurement window before and after the impulses (but don't include more than necessary because the captured "silence" will still contain noise).

EDIT: If you want to consider it stationary, then let the signal generator repeat the chirp periodically and proceed as you would with any periodic signal: Wait for the DUT to settle, then capture a couple of periods of the stationary state and apply a decent window function. Of course, unlike for a single pulse, the spectrum is not continuous then, but a comb formed by the fundamental+harmonics of the repetition rate.

EDIT: Corrected scale of phase plot. Corrected chirp and window lenght (missing "0").

[TheoB]

I cannot generate a single chirp with support. Otherwise I would have experimented with it.
Instead I tried a 4MHz FM modulated carrier (10kHz deviation, 625Hz modulation). That's a repeating signal, so I captured exactly one modulation period (Samplerate/Mdepth=625e6/1e6=625Hz).
The best and easiest result I still get by capturing a single rising or falling edge of a square wave signal.

[iMo]

I've been playing with my older python Bode experiments stuff on the 804 via usb, also tried with yours from github, but it works only from 3.10 upwards because of "match case".

Btw., the built in FFT in the scope - I can see a difference in reporting the amplitude of a sine fundamantel (in Vrms) for different windowing functions - like from 620mVrms to 707mVrms (at 50kHz, 10dBm on 50 ohm). Shouldn't we see identical values? They have to apply correction coefficients for each window, haven't they?

[gf]

Btw., the built in FFT in the scope - I can see a difference in reporting the amplitude of a sine fundamantel (in Vrms) for different windowing functions - like from 620mVrms to 707mVrms (at 50kHz, 10dBm on 50 ohm). Shouldn't we see identical values? They have to apply coefficients for each window, haven't they?

Only if the frequency of your sine wave is exactly one of sample_rate / num_fft_points * k (where k is an integer in the range 0...num_fft_points/2), then any window function will give the same result. For other frequencies, you will encounter more or less scalloping loss with different window functions. Use a flattop window to minimize scalloping loss, i.e., if you want to accurately measure the amplitude of sine waves of arbitrary frequency.

EDIT: Think of the windowing + FFT as a filter bank of num_fft_points bandpass filters, whose center frequencies have a spacing of sample_rate / num_fft_points, and a power and phase detector attached to the output of each filter. If you plot the (overlapping) frequency responses of these bandpass filters, you can see the magnitude response of each bandpass filter responds to a given signal frequency. Only if the shape of the frequency response of the bandpass filters has a flat top with a width >= sample_rate / num_fft_points, then at least one of the filters (FFT bins) can respond with full amplitude to a sine wave of arbitrary frequency.

[iMo]

Ok, with 31.25Msa/s and 1MPts, 93.750kHz sine (31.25*3000), DDS gen at 50ohm, FFT set center 93.750kHz:

Rectangular     242mVrms
Blackman        247
Hanning          247
Hamming        247
Flattop            257 (too wide?)
Triangle           232

[iMo]

Ok, now with python 3.11.9 (pyvisa-py, numpy, matplotlib, scipy, tkinter, tkinter.ttk installed)

Code: [Select]

C:\Users\Smith\MyCode\Python\dho-remote-main>python dho-remote.py 192.168.0.66
  File "C:\Users\Smith\MyCode\Python\dho-remote-main\dho-remote.py", line 373
    print(f"Sample rate={data["sr"] / 1e6:g} Ms")
                               ^^
SyntaxError: f-string: unmatched '['
It seems it does not like " ", works with ' '

[iMo]

With my 4kHz LC res circuit fed via 33k resistor (200Hz fundamental, sawtooth, HP3310B, 50ohm out, 500kSa/s, 10kpts, 10 precaptures)..
Nice!

PS: the channel colors in the graph do not match the scope channel's colors, it seems..

PPS: v2 With my 4kHz LC res circuit fed via 33k resistor (200Hz fundamental, sawtooth, HP3310B, 50ohm out, 31.25MSa/s, 1Mpts, 16 precaptures?)..

[gf]

Ok, with 31.25Msa/s and 1MPts, 93.750kHz sine (31.25*3000), DDS gen at 50ohm, FFT set center 93.750kHz:

Rectangular     242mVrms
Blackman        247
Hanning          247
Hamming        247
Flattop            257 (too wide?)
Triangle           232

1) Is the number of FFT points really 1,000,000, or 1,048,576 (2^20, power of two)?

2) Is the frequency accurate enough? (possible deviation between generator and scope time base?)
For 31Hz bin spacing, the deviation should not be much more than 1 Hz to consider it "essentially exact".

3) I remember older posts on EEVblog saying that the flattop window was broken in the DHO scopes. I don't have a DHO and I don't know if this applies to all DHO models and if it has been possibly fixed in the meantime. Perhaps one of the owners can say more about this. According to the cockroach theory, if one window function is broken, I would not rule out that others are broken, too.

[ If you calculate Windowing + FFT yourself from time domain samples, make sure that the window function is normalized such that mean(window_function)==1 in order to get comparable amplitudes to a FFT w/o window function (for a signal frequency matching exactly one of the bin frequencies). ]

[gf]

PPS: v2 With my 4kHz LC res circuit fed via 33k resistor (200Hz fundamental, sawtooth, HP3310B, 50ohm out, 31.25MSa/s, 1Mpts, 16 precaptures?)..

@TheoB, there still seems to be something wrong with the stairsteps in iMo's Bode plot, which are particularly visible in the 300...3000 Hz range. They do not exist in the real transfer function, and they are not noise either. IMO they are artifacts from the calculation.

I may be wrong, but I suspect that these stairsteps occur because the set of "Bode points" includes (many) frequencies that are not present in the stimulus signal. The stimulus signal is periodic, and a periodic signal does not have a continuous spectrum, but contains only the fundamental + harmonics, and no frequencies in between. Any leakage from these individual frequencies into neighboring bins in the FFT spectrum is not real - it does not exist in the original spectrum of the signal. Therefore, only the FFT bin closest to each fundamental or harmonic frequency should be considered a "Bode point" (and even that is still just an approximation, since the true frequency may lie somewhere between two adjacent bins). If the stimulus spectrum is sparse, then only sparse points of the transfer function can be measured. Only if the stimulus spectrum is continuous (-> non-periodic signal), then every FFT point is a potential "Bode point" candidate.

[iMo]

1) Is the number of FFT points really 1,000,000, or 1,048,576 (2^20, power of two)?

2) Is the frequency accurate enough? (possible deviation between generator and scope time base?)
For 31Hz bin spacing, the deviation should not be much more than 1 Hz to consider it "essentially exact".

3) I remember older posts on EEVblog saying that the flattop window was broken in the DHO scopes.

1) The user guide says "1Mpts" max..
2) The deviation could be way much larger than 1Hz, sure. Let say +/-10Hz at that 93.750KHz (50ppm crystal in the DDS and 25ppm in the scope afaik) ..
3) Yup, there is a patch file with the new window, afaik, not found it and not installed it yet.

[gf]

1) Is the number of FFT points really 1,000,000, or 1,048,576 (2^20, power of two)?

2) Is the frequency accurate enough? (possible deviation between generator and scope time base?)
For 31Hz bin spacing, the deviation should not be much more than 1 Hz to consider it "essentially exact".

3) I remember older posts on EEVblog saying that the flattop window was broken in the DHO scopes.

1) The user guide says "1Mpts" max..

It's not always clear what "1M" is intended to mean. Strictly speaking, the prefix "Mega" (M) means 1,000,000 and "Mebi" (Mi) means 1,048,576, but some people don't take that too seriously. Some scopes only support power of 2 FFT sizes (radix-2 FFT). If the Rigol DHO supports decimal FFT sizes (mixed radix-2 and radix-5 FFT), then this would certainly be convenient.

2) The deviation could be way much larger than 1Hz, sure. Let say +/-10Hz at that 93.750KHz (50ppm crystal in the DDS and 25ppm in the scope afaik) ..

What does the scope's frequency counter measure? Adjust the frequency until the scope sees the desired frequency when you perform this testcase.