Fork PyRPL on GitHub

# pyrpl.hardware_modules.iir package¶

## pyrpl.hardware_modules.iir.iir module¶

class pyrpl.hardware_modules.iir.iir.IIR(rp, name)[source]

Setup Attributes:

• loops: Decimation factor of IIR w.r.t. 125 MHz. Must be at least 5.
• zeros:
• poles:
• gain:
• on: IIR is on
• bypass: IIR is bypassed
• data_curve: NA curve id to use as a basis for the graphical filter design
bypass

IIR is bypassed

coefficients
complex_poles
complex_zeros
data_curve

NA curve id to use as a basis for the graphical filter design

end = 'zeros'
gain
iirfilter = None
internal_overflow

returns True if the IIR filter has experienced an internal overflow (leading to saturation) since the last reset

loops

Decimation factor of IIR w.r.t. 125 MHz. Must be at least 5.

on

IIR is on

output_saturation

returns True if the output of the IIR filter has saturated since the last reset

overflow

a string indicating the overflow status of the iir module

overflow_bitfield

Bitmask for various overflow conditions

poles
real_poles
real_zeros
sampling_time
select_pole_or_zero(value, logdist=True, search_in=['real_poles', 'real_zeros', 'complex_poles', 'complex_zeros'])[source]

selects the pole or zero closest to value

logdist=True computes the distance in logarithmic units search_in may be used to restrict the search to certain sublists

setup(**kwds)

Setup an IIR filter

the transfer function of the filter will be (k ensures DC-gain = g):

(s-2*pi*z[0])*(s-2*pi*z[1])...
H(s) = k*——————-
(s-2*pi*p[0])*(s-2*pi*p[1])...
coefficients data to be passed to iir.bodeplot to plot the
realized transfer function
input

selects the input signal of the module

loops

Decimation factor of IIR w.r.t. 125 MHz. Must be at least 5.

zeros
poles
output_direct

selects to which analog output the module signal is sent directly

inputfilter

Input filter bandwidths [Hz]. 0 = off, positive bandwidth <=> lowpass, negative bandwidth <=> highpass.

gain
on

IIR is on

bypass

IIR is bypassed

data_curve

NA curve id to use as a basis for the graphical filter design

start = 'complex'
transfer_function(frequencies, extradelay=0, kind='final')[source]

Returns a complex np.array containing the transfer function of the current IIR module setting for the given frequency array. The best-possible estimation of delays is automatically performed for all kinds of transfer function. The setting of ‘bypass’ is ignored for this computation, i.e. the theoretical and measured transfer functions can only agree if bypass is False.

Parameters: frequencies (np.array or float) – Frequencies to compute the transfer function for extradelay (float) – External delay to add to the transfer function (in s). If zero, only the delay for internal propagation from input to output_signal is used. If the module is fed to analog inputs and outputs, an extra delay of the order of 150 ns must be passed as an argument for the correct delay modelisation. kind (str) – The IIR filter design is composed of a number of steps. Each step slightly modifies the transfer function to adapt it to the implementation of the IIR. The various intermediate transfer functions can be helpful to debug the iir filter. kind should be one of the following (default is ‘implemented’): - ‘all’: returns a list of data to be passed to iir.bodeplot with all important kinds of transfer functions for debugging ‘continuous’: the designed transfer function in continuous time ‘before_partialfraction_continuous’: continuous filter just before partial fraction expansion of the coefficients. The partial fraction expansion introduces a large numerical error for higher order filters, so this is a good place to check whether this is a problem for a particular filter design ‘before_partialfraction_discrete’: discretized filter just before partial fraction expansion of the coefficients. The partial fraction expansion introduces a large numerical error for higher order filters, so this is a good place to check whether this is a problem for a particular filter design ‘before_partialfraction_discrete_zoh’: same as previous, but zero order hold assumption is used to transform from continuous to discrete ‘discrete’: the transfer function after transformation to discrete time ‘discrete_samplehold’: same as discrete, but zero delay between subsequent biquads is assumed ‘highprecision’: hypothetical transfer function assuming that 64 bit fixed point numbers were used in the fpga (decimal point at bit 48) ‘implemented’: transfer function after rounding the coefficients to the precision of the fpga tf (np.array(..., dtype=np.complex)) – The complex open loop transfer function of the module. If kind==’all’, a list of plotdata tuples is returned that can be passed directly to iir.bodeplot().
zeros
class pyrpl.hardware_modules.iir.iir.IirComplexListProperty(*args, **kwargs)[source]

slave property to store complex part of zeros and poles

validate_and_normalize_element(obj, val)[source]

real part should be strictly negative. imaginary part is in principle arbitrary, but will be kept positive for simplicity.

class pyrpl.hardware_modules.iir.iir.IirFloatListProperty(*args, **kwargs)[source]

slave property to store real part of zeros and poles

list_changed(module, operation, index, value=None)[source]

make sure that an element from one of the four lists is selected at once

validate_and_normalize_element(obj, val)[source]

makes sure that real poles are strictly positive. val=0 is turned into val=-1.

value_updated(obj, value=None, appendix=[])[source]
class pyrpl.hardware_modules.iir.iir.IirListProperty(min=<MagicMock name='mock.inf.__neg__()' id='139827380281936'>, max=<MagicMock name='mock.inf' id='139827380218960'>, increment=0, log_increment=False, **kwargs)[source]

master property to store zeros and poles

default = []
get_value(obj)[source]

the master’s getter collects its value from the real and complex list

set_value(obj, value)[source]

the master’s setter writes its value to the slave lists

validate_and_normalize(obj, value)[source]

Converts the value in a list of float numbers.

validate_and_normalize_element(obj, val)[source]
class pyrpl.hardware_modules.iir.iir.OverflowProperty(default=None, doc='', ignore_errors=False, call_setup=False)[source]
get_value(obj)[source]
set_value(obj, value)[source]
validate_and_normalize(obj, value)[source]
class pyrpl.hardware_modules.iir.iir.SignalLauncherIir(module)[source]
update_plot = <MagicMock name='mock.QtCore.Signal()' id='139827380610576'>

## pyrpl.hardware_modules.iir.iir_theory module¶

class pyrpl.hardware_modules.iir.iir_theory.IirFilter(zeros, poles, gain, loops=None, dt=8e-09, minloops=4, maxloops=1023, iirstages=16, totalbits=32, shiftbits=29, tol=0.001, frequencies=None, inputfilter=0, moduledelay=5)[source]

Bases: object

Computes coefficients and predicts transfer functions of an IIR filter.

Parameters: sys ((zeros, poles, gain)) – zeros: list of complex zeros poles: list of complex poles gain: DC-gain zeros/poles with nonzero imaginary part should come in complex conjugate pairs, otherwise the conjugate zero/pole will automatically be added. After this, the number of poles should exceed the number of zeros at least by one, otherwise a real pole near the Nyquist frequency will automatically be added until there are more poles than zeros. loops (int or None) – the number of FPGA cycles per filter sample. None tries to automatically find the value leading to the highest possible sampling frequency. If the numerical precision of the filter coefficients in the FPGA is the limiting, manually setting a higher value of loops may improve the filter performance. dt (float) – the FPGA clock frequency. Should be very close to 8e-9 minoops (int) – minimum number of loops (constant of the FPGA design) maxloops (int) – maximum number of loops (constant of the FPGA design) tol (float) – tolerancee for matching conjugate pole/zero pairs. 1e-3 is okay. intermediatereturn (str or None) – if set to a valid option, the algorithm will stop at the specified step and return an intermediate result for debugging. Valid options are coefficients coefficients is an array of float arrays of length six, which hold the filter coefficients to be passed directly to the iir module IirFilter.loops of the IirFilter instance is automatically corrected to a number of loops compatible with the implemented design.
coefficients

Computes and returns the coefficients of the IIR filter set by IirFilter.sys.

Parameters: sys ((zeros, poles, gain)) – zeros: list of complex zeros poles: list of complex poles gain: DC-gain zeros/poles with nonzero imaginary part should come in complex conjugate pairs, otherwise the conjugate zero/pole will automatically be added. After this, the number of poles should exceed the number of zeros at least by one, otherwise a real pole near the nyquist frequency will automatically be added until there are more poles than zeros. loops (int or None) – the number of FPGA cycles per filter sample. None tries to automatically find the value leading to the highest possible sampling frequency. If the numerical precision of the filter coefficients in the FPGA is the limiting, manually setting a higher value of loops may improve the filter performance. dt (float) – the FPGA clock frequency. Should be very close to 8e-9 minoops (int) – minimum number of loops (constant of the FPGA design) maxloops (int) – maximum number of loops (constant of the FPGA design) tol (float) – tolerancee for matching conjugate pole/zero pairs. 1e-3 is okay. intermediatereturn (str or None) – if set to a valid option, the algorithm will stop at the specified step and return an intermediate result for debugging. Valid options are coefficients coefficients is an array of float arrays of length six, which hold the filter coefficients to be passed directly to the iir module IirFilter.loops of the IirFilter instance is automatically corrected to a number of loops compatible with the implemented design.
coefficients_rounded
designdata
finiteprecision(coeff=None, totalbits=None, shiftbits=None)[source]
frequencies
minimize_delay(coefficients=None)[source]

Minimizes the delay of coefficients by rearranging the biquads in an optimal way (highest frequency poles get minimum delay.

Parameters: coefficients – new coefficients
prewarp(z, p, dt=8e-09)[source]

prewarps frequencies in order to correct warping effect in discrete time conversion

proper_sys

Makes sure that a system is strictly proper and that all complex poles/zeros have conjugate parters.

Parameters: zeros (list of zeros) – poles (list of poles) – loops (number of loops to implement. Can be None for autodetection.) – minloops (minimum number of loops that is acceptable) – maxloops (minimum number of loops that is acceptable) – iirstages (number of biquads available for implementation) – tol (tolerance for matching complex conjugate partners) – (zeros, poles, minloops) - the corrected lists of poles/zeros and the number of loops that are minimally needed for implementation
rescaled_sys

rescales poles and zeros with 2pi and returns the prefactor corresponding to dc-gain gain

rp2coefficients(r, p, c, tol=0)[source]

Pairs residues and corresponding poles into second order sections.

Parameters: r (array with numerator coefficient) – p (array with poles) – tol (tolerance for combining complex conjugate pairs.) – coefficients array((N, 6), dtype=np.float64) where N is number of biquads
sampling_rate
sys
tf_coefficients(frequencies=None, coefficients=None, delay=False)[source]

computes implemented transfer function - assuming no delay and infinite precision (actually floating-point precision) Returns the discrete transfer function realized by coefficients at frequencies.

Parameters: coefficients (np.array) – coefficients as returned from iir module frequencies (np.array) – frequencies to compute the transfer function for dt (float) – discrete sampling time (seconds) zoh (bool) – If true, zero-order hold implementation is assumed. Otherwise, the delay is expected to depend on the index of biquad. np.array(.., dtype=np.complex)
tf_continuous(frequencies=None)[source]

Returns the continuous transfer function of sys at frequencies.

Parameters: sys (tuple) – (zeros, poles, gain) zeros: list of complex zeros poles: list of complex poles gain: float frequencies (np.array) – frequencies to compute the transfer function for np.array(.., dtype=np.complex)
tf_discrete(rp_discrete=None, frequencies=None)[source]

Returns the discrete transfer function realized by coefficients at frequencies.

Parameters: rpz (np.array) – coefficients as returned from iir module (array of biquad coefficients) frequencies (np.array) – frequencies to compute the transfer function for dt (float) – discrete sampling time (seconds) delay_per_cycle (float) – the biquad at coefficients[i] experiences an extra delay of i*delay_per_cycle zoh (bool) – If true, zero-order hold implementation is assumed. Otherwise, the delay is expected to depend on the index of biquad. np.array(.., dtype=np.complex)
tf_final(frequencies=None)[source]

Returns the discrete transfer function realized by coefficients at frequencies.

Parameters: coefficients (np.array) – coefficients as returned from iir module frequencies (np.array) – frequencies to compute the transfer function for dt (float) – discrete sampling time (seconds) zoh (bool) – If true, zero-order hold implementation is assumed. Otherwise, the delay is expected to depend on the index of biquad. np.array(.., dtype=np.complex)
tf_inputfilter(inputfilter=None, frequencies=None)[source]
tf_partialfraction(frequencies=None)[source]

Returns the transfer function just before the partial fraction expansion for frequencies.

Parameters: sys ((poles, zeros, k)) – dt (sampling time) – continuous (if True, returns the transfer function in continuous) – time domain, if False converts to discrete one method (method for scipy.signal.cont2discrete) – alpha (alpha for above method (see scipy documentation)) – np.array(.., dtype=np.complex)
tf_rounded(frequencies=None, delay=False)[source]

Returns the discrete transfer function realized by coefficients at frequencies.

Parameters: coefficients (np.array) – coefficients as returned from iir module frequencies (np.array) – frequencies to compute the transfer function for dt (float) – discrete sampling time (seconds) zoh (bool) – If true, zero-order hold implementation is assumed. Otherwise, the delay is expected to depend on the index of biquad. np.array(.., dtype=np.complex)
pyrpl.hardware_modules.iir.iir_theory.bodeplot(data, xlog=False)[source]

plots a bode plot of the data x, y

Parameters: data (a list of tuples (f, tf[, label]) where f are frequencies and tf) – complex transfer data, and label the label for data xlog (sets xaxis to logscale) – figure –
pyrpl.hardware_modules.iir.iir_theory.cont2discrete(r, p, c, dt=8e-09)[source]

Transforms residue and pole from continuous to discrete time

Parameters: r (residues) – p (poles) – dt (sampling time) – (r, p) with the transformation applied
pyrpl.hardware_modules.iir.iir_theory.discrete2cont(r, p, c, dt=8e-09)[source]

Transforms residues and poles from discrete time to continuous

Parameters: r (residues) – p (poles) – c (constant term to be carried along) – dt (sampling time (s)) – r, p with the transformation applied
pyrpl.hardware_modules.iir.iir_theory.freqs(sys, w)[source]

This function computes the frequency response of a zpk system at an array of frequencies.

It loosely mimicks ‘scipy.signal.freqs’.

Parameters: system ((zeros, poles, k)) – zeros and poles both in rad/s, k is the actual coefficient, not DC gain w (np.array) – frequencies in rad/s np.array(.., dtype=np.complex) with the response
pyrpl.hardware_modules.iir.iir_theory.freqs_rp(r, p, c, w)[source]

same as freqs, but takes a list of residues and poles

pyrpl.hardware_modules.iir.iir_theory.freqz_(sys, w, dt=8e-09)[source]

This function computes the frequency response of a zpk system at an array of frequencies.

It loosely mimicks ‘scipy.signal.frequresp’.

Parameters: system ((zeros, poles, k)) – zeros and poles both in rad/s, k is the actual coefficient, not DC gain w (np.array) – frequencies in rad/s dt (sampling time) – np.array(.., dtype=np.complex) with the response
pyrpl.hardware_modules.iir.iir_theory.residues(z, p, k)[source]

this function uses the residue method (Heaviside Cover-up method) to perform the partial fraction expansion of a rational function defined by zeros, poles and a prefactor k. No intermediate conversion into a polynome is performed, which makes this function less prone to finite precision issues. In the current version, no pole value may occur twice and the number of poles must be strictly greater than the number of zeros.

Returns: np.array(dtype=np.complex128) containing the numerator array a of the expansion – product_i( s - z[i] ) a[j] k ————————- = sum_j ( ———- ) – product_j( s - p[j] ) s - p[j]
pyrpl.hardware_modules.iir.iir_theory.sos2zpk(sos)[source]

Return zeros, poles, and gain of a series of second-order sections

Parameters: sos (array_like) – Array of second-order filter coefficients, must have shape (n_sections, 6). See sosfilt for the SOS filter format specification. z (ndarray) – Zeros of the transfer function. p (ndarray) – Poles of the transfer function. k (float) – System gain.

Notes

New in version 0.16.0.

## Module contents¶

Sometimes it is interesting to realize even more complicated filters. This is the case, for example, when a piezo resonance limits the maximum gain of a feedback loop. For these situations, the IIR module can implement filters with ‘Infinite Impulse Response’ (https://en.wikipedia.org/wiki/Infinite_impulse_response). It is the your task to choose the filter to be implemented by specifying the complex values of the poles and zeros of the filter. In the current version of pyrpl, the IIR module can implement IIR filters with the following properties:

• strictly proper transfer function (number of poles > number of zeros)
• poles (zeros) either real or complex-conjugate pairs
• no three or more identical real poles (zeros)
• no two or more identical pairs of complex conjugate poles (zeros)
• pole and zero frequencies should be larger than :math:

rac{f_ m{nyquist}}{1000} (but you can optimize the nyquist frequency of your filter by tuning the ‘loops’ parameter) - the DC-gain of the filter must be 1.0. Despite the FPGA implemention being more flexible, we found this constraint rather practical. If you need different behavior, pass the IIR signal through a PID module and use its input filter and proportional gain. If you still need different behaviour, the file iir.py is a good starting point. - total filter order <= 16 (realizable with 8 parallel biquads) - a remaining bug limits the dynamic range to about 30 dB before internal saturation interferes with filter performance

Filters whose poles have a positive real part are unstable by design. Zeros with positive real part lead to non-minimum phase lag. Nevertheless, the IIR module will let you implement these filters.

In general the IIR module is still fragile in the sense that you should verify the correct implementation of each filter you design. Usually you can trust the simulated transfer function. It is nevertheless a good idea to use the internal network analyzer module to actually measure the IIR transfer function with an amplitude comparable to the signal you expect to go through the filter, as to verify that no saturation of internal filter signals limits its performance.

#reload to make sure settings are default ones
from pyrpl import Pyrpl
p = Pyrpl(hostname="192.168.1.100")

#shortcut
iir = p.rp.iir

#print docstring of the setup function
print(iir.setup.__doc__)

#prepare plot parameters
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10, 6)

#setup a complicated transfer function
zeros = [ 4e4j+300, +2e5j+1000, +2e6j+3000]
poles = [ 1e6, +5e4j+300, 1e5j+3000, 1e6j+30000]
iir.setup(zeros=zeros, poles=poles, loops=None, plot=True)
print("Filter sampling frequency: ", 125./iir.loops,"MHz")


If you try changing a few coefficients, you will see that your design filter is not always properly realized. The bottleneck here is the conversion from the analytical expression (poles and zeros) to the filter coefficients, not the FPGA performance. This conversion is (among other things) limited by floating point precision. We hope to provide a more robust algorithm in future versions. If you can obtain filter coefficients by another, preferrably analytical method, this might lead to better results than our generic algorithm.

Let’s check if the filter is really working as it is supposed:

# first thing to check if the filter is not ok
print("IIR overflows before:", bool(iir.overflow))
na = p.networkanalyzer

# measure tf of iir filter
iir.input = na.iq
na.setup(iq_name='iq1', start=1e4, stop=3e6, points = 301, rbw=100, avg=1,
amplitude=0.1, input='iir', output_direct='off', logscale=True)
tf = na.curve()

# first thing to check if the filter is not ok
print("IIR overflows after:", bool(iir.overflow))

# retrieve designed transfer function
designdata = iir.transfer_function(na.frequencies)

#plot with design data
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10, 6)
from pyrpl.hardware_modules.iir.iir_theory import bodeplot
bodeplot([(na.frequencies, designdata, "designed system"),
(na.frequencies, tf, "measured system")], xlog=True)


As you can see, the filter has trouble to realize large dynamic ranges. With the current standard design software, it takes some ‘practice’ to design transfer functions which are properly implemented by the code. While most zeros are properly realized by the filter, you see that the first two poles suffer from some kind of saturation. We are working on an automatic rescaling of the coefficients to allow for optimum dynamic range. From the overflow register printed above the plot, you can also see that the network analyzer scan caused an internal overflow in the filter. All these are signs that different parameters should be tried.

A straightforward way to impove filter performance is to adjust the DC-gain and compensate it later with the gain of a subsequent PID module. See for yourself what the parameter g=0.1 (instead of the default value g=1.0) does here:

#rescale the filter by 20 fold reduction of DC gain
iir.setup(zeros=zeros, poles=poles, g=0.1, loops=None, plot=False)

# first thing to check if the filter is not ok
print("IIR overflows before:", bool(iir.overflow))

# measure tf of iir filter
iir.input = na.iq
tf = na.curve()

# first thing to check if the filter is not ok
print("IIR overflows after:", bool(iir.overflow))

# retrieve designed transfer function
designdata = iir.transfer_function(na.frequencies)

#plot with design data
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10, 6)
from pyrpl.hardware_modules.iir.iir_theory import bodeplot
bodeplot([(na.frequencies, designdata, "designed system"),
(na.frequencies, tf, "measured system")], xlog=True)


You see that we have improved the second peak (and avoided internal overflows) at the cost of increased noise in other regions. Of course this noise can be reduced by increasing the NA averaging time. But maybe it will be detrimental to your application? After all, IIR filter design is far from trivial, but this tutorial should have given you enough information to get started and maybe to improve the way we have implemented the filter in pyrpl (e.g. by implementing automated filter coefficient scaling).

If you plan to play more with the filter, these are the remaining internal iir registers:

iir = p.rp.iir

# useful diagnostic functions
print("IIR on:", iir.on)
print("IIR bypassed:", iir.shortcut)
print("IIR copydata:", iir.copydata)
print("IIR loops:", iir.loops)
print("IIR overflows:", bin(iir.overflow))
print("Coefficients (6 per biquad):")
print(iir.coefficients)

# set the unity transfer function to the filter
iir._setup_unity()