Extract precise amplitude and phase from a frequency sweep (VNA from DSO+AWG)

[DiTBho]

I have no idea on is how the sweeping of frequency will upset things.

step transient

you have it even from a simple differential equation applied to a feedback (G_loop) system with one pole, and not in the origin ...

Such a system can triggers permanent oscillations as well as it can saturate as well as it can dampens early oscillations within a "step transient" time.

When the step transient time is unknown, and the system is assumed to have a nice transient response, I think you can monitor data, looking for stable results. If any.

[gf]

I have no idea on is how the sweeping of frequency will upset things.

A chirp of finite duration certainly has discontinuities at the beginning and the end. These discontinuities need to be tapered, for instance by applying a tukey window, having the consequence that only the center (say 99%) of the chirp is usable at the end. However, aside from the problem at the endpoints, a (linear) chirp is a spread spectrum signal, having a nice flat amplitude spectrum, which contains no discrete frequency lines.

But a spread spectrum stimulus is also a disadvantage, because a narrow bandpass filter cannot achieve processing gain. If the bandwidth is reduced, the filter captures less noise power, but also less signal power, and the SNR remains the same.

OTOH, if the stimulus signal energy is concentrated at a single frequency line (or at a limited number of discrete frequency lines), then the signal power passing through a bandpass centered at one of the frequency lines is independent of bandwidth, while the noise power passing through the bandpass still decreases with decreasing bandwidth, therefore the SNR can be improved by reducing the bandwidth, i.e. we can obtain processing gain by capturing a larger number of samples.

[Marco]

- the unknown circuit to measure can introduce phase shifts, changes in amplitude, and can produce harmonics

Unless the circuit is very fast and within the frequencies of interest only has a simple non-linear transfer curve without significant memory, you need to do something like Volterra series estimation to characterise it (if it is fast enough to have no significant memory, you just need a curve tracer, not this). At best this measurement will tell you what the output is for the exact chirp at the exact amplitude you applied and nothing much else.

Worrying about windowing for characterising a nonlinear circuit is also a bit silly, get all that linear thinking out of here.

[RoGeorge]

Did a brief read about the Hilbert Transform (HT) this Sunday, and it looks a lot like what I was heaving in mind.  Now I'm not sure if I should be embarrassed for not knowing what the HT is, or proud for rediscovering the wheel. 

Anyways a bummer, the fun was much bigger reinventing the wheel than learning somebody else's math. 

Just kidding, I'm grateful for being pointed to the HT, thank you all for that.  The math alone seems very intimidating to me, while in practice it looks trivially simple, which makes me think I still don't know what HT does, or what is it for. 

As I was saying, still not sure what HT is, but assuming HT only introduces (two?) 90* phase shift(s) for each frequency present in a signal, then what I was proposing was (probably) using only the real part of the HT from the measured signal.  Not sure.  Also, most of the remarks, advice and example posted here are not yet in my understanding range.

Tried to get the ADC samples from the DSO into PC, to test my idea on real data.  Had to dig into 100+ forks for Python-IVI drivers, and merge the relevant ones.  Then it still didn't work and had to debug it with Wireshark, and so on.  Still working on it, and already have communication with both the AWG and the DSO. 

Meanwhile, did a more realistic simulation in LTspice.  This time, emulated an AWG with a logarithmic chirp, then emulated a DSO, then plot the frequency response for two different DUT, an RC low pass, and an RLC circuit.  The observed amplitude Bode plots correspond with the theoretical response of the DUT.

This one is for RC LP:


Please note that no reference chirp was sampled, and only the measured signal is used to calculate (on-the-fly) the amplitude response of the DUT (the red trace).
 
The chirp is done with logarithmic frequency variation, though it chirps from high to low F, to minimize any errors caused by an eventual phase noise (no noise was added into simulation for now).

The next one is for a resonant RLC:


There are some artifacts at the ends of the spectrum
- some errors at the beginning of the chirp are caused because there is no data previous to t0.  These can be easily cut out
- other type of errors are caused when the frequency variation of the sweep is too fast, observed at the low frequencies end of the chirp.  Same as before, the interval can be cut away, or the frequency can be swiped slower in the reference chirp, so the frequency will appear as quasi-constant along a quarter of a wavelength.
- the most amplitude wiggle seen in the red trace was because of sweeping too fast, not because of ringing

[RoGeorge]

Adding this link just for the docs, because it tells what the Hilbert Transform is/does in a meaningful way:
https://www.comm.utoronto.ca/frank/notes/hilbert.pdf

I think I finally got it what the HT is about. 

What I was doing by reusing the same (measured) signal with a pi/2 wavelength lag, is the same as taking the HT at that particular wavelength only.  The pi/2 lag trick will work only assuming all the following conditions are true:

• 1. that particular wavelength is known (the frequency of the reference signal is known)
• 2. there are no other frequencies present in the measured signal (the reference signal was considered monochromatic and the DUT was considered LTI, so there is only one frequency at any given time along a slow chirp)
• 3. the frequency is constant (it is not, so instead, the reference frequency sweep has to be made slow enough so that the change in frequency is negligible when looking at pi/2 behind)
• 4. the amplitude is constant (it is not, same mitigation as for point 3, sweep slow enough)

That's a lot of assumptions/limitations there.
The advantages would be:

• computationally cheap when compared with FFT (best known FFT is O(n*log(n)), while this one is O(n), which might be a big advantage when the signal is of many millions of ADC samples)
• the time lag is very small (pi/2), can be computed on-the-fly, big advantage if the reaction time is critical
• doesn't need to know the reference chirp, it only needs to know the variation law of the reference chirp, plus a start-of-chirp marker
• if the lag is not critical, then the DC results can be further averaged to filter out noises, also on-the-fly moving average is possible

[hubi]

It took me a while to make sense of the plots until I figured out the chirp is from high to low frequency. Are you sure you used a pi/4 delay?

If you want to use a (time domain convolution) Hilbert filter for this, the maximum number of samples per period you intend to use will be important. If it is above, say 80 samples, generating a reasonably compact filter kernel is not going to be possible.

The optimized kernel I posted earlier can actually be improved a bit, I had an unnecessary constraint in the filter definition. But the size was already pretty big in order to match hamster's kernel. Processing the 8M data points from the example takes about 4 seconds on my very slow laptop without any optimization.

What programming language do you intend to use for this? If it's c++, I can help out with some code and optimized kernels for Hilbert filters.

[RoGeorge]

It took me a while to make sense of the plots until I figured out the chirp is from high to low frequency. Are you sure you used a pi/4 delay?

If you want to use a (time domain convolution) Hilbert filter for this, the maximum number of samples per period you intend to use will be important.

I was thinking all the time about "quarter of the circle" while writing everywhere pi/4 instead of pi/2... 

It's not that I didn't know the full circle is 2*pi.  Now I have to revisit all the code and simulations, to check if I didn't make the same "brain typo" of writing pi/4 (45°) while thinking about pi/2 (90°).   

Thank you very much for pointing me to that mistake, I wouldn't have spot it myself too soon.



The language of preference would be Python (with its 'numpy' module for math, which numpy is a wrapper for C/FORTRAN math libs, as far as I've heard).  Python is easier to use for a non programmer like me.  Could use C, too, if needed.  Can't say I know C++, but I might be able to use it code-monkey style, with stackoverflow googling all night long. 

So far still trying to figure which way to go, FFT or with this HT-related time processing, where a copy of the same signal at a distance of 90° in the past might be good enough, instead of a proper sin cos quadrature.

I'm not sure why I like this sin²(x)+cos²(x)=1 trick so much, maybe I should just use FFT.  Playing around with alternatives mostly because it's a hobby project, without fixed performance specs to met or deadlines, IDK.

Got some (buggy) IVI drivers for Python-IVI, so I hope these days to get ADC samples from a physical DUT (using the oscilloscope as an ADC, a Rigol DS1054Z DSO).  I remember you were asking if there is any data set to download, so to evaluate the performance for this method of yours described here:  https://www.eevblog.com/forum/programming/extract-precise-amplitude-and-phase-from-a-frequency-sweep-(vna-from-dsoawg)/msg4504627/#msg4504627 Will let you know as soon as I got the Python-IVI working.  If the IVI drivers won't work, will give up to the compatibility idea (between various models of instruments), and use low level SCPI for Rigol DS1000Z oscilloscopes only).

[nctnico]

If you are using FFT, you might as well send a burst of noise into the DUT and subtract the FFT of the incoming signal from the resulting signal.

Where it comes to talking to instruments using Python, I can highly recommend using Pyvisa with the internal Visa drivers.

[hubi]

The plots look like you used pi/2 where it counts and just mislabeled it.

You could implement the Hilbert transform with the FFT as well and compare the results. That is easy to do in python. I'm not very familiar with python signal processing libraries, I'll take a look and see if there is anything available for this. I'm sure there are convolution filters and you can always include the zeros in the filter mask if you're not concerned about efficiency, but I'm not so sure about the recursive filter and pyramids.

What is the maximum number of samples per period you intend to use?

[mawyatt]

Interesting thread, brings back memory's from ~80 when we tasked with implementing the Chirped Z Transform (CZT) based real time handheld battery powered SA.

The Chirps were created by I and Q custom multiplying DACs where the scaled input waveform was the DAC Reference and the Chirp coefficients were DAC digital words updated from SRAM. This created a Discrete Time Continuous Amplitude (DTCA) Chirp I & Q multiply result which was passed onto I & Q Convolvers created by very long custom Bucket Brigade CCD devices also operating in DTCA fashion. This program was awarded to us after two well known advanced technology companies failed, and like them we had in-house custom IC design and fabrication capability which was necessary with the multiple custom chips involved. The program was highly successful and put us "On the Map" in the advanced technology arena for future programs.

Since all the "computations" were done in custom hardware, the end result back then was very efficient from a SWAP standpoint.

Surpised the CZT hasn't found more appeal as its seems to be a hidden gem within the DSP domains, and certainly interesting to see it discussed here!!

Best,

[adeuring]

The language of preference would be Python (with its 'numpy' module for math, which numpy is a wrapper for C/FORTRAN math libs, as far as I've heard).  Python is easier to use for a non programmer like me.  Could use C, too, if needed.  Can't say I know C++, but I might be able to use it code-monkey style, with stackoverflow googling all night long. 

If you are considering to work with Python and NumPy, I'd strongly recommend to also look at SciPy, a huge library based on NumPy. Especially its "signal" sub-package provides a lot of functions you often use for signal processing.   https://docs.scipy.org/doc/scipy/reference/signal.html shows you what is available.

[RoGeorge]

What is the maximum number of samples per period you intend to use?

- The oscilloscope I have can not stream the ADC data, it needs to store a trace in its memory first.
- It can capture maximum 24 million ADC samples of 8 bits.  The maximum sample rate is 1GSa/s.
- The oscilloscope analog bandwidth is 100MHz.  This means 10 samples per period at 100MHz in.

In case of too few points/period, or when the 1/4 wavelength delay falls at a moment between two samples, I was thinking of doing interpolation, linear at first, to keep it simple.

Chirped Z Transform (CZT)

Took a brief look at CZT yesterday, very interesting, thank you.

[DiTBho]

USB DSOs look more interesting here.
One bulk command, and they return samples.

my old RIGOL DSO doesn't have the ethernet port, it has a nice RS232 port.
I can ask the DSO to export samples.

[DiTBho]

There are DDS chips able to generate A*sin(wt), A = const, w = f(freq ), freq = { 100Hz ... 20Mhz }.
SPI protocol.

[hubi]

What is the maximum number of samples per period you intend to use?

- The oscilloscope I have can not stream the ADC data, it needs to store a trace in its memory first.
- It can capture maximum 24 million ADC samples of 8 bits.  The maximum sample rate is 1GSa/s.
- The oscilloscope analog bandwidth is 100MHz.  This means 10 samples per period at 100MHz in.

In case of too few points/period, or when the 1/4 wavelength delay falls at a moment between two samples, I was thinking of doing interpolation, linear at first, to keep it simple.

You should not need any interpolation. I don't think I explained very clearly what I meant when I asked about the samples per period. The Hilbert filter transfer functions I posted previously use normalized wave numbers on the x-axis. The normalized wave number 1 corresponds to the Nyquist frequency, i.e., two samples per period. By design, these filters are symmetric around the normalized wave number 0.5, and have an upper and lower cut-off where the transfer function rapidly goes to zero. The upper cut-off is somewhere just over two samples per period, so should not really be a concern since you would normally over sample. The lower cut-off, however, might be more problematic. In the example below, it is around a normalized wave number of 0.033, which corresponds to about 60 samples per period. If you sample with, say, 200 samples per period, this filter will give you large errors. The problem is that you cannot easily reduce the lower cut-off frequency without using a large number of filter coefficients, aside from using pre-filters that amplify the signal at wave numbers 0 and 1. But even using pre-filters, the filter kernel size becomes large for low cut-off frequencies. There are other, more complicated ways to deal with this. The reason I was asking about the maximum number of samples per period you intend to use was to figure out whether you would need the more complicated algorithms or not. If you can keep the maximum number of samples per period below 80, you can probably use just a simple filter. I'll try to implement it with pyramid decomposition once I get over my jet lag. It might take me a little while, it's been over twenty years since I last worked on this type of stuff.

[RoGeorge]

That's for sure I didn't understand what yo meant. 

Half of the terms you used to describe what you did there are brand new to me.  Didn't lookup those yet because I'm a beginner in DSP, and the things you mentioned seemed to me like very specialized terms, and that might take many hours or days of reading before understanding.

Will search them one by one these days, very curious to learn about, just that I've sunk all my time from the last couple of days by digging into 108 forks of python-IVI (almost abandonware) then merging and debugging that, and learning how to download the ADC samples from the oscilloscope to the PC.  It works now, just that the download is too slow, almost unusable slow (20 minutes to download 24 million bytes of ADC samples over a 100Mbps LAN).

So far it seems the bottleneck is the oscilloscope itself, and asked about the download speed in another topic so to keep this one about data post processing only:  https://www.eevblog.com/forum/testgear/download-speed-from-rigol-ds1054z-or-similar-oscilloscope-to-a-pc/

[DiTBho]

download is too slow, almost unusable slow

that's why I suggested a usb-dso.

[gf]

download is too slow, almost unusable slow

that's why I suggested a usb-dso.

With a binary bulk transfer over Gigabit Ethernet it should be basically be possible to transfer 24MBytes in ~0.25 seconds. Or with 100Mbit/s Ethernet, in ~2.5 seconds. Obviously not the medium is the bottleneck here.

[2N3055]

download is too slow, almost unusable slow

that's why I suggested a usb-dso.

With a binary bulk transfer over Gigabit Ethernet it should be basically be possible to transfer 24MBytes in ~0.25 seconds. Or with 100Mbit/s Ethernet, in ~2.5 seconds. Obviously not the medium is the bottleneck here.

Nobody said Ethernet was the bottleneck. Pretty much all embedded scopes (those not made on x86 platform and windows) have slow ethernet speeds because of how it is implemented).
USB scopes (like Picoscope) are generally designed to move data as fast as possible to PC because that is how they work in a first place..  So in case of newer USB3 picos, they shovel data very fast to PC:

[DiTBho]

USB scopes (like Picoscope) are generally designed to move data as fast as possible to PC

Yup. Even cheap ones (40Mhz bandwidth) are able to move up to 20Mbyte/sec.

[DiTBho]

There are also DSOs and MSOs able to export the amplitude (peak-peak) and phase measures.
So, this could be a working environment:

Code: [Select]

.          ___________        ________        ________
          |           |      |        |      |       |
freq ---> |           | ---> | system |---> ch0      | ---> { amplitude, phase }
          |    DDS    |      |________|      |       |
          | sin(freq) |                      |  DSO  |
          |           |                      |       |
          |           | ------------------> ch1      | ---> { amplitude, phase }
          |___________|                      |_______|

- SPI, change DDS freq
- loop: RS232/usb/eth, read measures
- goto loop, until stable results
- eval: { gain = (amplidute.ch0 / amplitude.ch1), diff_phase = (phase.ch0 - phase.ch1) }

[nctnico]

It works now, just that the download is too slow, almost unusable slow (20 minutes to download 24 million bytes of ADC samples over a 100Mbps LAN).

That is because you are doing this wrong. Either the libraries you are using are faulty or you are requesting data in a format that needs to be pro-processed (like CSV). The fastest way for the oscilloscope to serve you data is in binary format so use that. Even then there are large differences between DSOs. From testing a couple of years ago, GW Instek is among the fastest DSOs to get data from.

[DiTBho]

GW Instek is among the fastest DSOs to get data from.

which model is? usb? ethernet? is the protocol documented? 
(second hand purchase possible, for me, just "because science")

[RoGeorge]

Dealt with the speed issue, the download went from 23 minutes to only 30 seconds for the entire 24MSa of ADC data.

There are many aspects about the download, so to not hijack this topic from its DSP theme, I've opened another topic about the download performance only, please let's keep this one about the DSP processing:
https://www.eevblog.com/forum/testgear/download-speed-from-rigol-ds1054z-or-similar-oscilloscope-to-a-pc/

[hubi]

I have been trying to get the Hilbert transform to work on a pyramid, but the error is very high (1-2%). Unfortunately, I cannot find my old code, so it might take me a while to figure this out. Here is what I'm getting for a geometric chirp, the recovered amplitude should be one everywhere. The gray line shows the input wave form. I'll try with some better interpolation filters for pyramid merging and see if that helps.