# Mathematics

## C++ equivalent of NumPy/PHP rad2deg

In order to convert radians to degrees, PHP provides the rad2deg function whereas Python’s NumPy library provides np.rad2deg.

In C++ there is no standard function to convert radians to degrees.

However, you can use this simple snippet:

#define _USE_MATH_DEFINES #include <cmath> /** * Convert the angle given in radians to degrees. */ template<typename F> F rad2deg(F angle) { return angle * 180.0 / M_PI; }

This will work for `double`

, `float`

etc. and returns the angle in degrees.

The formula is pretty simple:

\text{Degrees} = \frac{\text{Radians} \cdot 180°}{\pi}## How to use Z0 (characteristic impedance of vacuum) constant in SciPy/NumPy

If you want to use the `Z0`

constant (characteristic impedance of free space) in Python, use this snippet:

import scipy.constants Z0 = scipy.constants.physical_constants['characteristic impedance of vacuum'][0] print(Z0) # 376.73031346177066

In contrast to other constants, `Z0`

is not available directly like `scipy.constants.pi`

but you need to use the `scipy.constants.physical_constants`

dict in order to access it.

## Easily compute & visualize FFTs in Python using UliEngineering

UliEngineering is a mixed data analytics library in Python – one of the utilities it provides is an easy-to-use package to compute FFTs. In contrast to other package, this library is oriented towards practical usecases and allows you to do the FFT in only one line of code! No knowledge of Math required.

### How to install UliEngineering

UliEngineering is a Python 3 only library. Install using `pip`

:

sudo pip3 install -U UliEngineering

### Getting started

First, we’ll generate some test data. See this previous post for more details on how to generate sinusoid test data:

from UliEngineering.SignalProcessing.Simulation import * # Generate test data: 100 Hz + 400 Hz tone data = sine_wave(frequency=100.0, samplerate=1000, amplitude=1.) \ + sine_wave(frequency=400.0, samplerate=1000, amplitude=0.5)

The test data consists of a `100 Hz`

sine wave plus a `400 Hz`

sine wave (with half the amplitude). That signal is sampled with a sample rate of `1000 Hz`

.

Now we can compute & visualize the FFT using matplotlib:

# Compute FFT. Use the same samplerate # NOTE: Window defaults to "blackman"! from UliEngineering.SignalProcessing.FFT import compute_fft fft = compute_fft(data, samplerate=1e3) # Plot from matplotlib import pyplot as plt plt.style.use("ggplot") plt.gcf().set_size_inches(10, 5) # Use (20, 10) to get a larger plot plt.plot(fft.frequencies, fft.amplitudes) plt.xlabel("Frequency") plt.ylabel("Amplitude")

`compute_fft(data, samplerate=1e3)`

returns a tuple `(fftx, ffty)`

. It performs an FFT the size of the input (i.e. since data is a length-`1000`

array, the FFT will be of size `1000`

).

`fft.frequencies`

is an array of frequencies (in Hz), corresponding to the values in `fft.amplitudes`

. You can also use `fft.angles`

to get relative angles in degrees, but that is not being covered in this blogpost.

As you can see in the plot shown above, the maximum frequency that can be detected is always **half the samplerate**, i.e. for our samplerate of f_s = 1000\,\text{Hz} it is 500\,\text{Hz}. See this more math-centric FFT explanation if you want to know more details.

Internally, `compute_fft()`

performs this computation:

- 2 is a correction factor that takes into account that we throw away the latter half of the raw FFT result (since we’re doing a real FFT)
- \frac{1}{\text{len(data)}} normalizes the FFT results so they are indepedent of the data length (i.e. if you pass a longer sample of the same sine wave you will still get the same result
- \text{abs}\left(\cdots\right) Converts the complex phase-aware result of the FFT to a spectrum which is easier to read & visualize.
- Window (which defaults to
`blackman`

) is the window which is applied to the data to alleviate some mathematical effects at the beginning and the end of the dataset. See wikipedia on window functions for more details.*UliEngineering*currently offers this list of window functions:`blackman`

`bartlett`

`hamming`

`hanning`

`kaiser`

(Parameter is fixed to`2.0`

)`none`

### Selecting frequency ranges

Using the UliEngineering API, selecting a frequency range of the FFT is trivially easy: Just use `fft[lowfreq:highfreq]`

. You can use `fft[lowfreq:]`

to select everything starting from `lowfreq`

or use `fft[:highfreq]`

to select everything up to `highfreq`

.

from UliEngineering.SignalProcessing.FFT import compute_fft fft = compute_fft(data, samplerate=1e3) # Select frequency range: 50 to 200 Hz fft = fft[50.0:200.0] # Plot from matplotlib import pyplot as plt plt.style.use("ggplot") plt.gcf().set_size_inches(10, 5) # Use (20, 10) to get a larger plot plt.plot(fft.frequencies, fft.amplitudes) plt.xlabel("Frequency") plt.ylabel("Amplitude") plt.savefig("/ram/fft-frequency-range.svg")

### Extract amplitude & angle at a certain frequency

By using `[frequency]`

, i.e. getitem operator with a single value, you get an FFTPoint() object containing the frequency, amplitude and relative angle for a given frequency. The library automatically selects the closest FFT bucket, so even if your FFT does not have a bucket for that specific frequency, you will get sensible results.

from UliEngineering.SignalProcessing.FFT import compute_fft fft = compute_fft(data, samplerate=1e3) # Show value at a certain frequency print(fft[30]) # FFTPoint(frequency=30.0, value=2.91e-08, angle=0.0) print(fft[100]) # FFTPoint(frequency=100.0, value=0.419, angle=0.0)

### Short FFTs for long data

Using `compute_fft()`

, if we have an extremely long data array, this means we’ll compute an extremely long FFT. In many cases, this is not desirable and you want to compute a fixed-size FFT (most often a power-of-two FFT, e.g. `1024`

, `2048`

, `4096`

etc).

Let’s generate some long test data and assume we want to compute a size-`1024`

FFT on it

from UliEngineering.SignalProcessing.FFT import * fftx, ffty = simple_serial_fft_reduce(data, samplerate=1e3, fftsize=1024)

Plotting the data like above yields

which looks almost exactly like our `compute_fft()`

plot before – just as we expected.

`simple_serial_fft_reduce()`

takes care of all the magic and normalization for us, including partitioning the data into overlapping chunks, adding the FFTs and properly normalizing the result.

The naming convention is significant here:

`simple_...._reduce`

means this is a variant with sensible defaults (`simple`

) for using a`reduction`

function (default:`sum`

) on multiple FFTs.`serial`

means the individual FFTs

### Parallelizing the FFTs

If you have a huge dataset, you can use `simple_parallel_fft_reduce()`

identically to`simple_serial_fft_reduce()`

:

from UliEngineering.SignalProcessing.FFT import * fftx, ffty = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024)

However in most cases you want to initialize the executor manually so you can re-use it later:

from UliEngineering.SignalProcessing.FFT import * from concurrent.futures import ThreadPoolExecutor executor = ThreadPoolExecutor() # No argument => use num_cpus threads fftx, ffty = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024, executor=executor)

We can use a `ThreadPoolExecutor()`

since `scipy.fftpack`

(which `UliEngineering`

uses to do the hard math) unlocks the Python GIL.

Note that due to the need to do a lot of housekeeping tasks, `simple_parallel_fft_reduce()`

is much slower than `simple_serial_fft_reduce()`

if you have a dataset so small that parallelization is not effective. My initial recommendation is to consider using the parallel variant if the total execution time of the serial variant is larger than 0.5\,s

## Easily generate square/triangle/sawtooth/inverse sawtooth waveform data in Python using UliEngineering

In a previous post, I’ve detailed how to generate sine/cosine wave data with frequency, amplitude, offset, phase shift and time offset using only a single line of code.

This post extends this approach by showing how to generate **square wave**, **triangle wave**, **sawtooth wave** and **inverse sawtooth wave** data – still in only one line of code.

All of the parameters, including **frequency, amplitude, offset, phase shift and time offset** apply for those functions as well – see the previous post for details and examples for those parameters.

We are using the UliEngineering library, more specifically the `UliEngineering.SignalProcessing.Simulation`

package:

### How to install UliEngineering

UliEngineering is a Python 3 only library. Install using `pip`

:

sudo pip3 install -U UliEngineering

### Square wave

from UliEngineering.SignalProcessing.Simulation import square_wave data = square_wave(frequency=10.0, samplerate=10e3)

### Triangle wave

from UliEngineering.SignalProcessing.Simulation import triangle_wave data = triangle_wave(frequency=10.0, samplerate=10e3)

### Sawtooth wave

from UliEngineering.SignalProcessing.Simulation import sawtooth data = sawtooth(frequency=10.0, samplerate=10e3)

### Inverse sawtooth wave

from UliEngineering.SignalProcessing.Simulation import inverse_sawtooth data = inverse_sawtooth(frequency=10.0, samplerate=10e3)

### Plotting code

This code was used to generate the plots for this post in Jupyter:

%matplotlib inline from matplotlib import pyplot as plt plt.style.use("ggplot") from UliEngineering.SignalProcessing.Simulation import square_wave data = square_wave(frequency=10.0, samplerate=10e3) # set_size_inches(20, 10) to make it even larger! plt.gcf().set_size_inches(10, 5) plt.plot(data, label="original") plt.savefig("/dev/shm/square-wave.svg")

## Easily generate sine/cosine waveform data in Python using UliEngineering

In order to generate sinusoid test data in Python you can use the UliEngineering library which provides an easy-to-use functions in `UliEngineering.SignalProcessing.Simulation`

:

### How to install UliEngineering

UliEngineering is a Python 3 only library. Install using `pip`

:

sudo pip3 install -U UliEngineering

### Basic example

from UliEngineering.SignalProcessing.Simulation import sine_wave # Default: Generates 1 second of data with amplitude = 1.0 (swing from -1.0 ... 1.0) sine = sine_wave(frequency=10.0, samplerate=10e3) cosine = cosine_wave(frequency=10.0, samplerate=10e3)

### Amplitude & offset

Use `amplitude=0.5`

to specify that the sine wave should swing between `-0.5`

and `0.5`

.

Use `offset=2.0`

to specify that the sine wave should be centered vertically around `2.0`

:

from UliEngineering.SignalProcessing.Simulation import sine_wave data = sine_wave(frequency=10.0, samplerate=10e3, amplitude=0.5, offset=2.0)

### Phaseshift example

You can specify the phaseshift to use – to use a `180°`

phaseshift, just use `phaseshift=180.0`

from UliEngineering.SignalProcessing.Simulation import sine_wave original = sine_wave(frequency=10.0, samplerate=10e3) shifted = sine_wave(frequency=10.0, samplerate=10e3, phaseshift=180.)

### Time delay example

While this is functionally equivalent to phase offset, it is often convenient to specify the *time delay* in seconds instead of the phase shift in degrees:

from UliEngineering.SignalProcessing.Simulation import sine_wave original = sine_wave(frequency=10.0, samplerate=10e3) # Shifted signal is delayed 0.01s = 10 milliseconds compared to original shifted = sine_wave(frequency=10.0, samplerate=10e3, timedelay=0.01)

Note that in the time domain, the signals appear to be shifted backwards when you use a positive timedelay value. This is in accordance with the *delay* naming, implying that the signal is *delayed* by that amount.

You can specify both *phase shift* and *time delay*, meaning that both will be applied (the offset is added)

### Plotting

If you want to debug your signals visually, this is the code that was used inside Jupyter to generate the plots shown above:

%matplotlib inline from matplotlib import pyplot as plt plt.style.use("ggplot") # Generate data from UliEngineering.SignalProcessing.Simulation import sine_wave original = sine_wave(frequency=10.0, samplerate=10e3) shifted = sine_wave(frequency=10.0, samplerate=10e3, timedelay=0.01) # set_size_inches(20, 10) to make it even larger! plt.gcf().set_size_inches(10, 5) plt.plot(original, label="original") plt.plot(shifted, label="shifted") plt.savefig("/dev/shm/timedelay.svg") plt.legend(loc=1) # Top right

## Easy zero crossing detection in Python using UliEngineering

In order to perform zero crossing detection in NumPy arrays you can use the UliEngineering library which provides an easy-to-use `zero_crossings`

function:

### How to install UliEngineering

UliEngineering is a Python 3 only library. Install using `pip`

:

sudo pip3 install -U UliEngineering

### Usage example:

from UliEngineering.SignalProcessing.Utils import zero_crossings # To generate test data from UliEngineering.SignalProcessing.Simulation import sine_wave # Generate test data: 10 Hz 10ksps, 1 second (= 1D 10000 values NumPy array) data = sine_wave(frequency=10.0, samplerate=10e3) # Prints: [0, 500, 1000, 1500, ...] print(zero_crossings(data))

Thanks to Jim Brissom on StackOverflow for the original solution!

## Iterating prime numbers in Python

### Problem:

Using Python you want to iterate through consecutive prime numbers

### Solution:

You can use gmpy2‘s `next_prime()`

to do this:

import gmpy2 def primes(start=2): n = start while True: n = gmpy2.next_prime(n) yield n # Usage example 1 for prime in primes(): print(prime) # Usage example 2 for i, prime in enumerate(primes()): print("The {}th prime number is {}".format(i, prime))

## Accurate calculation of PT100/PT1000 temperature from resistance

### TL;DR for impatient readers

PT100/PT1000 temperatures calculation suffers from accuracy issues for large sub-zero temperatures. `UliEngineering`

implements a polynomial-fit based algorithm to provide 58.6 \mu{\degree}C peak-error over the full defined temperature range from -200 {\degree}C to +850 °C.

Use this code snippet (replace `pt1000_`

by `pt100-`

to use PT100 coefficients) to compute an accurate temperature (in degrees celsius) e.g. for a resistane of 829.91 Ω of a PT1000 sensor.

from UliEngineering.Physics.RTD import pt1000_temperature # The following calls are equivalent and print -43.2316359463 print(pt1000_temperature("829.91 Ω")) print(pt1000_temperature(829.91))

You install the library (compatible to Python 3.2+) using

$ pip3 install -U UliEngineering