Measure Total Harmonic Distortion (THD) on Rigol DS1000Z Series Oscilloscope

Rigol FFT THD

In MrCarlsonsLab's electronics tutorials, the detailed design and build of a simple oscillator circuit is shown. To assess the quality of the oscillator, the Total harmonic distortion is measured using a Spectrum analyser... I don't have one... But on my bench there is a Rigol DS1054Z (opened up to 100MHz). It can do FFT. It's rubbish, but it can do it. However measuring THD based on the FFT is to my knowledge not possible on the oscilloscope directly.

However, if the data is transferred and analysed, distortion (and frequency) can be calculated with quite reasonable accuracy.

Here a python script that does exactly that.

Rigol DS1000Z - Total Harmonic Distortion Analysis

Requirements:
# python 3.x.x
# pyvisa (http://pyvisa.sourceforge.net/) pip install pyvisa
# numpy (http://numpy.scipy.org/) pip install numpy
# MatPlotLib (http://matplotlib.sourceforge.net/) pip install matplotlib
# peakutils (https://pypi.python.org/pypi/PeakUtils) pip install peakutils
# scipy (https://www.scipy.org/) pip install scipy
# National Instruments VISA Driver: http://www.ni.com/visa/ (install with default options)

Preparation

  1. Install the above packages and driver
    a) packages: pip install pyvisa numpy matplotlib peakutils scipy
    b) download and install: http://www.ni.com/visa/
  2. Connect your Rigol DS1000Z series scope via to the lan network using a lan cable.
  3. On the oscilloscope, go to UTILITY-->IO Setting-->LAN Conf.
  4. Read the "IP address". It will probably be something like 192.168.XXXX.XXX)
  5. Enter it below in the variable scopeIP.

 MEASURE

  1. Enter the remaining variables, such as CHANNEL.
  2. You can choose to do averaging and oversampling to reduce random noise.
  3. Make sure to enter the reasonable plot limits for the X axis. The default values will show a window between 0 and 12 KHz
  4. To measure FFT, set up the oscilloscope (no need to enable FFT), so that you can see the signal which you want to do an FFT on. The more data points in the transient the better. Once you see it, either leave the oscilloscope running or stop it. Your choice. Won't make a difference.
  5. Run this script. (The oscilloscope is stopped automatically before data transfer).
  6. After that the oscilloscope will be set back into RUN mode.

Data is transferred in chunks pretty much as described in the Rigol Programming Manual.
(Hence the data acquisition section below is a bit bloated. Can't do anything about it.)

Next the data is averaged and oversampled to reduce random noise.

The data is then Fourier transformed in two different ways and converted to dBV.

The frequency of the transient is calculated based on the oversampled data.

Total harmonic distortion (THD+N) is calculated

The Voltage data is plotted for reference.

The THD and frequency values are displayed in the console and on the plots, i.e. this script will also work, if you comment out the plot section.

# -*- coding: utf-8 -*- 
#!/usr/bin/env python3

"""

THD and FFT analysis v1.0

Release notes
v0.9: beta release
v1.0: Important Bugfix in the data transfer section! There were indents missing after pasting this code on the website. The error manifests itself in unusually noisy FFT (rectangle) data.

Created on Sun Aug 6 18:02:16 2017 @author: Peter Lachmann THD and frequency analysis functions are adapted from https://gist.github.com/endolith The software is licensed under The MIT License (MIT) Copyright (c) 2017 Peter Lachmann, existingwheel@gmail.com Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE Requirements: # python 3.x.x # pyvisa (http://pyvisa.sourceforge.net/) pip install pyvisa # numpy (http://numpy.scipy.org/) pip install numpy # MatPlotLib (http://matplotlib.sourceforge.net/) pip install matplotlib # peakutils (https://pypi.python.org/pypi/PeakUtils) pip install peakutils # scipy (https://www.scipy.org/) pip install scipy # National Instruments VISA Driver: http://www.ni.com/visa/ (install with default options) How to: ##### INSTALL AND SETUP ##### 1. Install the above packages and driver a) packages: pip install pyvisa numpy matplotlib peakutils scipy b) download and install: http://www.ni.com/visa/ 2. Connect your Rigol DS1000Z series scope via to the lan network using a lan cable. 3. On the oscilloscope, go to UTILITY-->IO Setting-->LAN Conf. 4. Read the "IP address". It will probably be something like 192.168.XXXX.XXX) 5. Enter it below in the variable scopeIP. ##### MEASURE ##### 1. Enter the remaining variables, such as CHANNEL. 2. You can choose to do averaging and oversampling to reduce random noise. 3. Make sure to enter the reasonable plot limits for the X axis. The default values will show a window between 0 and 12 KHz. 4. To record the data, set up the oscilloscope (no need to enable FFT), so that you can see the signal which you want to do an FFT on. The more data points in the transient the better. Once you see it, either leave the oscilloscope running or stop it. Your choice. Won't make a difference. 5. Run this script. (The oscilloscope is stopped automatically before data transfer). 6. After that the oscilloscope will be set back into RUN mode. 7. Data is transferred in chunks pretty much as described in the Rigol Programming Manual. (Hence the data acquisition section below is a bit bloated. Can't do anything about it.) Next the data is averaged and oversampled to reduce random noise. The data is then Fourier transformed in two different ways and converted to dBV. The frequency of the transient is calculated based on the oversampled data. Total harmonic distortion (THD+N) is calculated The Voltage data is plotted for reference. The THD and frequency values are displayed in the console and on the plots, i.e. this script will also work, if you comment out the plot section. """ import visa import numpy as np import matplotlib.pyplot as plt import scipy.fftpack from scipy.signal import blackmanharris, fftconvolve from scipy.stats import binned_statistic from scipy.optimize import leastsq import peakutils from peakutils.plot import plot as pplot import time from matplotlib.mlab import find #################################################### ### Set these variables prior to acquiring data ### #################################################### #plot limits (frequency in Hz) minPlotLimit = 0 #Hz maxPlotLimit = 1.2e4 #Hz #oversampling? If not, set to 1 window = 32 #bin size. #averages #number of readouts from the oscilloscope averaged on the PC. #Not the number of averages set on the scope. average = 1 #Oscilloscope channel channel = "CHANNEL1" #Oscilloscope IP scopeIp = "192.168.XXX.XXX" #see "how to" info above. #this uses port 618 and 111. So, if this has to go through a firewall (router), #then these ports should be opened. #Generate resource address resourceAddress = "TCPIP0::"+scopeIp+"::INSTR" #Rigol ## def sampleRate(transient): #return samplerate as float. This assumes data on a linear x axis. return round(float(1/(transient[0][1]- transient[0][0])),5) def calculateFFT(transient, minPlotLimit, maxPlotLimit): #here we do the DFT (Rectangle) and convert to dbV #xinc is the time between two points. We assume that we have a linear time axis xinc = float(transient[0][1] - transient[0][0]) sample_rate = 1/xinc ############fourier transform y axis############### f = np.fft.rfft(transient[1]) #blackmann-harris dft i = np.argmax(abs(f)) #convert to db yf = 20 * np.log10(f) #generate x axis from time per point and number of points. also window average xf = np.linspace(0.0, int((sample_rate/2)+0.5), len(f)) #find fundamental frequency fundamental = sample_rate * (i / len(transient[1])) #print('Frequency: %f Hz' % fundamental) # Not exac print('Approximate Frequency: %f Hz' % fundamental) # Not exac #why is the spectrum twice as long as expected? print("length xf: " + str(len(xf))) print("length yf: " + str(len(yf))) return [xf[1:], yf[1:]] def dataTransfer(channel): time.sleep(5) #wait until acquisition has finished... I should really predict the wait here. But this works ok for transient <=2s scope.write(":STOP") #stop the scope to transfer data horizontalOffset = round(float(scope.query(":TIMebase:MAIN:OFFSet?")),15) #measure the timebase offset scope.write(":WAVeform:SOURce " + str(channel)) #set waveform source (channel) scope.write(":WAV:MODE RAW") #define data transfer mode scope.write(":WAV:FORM BYTE") #define data transfer format pre = scope.query(":waveform:preamble?").split(',') #read info that allows us to resote, x and y axis units. wFormat = int(pre[0]) #waveform format wType = int(pre[1]) #waveform type wPoints = int(pre[2]) #number of points wCount = int(pre[3]) #number of averages xincrement = float(pre[4]) #time between two neighboring points xorigin = float(pre[5]) #start time of the waveform data xreference = int(float(pre[6])) #reference time of data yincrement = float(pre[7]) #voltage between two points yorigin = float(pre[8]) #vertical offset relative to the reference position yreference = int(float(pre[9])) #vertical reference posisiton #unbelievable, the scope only allows us to transfer 250000 points at a time. #so we have to faff about a bit. pointsLeft = wPoints askTimes = round((wPoints/250000)+0.5) moduloPoints = wPoints % 250000 n = 0 raw_data=0 firstRun = 1 print('Number of Points: %s' % wPoints) print('Number of Queries: %s' % askTimes) print('Downloading data.') while askTimes: n = n + 1 scope.write(":WAV:STAR %s" % n) #start the data transfer if pointsLeft > moduloPoints: print('.',end='') n=n+249999 else: print('-',end='') n=n+moduloPoints-1 ###MOPET XXX -1? are you sure? scope.write(":WAV:STOP %s" % n) pointsLeft = wPoints - n # Actually read waveform data block if firstRun == 1: raw_dataList = scope.query_binary_values('waveform:data?', datatype='B', is_big_endian=True) raw_data = np.array(raw_dataList) firstRun = 0 else: raw_dataList = scope.query_binary_values('waveform:data?', datatype='B', is_big_endian=True) raw_data = np.append(raw_data, np.array(raw_dataList)) askTimes = askTimes - 1 print("") #carriage return #convert raw data to voltage Volts = ((raw_data - yorigin - yreference) * yincrement) #assemble time axis Time = np.arange(xorigin-horizontalOffset, ((xincrement * len(Volts)) + xorigin) - horizontalOffset , xincrement) scope.write(":RUN") #start the scope again return [Time, Volts] def calculateXinc2D(data): return data[0][1]-data[0][0] def rms_flat(a): """ Nicked from Endolith https://gist.github.com/endolith/246092 Return the root mean square of all the elements of *a*, flattened out. """ return np.sqrt(np.mean(np.absolute(a)**2)) def find_range(f, x): """ Nicked from Endolith https://gist.github.com/endolith/246092 Find range between nearest local minima from peak at index x """ for i in np.arange(x+1, len(f)): if f[i+1] >= f[i]: uppermin = i break for i in np.arange(x-1, 0, -1): if f[i] <= f[i-1]: lowermin = i + 1 break return (lowermin, uppermin) def THDN(signal, sample_rate, minPlotLimit, maxPlotLimit): """ Adapted from Endolith https://gist.github.com/endolith/246092 Measure the THD+N for a signal and print the results Prints the estimated fundamental frequency and the measured THD+N. This is calculated from the ratio of the entire signal before and after notch-filtering. Currently this tries to find the "skirt" around the fundamental and notch out the entire thing. """ signal[1] -= np.mean(signal[1]) # Get rid of DC. The signal will now be centered around 0V. windowed = signal[1] * blackmanharris(len(signal[1])) # window the signal using a blackman-harris window total_rms = rms_flat(windowed) # Measure the total signal before filtering but after windowing ############fourier transform y axis############### f = np.fft.rfft(windowed) #blackmann-harris dft i = np.argmax(abs(f)) # Find the highest peak of the spectrum #convert to db yf = 20 * np.log10(f) #generate x axis from time per point and number of points. also window average xf = np.linspace(0.0, int(sample_rate/2+0.5), len(f)) fundamental = sample_rate * (i / len(windowed)) #find fundamental frequency print('Frequency: %f Hz' % fundamental) # Not exact try: lowermin, uppermin = find_range(abs(f), i) f[lowermin: uppermin] = 0 except: print("###########FAILFAILFAILFAILFAILFAILFAILFAILFAILFAILFAILFAILFAILFAILFAIL#############") print("This will fail. I blame the oscilloscope. (Although I could try to catch this earlier and do some sort of magic in the data transfer section") print("http://www.eevblog.com/forum/testgear/rigol-ds1000z-series-buglist-continued-(from-fw-00-04-04-03-02)/msg1263237/#msg1263237") print("Simply run again.") print("Chances are, it will work.") print("####################################################################################") # Transform noise back into the signal domain and measure its rms noise = np.fft.irfft(f) THDN = rms_flat(noise) / total_rms thdnPerCent = THDN * 100 #convert to % Total harmonic distortion thdnDb = 20 * np.log10(THDN) #convert to dB print("THD+N: %.4f%% or %.1f dB" % (THDN * 100, 20 * np.log10(THDN))) print("Total RMS: "+ str(total_rms)) print("RMS flat: " + str(rms_flat(noise))) return [xf, yf, fundamental, thdnPerCent, thdnDb] def freqFromAutocorr(sig): """Adapted from Endolith https://gist.github.com/endolith/246092 Estimate frequency using autocorrelation """ sampleRate = 1/calculateXinc2D(sig) # Calculate autocorrelation (same thing as convolution, but with # one input reversed in time), and throw away the negative lags corr = fftconvolve(sig[1], sig[1][::-1], mode='full') corr = corr[len(corr)//2:] # Find the first low point d = np.diff(corr) start = find(d > 0)[0] # Find the next peak after the low point (other than 0 lag). This bit is # not reliable for long signals, due to the desired peak occurring between # samples, and other peaks appearing higher. # Should use a weighting function to de-emphasize the peaks at longer lags. peak = np.argmax(corr[start:]) + start px, py = parabolic(corr, peak) return sampleRate / px def parabolic(f, x): """ Nicked from Endolith https://gist.github.com/endolith/246092 Quadratic interpolation for estimating the true position of an inter-sample maximum when nearby samples are known. f is a vector and x is an index for that vector. Returns (vx, vy), the coordinates of the vertex of a parabola that goes through point x and its two neighbors. Example: Defining a vector f with a local maximum at index 3 (= 6), find local maximum if points 2, 3, and 4 actually defined a parabola. In [3]: f = [2, 3, 1, 6, 4, 2, 3, 1] In [4]: parabolic(f, argmax(f)) Out[4]: (3.2142857142857144, 6.1607142857142856) """ xv = 1/2. * (f[x-1] - f[x+1]) / (f[x-1] - 2 * f[x] + f[x+1]) + x yv = f[x] - 1/4. * (f[x-1] - f[x+1]) * (xv - x) return (xv, yv) def oversample(transient, window): #binning, averaging and decimation to clean up the signal x = binned_statistic(transient[0][:len(transient[0])], transient[1][:len(transient[0])], bins=int(len(transient[0])/window))[1] y = binned_statistic(transient[0][:len(transient[0])], transient[1][:len(transient[0])], bins=int(len(transient[0])/window))[0] return [x[:int(len(transient[0])/window)],y[:int(len(transient[0])/window)]] def zoom(transient, fundamental, nPeriods): #find the first n periods of the signal minLim = transient[0][0] maxLim = (1/fundamental)*nPeriods return [minLim, maxLim] def prepareTimePlot(title=""): #plot plt.ylabel("Volts") plt.xlabel("Time (s)") plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0)) plt.title(title) def prepareFftPlot(title=""): plt.ylabel("dBV") plt.xlabel("Hz") plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0)) plt.title(title) ############ ### MAIN ### ############ #Initialise Oscilloscope Connection resourceManager = visa.ResourceManager() scope = resourceManager.open_resource(resourceAddress) #Verify communications idn = scope.query("*IDN?") print("Connected to " + idn) #transfer data and average and oversample averageTransient = dataTransfer(channel) transient = averageTransient for n in range(1,average): transient = dataTransfer(channel) averageTransient[1] = transient[1] + averageTransient[1] transient[1] = averageTransient[1] / average #now that the averaging is done, do the oversampling. transientOvs = oversample(transient, window) #calculate Discrete FT (rectangle), convert to dbV and plot fft = calculateFFT(transientOvs, minPlotLimit, maxPlotLimit) #Calculate total harmonic distortion using discrete FT (blackmann-harris window) rfft = THDN(transientOvs, sampleRate(transientOvs), minPlotLimit, maxPlotLimit) #Find frequency: frequency = freqFromAutocorr(transientOvs) #############PLOT################ try: plt.close() except: print("") ax1 = plt.subplot(2,2,1) #plt.xlim(minPlotLimit, maxPlotLimit) prepareFftPlot("DFT Rectangle") plt.plot(fft[0], fft[1]) props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) text = "Fundamental: " + str(frequency) + "Hz\nTHD+N: " + str(rfft[3]) + "%" ax1.text(0.05, 0.1, text, transform=ax1.transAxes, fontsize=7, verticalalignment='top', bbox=props) #set plot limits ax2 = plt.subplot(2,2,3, sharex=ax1) plt.xlim(minPlotLimit, maxPlotLimit) prepareFftPlot("DFT Blackman-Harris") plt.plot(rfft[0], rfft[1]) #Plot Original Data minLim, maxLim = zoom(transient, rfft[2], 2) #find plot limits (zoom in) ax3 = plt.subplot(2,2,2) #plt.xlim(minLim, maxLim) prepareTimePlot("Original Data") plt.plot(transient[0], transient[1]) #Plot Oversampled Data minLim, maxLim = zoom(transientOvs, rfft[2], 2) #find plot limits (zoom in) ax4 = plt.subplot(2,2,4, sharex=ax3) #plt.xlim(minLim, maxLim) prepareTimePlot("Oversampled Data") plt.plot(transientOvs[0], transientOvs[1]) plt.tight_layout() plt.show()

Posted in Code, Misc Electronics and tagged , , , , , .

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.