At the maximum 3.579225 (DDS freq) the scope's counter shows 3.58032MHz.
The DDS freq used in the math and in the graph.
With 5000 point records, 300ppm deviation is OK for the detector if the sample rate is > 5.5MSa/s.
The vout/vin measurement mostly cancels out magnitude and phase errors due to a small frequency offset.
Only SNR suffers, when the filter is not centered at the frequency.
The plot still reflects the given frequency, of course.
Of course you could measure the ratio f_scope_counter/f_DDS once (at a high enough frequency) and use f_DDS*ratio for the detector calculation. Then the residual deviation is mostly a matter of the short term stability (during the sweep) of the DDS and scope oscillators.
There are small periodic wiggles on the MAG_DB visible, it may come from some effects caused by the differences in DDS generator's and scope's sampling frequencies.
Attenuator mismatch can also introduce "steps" where V/div changes (and I guess V/div was not the same for the whole sweep).
What is the frequency resolution of the used DDS (usually it is sample_rate/2^N, where N is an integer, so the resolution is usually not a decimal number). If you set a DDS to a certain frequency, it will round the given number to an integer multiple of the resolution, and the actual generated frequency will in most cases be different from the given number. Unless the resolution is significantly better than 1Hz, it is better to use the actually generated, rounded frequency for calculation and plotting.
For example: If the resolution is e.g. 0.9Hz, and you make 1Hz steps (nominal), then you actually get [0, 0.9, 1.8, 2.7, 3.6, 5.4, 6.3, 7.2,...]+f_start Hz, which is no longer evenly spaced, but every ~9th step has twice the size. Obviously, if you plot these points with even spacing at [0,1,2,3,4,5,6,7...]+f_start, you'll get a discontinuity at every ~9th frequency point.
Since your record length is limited, you can try to average the complex vout/vin quotient of N (say 10...100) measurements at the same frequency point before calculating magnitude and phase -- to reduce noise. Only the first measurement needs to do the setup, while the following N-1 measurements just capture a record and calculate the gain vout/vin. The complex sine wave exp(...) and window function do not need to be recalculated either.
Also the phase differences in coaxes, delay/skew between the channels, impedance matching, etc. might be the cause.
Thus the Quality Bode shall provide calibration too..
Of course, you can perform a response normalization by running a sweep with both probes connected to the reference signal (at the same point). This will cancel out the frequency and phase response differences between the two scope channels, including the probes. However, you are still only measuring the ratio between the voltages at the probe tips, not the full network parameters of the DUT. So the Bode plot cannot predict how the same DUT would behave in a different environment, with different source and load impedances. And it cannot take into account the loading effects of the probes. The Bode plot is not a network analyzer. The measurement setup is still potentially problematic at high frequencies.