# Measure Total Harmonic Distortion (THD) on Rigol DS1000Z Series Oscilloscope 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.0Release notesv0.9: beta releasev1.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- transient)),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 - transient)
sample_rate = 1/xinc

############fourier transform y axis###############
f = np.fft.rfft(transient) #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))

#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) #waveform format
wType = int(pre) #waveform type
wPoints = int(pre) #number of points
wCount = int(pre) #number of averages
xincrement = float(pre) #time between two neighboring points
xorigin = float(pre) #start time of the waveform data
xreference = int(float(pre)) #reference time of data
yincrement = float(pre) #voltage between two points
yorigin = float(pre) #vertical offset relative to the reference position
yreference = int(float(pre)) #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-data

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 -= np.mean(signal) # Get rid of DC. The signal will now be centered around 0V.
windowed = signal * blackmanharris(len(signal)) # 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, sig[::-1], mode='full')
corr = corr[len(corr)//2:]

# Find the first low point
d = np.diff(corr)
start = find(d > 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 : f = [2, 3, 1, 6, 4, 2, 3, 1]

In : parabolic(f, argmax(f))
Out: (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[:len(transient)], transient[:len(transient)], bins=int(len(transient)/window))
y = binned_statistic(transient[:len(transient)], transient[:len(transient)], bins=int(len(transient)/window))
return [x[:int(len(transient)/window)],y[:int(len(transient)/window)]]

def zoom(transient, fundamental, nPeriods):
#find the first n periods of the signal
minLim = transient
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 = transient + averageTransient

transient = averageTransient / 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, fft)
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
text = "Fundamental: " + str(frequency) + "Hz\nTHD+N: " + str(rfft) + "%"
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, rfft)

#Plot Original Data
minLim, maxLim = zoom(transient, rfft, 2) #find plot limits (zoom in)
ax3 = plt.subplot(2,2,2)
#plt.xlim(minLim, maxLim)
prepareTimePlot("Original Data")
plt.plot(transient, transient)

#Plot Oversampled Data
minLim, maxLim = zoom(transientOvs, rfft, 2) #find plot limits (zoom in)
ax4 = plt.subplot(2,2,4, sharex=ax3)
#plt.xlim(minLim, maxLim)
prepareTimePlot("Oversampled Data")
plt.plot(transientOvs, transientOvs)

plt.tight_layout()

plt.show()```

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

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