This notebook presents an algorithm to extract frequency, amplitude and phase from a sampled sine wave signal. The performance of the algorithm is evaluated off-line on a previously captured signal.

Sine fitting algorithm

The algorithm presented in [1][2] uses a method of least squares to fit the sine wave to data points and thus obtaing the parameters of the sine wave with lowest error between reconstructed and measured signal.

The signal can be modeled with the following equation:

$y[n] = M * sin(\omega t_n + \phi) + C = A cos(\omega t_n) + B sin(\omega t_n) + C$

where $t_n = n*Ts$ and $M$ and $\phi$ can be expressed as:

$M = \sqrt{A^2 + B^2}$, $\phi = atan(B/A)$

When frequency is unknown, we can start with an initial estimate of the frequency and add a fourth parameter (difference of frequency - $\Delta \omega$) to the set of unknowns.

We can iterativelly calculate the frequency difference between the measured frequency and our estimation and then add the frequency difference to our previous estimate. Instead of using assuming fixed frequency $\omega$, we assume that a small frequency difference $\Delta \omega$ will be present in our model.

We can then develop our initial model to also incrporate $\Delta \omega$ in its parameters:

$y(n, w + \Delta \omega) \approx y(n,w) + \Delta \omega \frac{\delta y(n,w)}{\delta \omega}$

$y(n, w + \Delta \omega) \approx A cos(\omega t_n) + B sin(\omega t_n) + C + \Delta \omega (-A sin(wt) t + B cos (wt)t)$

The equations can be written in matrix form:

$\mathbf{D} \cdot \mathbf{x} = \mathbf{y}$

$\mathbf{D} = \left[ \begin{array}{cccc} A cos(\omega 0 T_s) & B sin(\omega 0 T_s) & 1 & -A sin(w0 T_s) 0 T_s + B cos (w0 T_s)0 T_s \\ A cos(\omega 1 T_s) & B sin(\omega 1 T_s) & 1 & -A sin(w0 T_s) 0 T_s + B cos (w0 T_s)0 T_s \\ \ldots & \ldots & \ldots & \ldots \\ A cos(\omega (N-1) T_s) & B sin(\omega (N-1) T_s) & 1 & -A sin(w(N-1) T_s) (N-1) T_s + B cos (w(N-1) T_s)(N-1) T_s \\ \end{array} \right]$

$\mathbf{x} = \left[ \begin{array}{c} A \\ B \\ C \\ \Delta \omega \\ \end{array} \right]$

and $\mathbf{y}$ is the vector of measurements of the unknown signal.

The optimal solution of overdetermined set of equation can thus be obtained by solving the equation:

$\hat{x} = (\mathbf{D}^T \cdot \mathbf{D})^{-1} \cdot \mathbf{D}^T \cdot \mathbf{y}$

Some comments regarding the implementation of the algorithm in an FPGA or in a processor:

1. The algorithm needs an initial estimation of frequency, amplitude and phase. When the algorithm is being used repeatedly, the results from previous run can be used as a starting point for the next run. Since the difference in parameters between consecutive runs should be small, this will ensure faster convergence.

2. The iterative nature of the algorithm, or more precisely the iterative nature of the frequency matching part of the algorithm, could in some unexpected situation cause problems with the convergence and thus with the performance of the entire system. It is recommended that the number of iterations is bound in the real-time application. If the algorithm does not converge in a predefined number of loops, an error should be signaled, post-mortem data should be collected and the system should be stopped.

Demonstration of the algorithm

In the following example, a sine wave signal with "unknown" amplitude, frequency and phase is generated (function gen_sine). Some noise is added to the "unknown" signal to simulate the noise obtained from the measurement of the real signals. The noise of amplitude of 0.01 with the signal of amplitude of 2.23 results in SNR of -46 dB ($20 log_{10}(0.01/2.23)$. The "unknown" signal is then feed to the algorithm which calculates the 4 unknown variables of the signal. Using the calculated parameters, a new sine wave is reconstructed and shown together with the "unknown" signal.

The results of the algorithm are printed and the plots of original signal and reconstructed signal are shown. The results are really promising, and the algorithm is able to match the original waveform in spite of added noise.

The calculated amplitude error is 250 ppm, phase error is 250 ppm, frequency error is 2ppm and offset error is 70 ppm.

Measurements with CAEN ELS TetrAMM

After the performance of proposed algorithm has been verified on an artificial data it is time to evaluate the algorithm on the real world data. The measurements were done with CAEN ELS TetrAmm and the data was stored in csv files for later (off-line) elaboration. The TetrAmm front end was modified to allow voltage inputs of $\pm 1.25V$ with analog bandwidth limited at 15 kHz. A Python script was used to capture 100 samples at 100 kHz from TetrAMM (using the FASTNAQ command) and the data was stored in csv files.

A Tektronix arbitrary function generator was used to produce sine waveform with different frequency and between-channel phase. Bellow is a list of all phases and frequencies tested. In each point 5 measurement were performed, producing a total of 130 files.

Captured data

Shown bellow is a Python script which opens one of csv files, parses the data in it and plots the captured data. It can be noted that because we are sampling roughly at 20-times the sine wave still looks a little bit jagged.

Algorithm in action

Now we can run our algorithm on previously captured data. The param_est function is called for each of two channels and the calculated parameters are printed out. The phase difference between the two channels is also calculated. Lastly, a sine wave is recreated (reconstructed) from the calculated parameters and a plot with both captured and reconstructed sine wave is shown.

After we have confirmed that we are able to successfully extract amplitude, frequency and phase from a sine captured with TetrAMM unit, we can run the same algorithm on each of the 130 measurements. This would allows us to better understand performance of the algorithm in various condition.

A table with all measurements (frequency from both signals, and phase) are shown in the table bellow.

Measurement errors

The previous script performed the calculation of the algorithm on various data sets with different frequencies and phase between channels. To understand the precision of algorithm, we can compare the calculated (measured) values to the original values - the values which were set on function generator. For each data set we can calculate the difference between the parameters and plot the difference.

From the two Figures above, we can see that using the previously described algorithm and performing the measurements with CAEN ELS TetrAMM unit with sampling at 100 kHz we can obtain the frequency of signal with precision of 3 Hz and phase between two sine waves with precision of 0.025 deg.

Implementation of the algorithm in FPGA

We have presented and evaluated the sine-fitting algorithm which is able to measure frequency, phase and amplitude of signal in high-level language (Python) using external library (Numpy). Because the algorithm would be used in FPGA signal-processing chain, some implementational details should be also considered.

Execution time

Let's estimate the execution time of the algorithm. The single columns in D matrix can be calculated on-the-fly when the matrices D.T and D are multiplied together. Since the floating-point multipliers in FPGA are pipelined, the calculation of $D^T \cdot D$ would take approximatelly 1600 clock cycles (100 clock cycles for each element in the resulting 4x4 matrix).

The obtain the execution time for matrix inverse Vivado High Level Synthesis (HLS) tool was used. Vivado HLS compiles the C and C++ programs (only subset of the language is supported) to Verilog, which can then be compiled to FPGA bitstream.

The following is a source code to calculated 4x4 matrix inverse using Cholesky decomposition.

#include <hls_linear_algebra.h>

void matrix_inverse(
	float mat_input[4][4],
	float mat_output[4][4]
){
#pragma HLS INTERFACE s_axilite port=mat_input
#pragma HLS INTERFACE s_axilite port=mat_output

	int success;
	hls::cholesky_inverse<4,float,float>(mat_input, mat_output, success);
}

For the above algorithm, Vivado HLS produces this report:

We can see that a matrix can be inverted in roughly 1000 clock cycles.

To multiply the 4x4 matrix with $D^T$ matrix of dimensions 4x100 we can use 4 multipliers and 3 binary adders to get one element of the matrix $(\mathbf{D}^T \cdot \mathbf{D})^{-1} \cdot \mathbf{D}^T$ in a clock cycle, requiring a total of 400 clock cycles. The final multiplication between matrix $(\mathbf{D}^T \cdot \mathbf{D})^{-1} \cdot \mathbf{D}^T$ and vector $\mathbf{y}$ would again 400 clock cycles (using only one multiplier).

A rought estimation of execution time of a single iteration of algorithm is therefore 1600+1000+400+400 = 3400 clock cycles. With operating clock frequency of 200 MHz, each iteration would take roughly 20 us. The execution time can be significantly reduced by performing several operations in parallel.

FPGA resource usage

There are several matrices which must be held in memory when performing the calculations. The largest matrices would have dimensions 100x4 or 100x4, with total of 400 elements. If each element is stored as a single-precision float the total memory usage would be 1600 bytes (12800 bits). The FPGA in CAEN ELS BEST Central Unit (Altera Cyclone V 5CGTFD9E5F35C7N [3]) has a total of 12200 kbits of memory, which is more than enough to store a couple of matrices

The single-precision floating-point multiplier would require a little less than 200 ALMs and 4 DSP blocks in the FPGA, while the adder would require a little less than 500 ALMs [4]. With 110000 ALMs and 684 DSP blocks in Cyclone V GT FPGA, there is also more than enough space to implement the algorithm in widely parallelized form.

References

[1] IEEE Standard for Digitizing Waveform Recorders - Redline," in IEEE Std 1057-2007(Revision of IEEE Std 1057-1994) - Redline , vol., no., pp.1-210, April 18 2008

[2] P. Handel, Evaluation Of A Standardized Sine Wave Fit Algorithm, http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.35.3486

[3] Cyclone V Device Overview, https://www.altera.com/content/dam/altera-www/global/en_US/pdfs/literature/hb/cyclone-v/cv_51001.pdf

[4] Floating-Point IP Cores User Guide, https://www.altera.com/content/dam/altera-www/global/en_US/pdfs/literature/ug/ug_altfp_mfug.pdf