text:
# Ways of estimating phase shift

I have a system that's sampling two 50Hz sine waves at 1kHz. What's the best way to estimate the phase shift between the two?

The measurements are somewhat noisy and the sampling period is not very even. I know the sample times to a high resolution but can't control them to be evenly spaced. The phase shift needs to be estimated on chunks of data covering approximately one second.

The measurement noise makes the zero-crossing method unreliable; if the noise causes a measured signal to cross zero twice at one "real" zero-crossing then the algorithm breaks down somewhat.

At present I'm doing something like this (python code):

t, v, a = [NumPY arrays of samples]
p_inst = v * a
s = numpy.sign(p_inst)
s_avg = numpy.trapz(s, x=t) / (t[-1] - t[0])
shift_fraction = s_avg / 2 + 0.5
shift_angle = (1 - shift_fraction) * numpy.pi

Here s is the sign of the product of the two sine waves. Finding the average of this using numpy.trapz gives a value representing how often the two signals have the same sign, with 1 meaning "all of the time", -1 meaning "none of the time" and zero meaning "half of the time". This is converted to a fraction of the time the signals have the same sign (shift_fraction) and this is converted to a phase-shift angle.

The biggest problem with this is that it doesn't tell you whether the phase difference is leading or lagging.

Does anyone have a better method?

Edit For those suggesting finding the correlation peak, am I right in thinking that, with 20 samples per cycle, that method will necessarily limit the resolution of the estimate to 18 degrees? What would you think of instead calculating numpy.correlate(v, -a), finding the first zero crossing and then interpolating to estimate the actual zero point? I think that, for sinusoidal signals, the correlation should also be sinusoidal, and so the small angle approximation should make the interpolation reasonably good.

Can anyone compare the correlation and demodulation methods?

Further edit And, for those giving the demodulation method, am I right in thinking that if I don't care about the actual phase difference, just the power factor, then I can do this:

s = sin(t * 2 * pi * f)
c = sin(t * 2 * pi * f)
a_s_d = s * a
v_s_d = s * v
a_c_d = c * a
v_c_d = c * v
s_a = numpy.trapz(a_s_d, t)
c_a = numpy.trapz(a_c_d, t)
s_v = numpy.trapz(v_s_d, t)
c_v = numpy.trapz(v_c_d, t)
# cos ( a - b ) = cos a * cos b + sin a * sin b
power_factor = c_a * c_v + s_a * s_v

Or can I not rely on 0 < s_a, c_a, s_v, c_v < 1?

• Have you considered taking a cross-correlation of the two signals, using numpy's correlate or a library such as scipy (i.e., scipy.signal.signaltools.correlate(A, B))? The time-shift of the peak of the cross-correlation captures the phase difference. There's a good answer on StackOverflow that discusses this (it's more complex though). – Vladislav Martin Jun 15 '17 at 14:31
• cross correlation, but you will get peaks every 2*PI if signals are sine. – Marko Buršič Jun 15 '17 at 15:11
• Why is this tagged, "power factor correction"? Power factor correction has to do with high voltage lines and not sampling. – KingDuken Jun 15 '17 at 15:46
• @KingDuken, he's estimating the phase difference between voltage and current waveforms (based on the names of the variables in the pseudo-code). Probably for the purpose of estimating the power factor. – The Photon Jun 15 '17 at 16:11
• Filtering could help with the noise. Maybe zero crossing would then be more reliable. – mkeith Jun 16 '17 at 6:27

The easy to follow but cpu intensive way

If you want something which is obvious how it works and simple to understand, and you don't mind throwing lots of CPU cycles at it (sounds like you don't, if you're using python), then you could just fit sine waves to the two signals and read the phase out from the fit function.

You'd probably want to use scipy.optimise.curve_fit() to fit something like $V=A \times \sin \left(2 \pi f t+\phi \right)$ with $\phi$ as the only free parameter. This could converge on several different values of $\phi$, so take it modulo $2\pi$, then do the same for the current and subtract the two values of $\phi$, taking care to add or subtract $2\pi$ to get it in the meaningful range.

This will be much, much more CPU intensive than filtering, or doing realtime homodyne detection, but it's easy to follow, only a few lines long, and CPU cycles are cheap.

If you can't afford those cpu cycles

Then you might want to try generating two quaderature sine waves in software, then using them to do IQ demodulation of the two signals. Then calculate the phases and subtract. This will not be as easy to follow, but could be made very fast (even in python) and will be very, very accurate.

Take a chunk of data which is an integer number of cycles (assuming your 50Hz grid is nice and stable, just take 1s of data), and multiply it point-by-point with each of the generated sine waves. Then integrate each one using numpy.trapz() to get two scalars. Take these as the arguments to numpy.arctan2(), and voila, you have the phase of that signal relative to the generated sine wave. Do the same for the the second signal and subtract the two to get the difference. As above, add or subtract $2\pi$ to get it in the meaningful range.

• I ended up using your quadrature demodulation method - sorry it took so long to credit! – Tom Jul 10 '18 at 15:59

You do not mention it, but based on your variables it appears that you are attempting to measure the phase shift between current and voltage on a 50 Hz line.

To do proper phase measurements it is essential that the signal be clean and you specifically mentioned it is not. To get a clean signal, you can do low pass filtering either in an analog fashion prior to the ADC or digitally post the ADC. Naturally any phase shift from filtering must be equal for both signals or the introduced phase error must be predictable.

Two common means of implementing a digital low pass filter are IIR and FIR algorithms. Fortunately, the scipy library has them built in. For example, look at scipy ifilter functions. Be aware however that your data is considered nonuniform time data so you will first need to first resample your data using a spline or other appropriate interpolation algorithim.

Once the signals have gone through the LP filter then the most common technique for finding the phase difference is to find the peak of each signal, compute the time difference between peaks, and convert this to degrees or radians.

• Actually it'd be more accurate to say that I'm not very confident about the noise in the signals, rather than that I have any hard evidence to say it's bad. As for your method of calculating the phase shift, how is it better than the zero-crossing method? It seems to me that the small slope at the peak of a sinusoid would make estimating the peak worse than estimating the zero crossing. – Tom Jun 16 '17 at 8:23
• I have used it before quite successfully. With your high sample rate, it shouldn't be s problem. – Glenn W9IQ Jun 16 '17 at 11:08

I'm going to give two options. I don't know which will give better results.

1. Interpolate your data onto a regular time grid. Then use all the various well-known methods for phase estimation (zero-crossings, cross-correllation, ...) that you and other answers mentioned.

2. Construct a new sine and cosine array on your time grid. Estimate

$$A = \int_0^T f(t)\sin\left(2\pi f t\right){\rm d}t$$ and $$B = \int_0^T f(t)\cos\left(2\pi f t\right){\rm d}t$$

for each of your two signals, using trapz as you did to estimate average power.

Now the phase of each signals relative to the reference cosine can be estimated as $\tan^{-1}\left(\frac{A}{B}\right)$, or in Python numpy.arctan2(B, A). And the phase difference between the two signals can be obtained by subraction.

Caveat: There's a chance I made an error in the above that will give an error of $\pi/2$ in the individual signals' phase estimates. Also be careful of possible $2\pi$ wrap-around in the final subtraction. Check your results for reasonableness before using them.

One way is to use phase-sensitive detection. At 50 Hz with 1000 samples, you are getting 20 samples per cycle. Use 20-sample blocks, and subtract the sum of the second 10 samples from the sum of the first 10 samples. This will give you a value that is related to the phase relationship between the boundary of the block and the signal by the sine function (you will have to scale for amplitude-you can use an average of the P-P voltage). It will be at a positive peak when exactly in phase and a negative peak when you are 180 out. The approach removes any DC component. You can then average this value over many cycles to mitigate your noise. If you measure both signals with the same block timing, you get the phase relationship between each signal and your sampling system, which gets you an answer. Alternatively, you can build a loop that moves the block boundary (by shifting samples) and keeps the output at zero, but if you take this approach your resolution is 1/20th of a full cycle which is 18 degrees (assuming your sample times are fixed at 1 kHz).

• Does this method work with irregular sample intervals? – The Photon Jun 15 '17 at 16:15
• No. You are relying on two things: the amount of time for the positive samples and the amount of time for the inverted samples are exactly the same, and that they are at the frequency of interest. I think you can see that you are basically taking the data and subtracting one half of the samples from the other half, which should yield zero for any random stuff. Only data that has the characteristic of having a value of different polarity in each half cycle at your block rate will yield any output over the long haul. – John Birckhead Jun 15 '17 at 20:11

In the world of signal processing, phase-difference is commonly estimated with the cross-correlation of two signals; let's call them a(t) and b(t). Given a(t) and b(t) are similarly shaped wave-forms, cross-correlation should be a reliable approach.

Here's a Python script (part of which was inspired by this SO answer) that uses the correlate function found in the numpy library to calculate the phase shift between two sine waves, a(t) and b(t). You should be able to extrapolate the contents of this script for your own application:

import numpy as np

# Create the time axis (seconds)
num_samples = 10001
samples_per_second = 1000
freq_Hz = 50.0
t = np.linspace(0.0, ((num_samples - 1) / samples_per_second), num_samples)
# Create a sine wave, a(t), with a frequency of 1 Hz
a = np.sin((2.0 * np.pi) * freq_Hz * t)
# Create b(t), a (pi / 2.0) phase-shifted replica of a(t)
b_shift = (np.pi / 2.0)
b = np.sin((2.0 * np.pi) * freq_Hz * t + b_shift)

# Cross-correlate the signals, a(t) & b(t)
ab_corr = np.correlate(a, b, "full")
dt = np.linspace(-t[-1], t[-1], (2 * num_samples) - 1)
# Calculate time & phase shifts
t_shift_alt = (1.0 / samples_per_second) * ab_corr.argmax() - t[-1]
t_shift = dt[ab_corr.argmax()]
# Limit phase_shift to [-pi, pi]
phase_shift = ((2.0 * np.pi) * ((t_shift / (1.0 / freq_Hz)) % 1.0)) - np.pi

manual_t_shift = (b_shift / (2.0 * np.pi)) / freq_Hz

# Print out applied & calculated shifts
print ("Manual time shift: {}".format(manual_t_shift))            --> 0.005
print ("Alternate calculated time shift: {}".format(t_shift_alt)) --> 0.005000000000000782
print ("Calculated time shift: {}".format(t_shift))               --> 0.005000000000000782
print ("Manual phase shift: {}".format(b_shift))                  --> 1.5707963267948966, b(t) relative to a(t)
print ("Calculated phase shift: {}".format(phase_shift))          --> -1.570796326794651, a(t) relative to b(t)

When phase_shift is negative, then you know that the 1st signal input in correlate() is lagging behind the 2nd signal input. In our case, this means that a(t) is lagging behind b(t) by (pi / 2). When phase_shift is positive, then you know that the 1st signal input is leading ahead of the 2nd signal input.

As you can see, the cross-correlation of two signals can be simply used to detect the phase difference between two signals. The error vector between the actual & calculated time / phase shifts will remain low as long as your signals are linearly related, sampled at the a rate above their Nyquist Limit, and have a decent SNR.

Fun fact, you can still use the cross-correlation for this application even if the signals are a-periodic! Take caution when taking the cross-correlation between two signals that are related non-linearly.

• Thanks for a detailed answer. Could you have a look at the edit to the question and respond? In your example, you have 1000 samples per cycle, while I have 20 - doesn't this rather limit the resolution of this method? – Tom Jun 16 '17 at 8:26
• Also, isn't generating an array of [-t:t/num_samples:t] and using an index lookup rather a wasteful way of calculating the time shift? What's wrong with 2 * t * (ab_corr.argmax() - num_samples/2) / num_samples? – Tom Jun 16 '17 at 8:29
• @Tom See my results; they match the (pi / 2) phase shift we expect well. I don't believe that cross-correlation is limited by the sampling rate, as long as you're obeying Nyquist's Theorem (i.e., not under-sampling your signals). Regarding the dt array; I've included an alternative method for calculating time shift that avoids an index look-up (producing same time shift as with the dt index look-up method). One thing that is wrong with your suggestion is that it will also generate an array (remember, t is an array)... – Vladislav Martin Jun 16 '17 at 16:48

If your only problem with the measurement is the ambiguity, you can resolve it by going back to the data, using some method to harmonize zero crossings too close to one another (take their mean arrival time, maybe), and then calculating the angle using the zero crossing method. The result you want will be whichever of the two possible phase shifts is closer to the zero crossing result.

While I'm on the subject, one method to find the phase angle that I've never tried, but should work, which also gives ambiguous results, is to measure the voltage of each sine wave, then measure the voltage difference between them. That would give you the 3 sides of a triangle and you can find the angles with law of sines/law of cosines.


text:
Negative Quarks
 Characteristics of negative quarks total angular momentum zero charge zero lepton number ±1/8 baryon number zero temperature 676 (℃) internal energy -33 (MeV)

A seed-aggregate composed from one negative seed and one ordinary seed is called an ordinary negative quark. It is represented symbolically using the letter e. Here is a mathematical statement of the definition along with an icon to illustrate the quark.

 $\sf{ e \equiv \left\{ E, O \right\} }$

You can usually click on quark icons to come back to a page like this for reference. The negative anti-quark is defined as the union of one negative seed with an

 $\sf{ \bar{e} \equiv \left\{ E, \overline{O} \right\} }$

These definitions can also be stated using seed icons as

 + →
 + →
Sensory interpretation: Ordinary negative quarks represent yellow Anaxagorean sensations on the right-side, and negative anti-quarks are objectified from yellow perceptions felt on the left.

 Summary
 Noun Definition Negative Quark $\sf{e} \equiv \{ \sf{E}, \sf{O} \}$ 3-44
 Noun Definition Negative Anti-quark $\sf{\overline{e}} \equiv \{ \sf{E}, \sf{\overline{O}} \}$ 3-43

Note to editors: the layout for this page is determined by this template.

page revision: 12, last edited: 24 Aug 2019 21:48

format: plain


text:
zbMATH — the first resource for mathematics

Examples
 Geometry Search for the term Geometry in any field. Queries are case-independent. Funct* Wildcard queries are specified by * (e.g. functions, functorial, etc.). Otherwise the search is exact. "Topological group" Phrases (multi-words) should be set in "straight quotation marks". au: Bourbaki & ti: Algebra Search for author and title. The and-operator & is default and can be omitted. Chebyshev | Tschebyscheff The or-operator | allows to search for Chebyshev or Tschebyscheff. "Quasi* map*" py: 1989 The resulting documents have publication year 1989. so: Eur* J* Mat* Soc* cc: 14 Search for publications in a particular source with a Mathematics Subject Classification code (cc) in 14. "Partial diff* eq*" ! elliptic The not-operator ! eliminates all results containing the word elliptic. dt: b & au: Hilbert The document type is set to books; alternatively: j for journal articles, a for book articles. py: 2000-2015 cc: (94A | 11T) Number ranges are accepted. Terms can be grouped within (parentheses). la: chinese Find documents in a given language. ISO 639-1 language codes can also be used.

Operators
 a & b logic and a | b logic or !ab logic not abc* right wildcard "ab c" phrase (ab c) parentheses
Fields
 any anywhere an internal document identifier au author, editor ai internal author identifier ti title la language so source ab review, abstract py publication year rv reviewer cc MSC code ut uncontrolled term dt document type (j: journal article; b: book; a: book article)

Found 2931 documents (Results 1–100)

Münster: Univ. Münster, Mathematisch-Naturwissenschaftliche Fakultät, Fachbereich Mathematik und Informatik (Diss.). vi, 92 p. (2013).
Kielanowski, Piotr (ed.) et al., Geometric methods in physics. XXX workshop, Białowieża, Poland, June 26 – July 2, 2011. Selected papers based on the presentations at the workshop. Basel: Birkhäuser (ISBN 978-3-0348-0447-9/hbk; 978-3-0348-0448-6/ebook). Trends in Mathematics, 303-310 (2013).
MSC:  46L05 46L55
Benguria, Rafael (ed.) et al., Spectral analysis of quantum Hamiltonians. Spectral Days 2010. Selected papers of the conference, Santiago, Chile, September 2010. Basel: Birkhäuser (ISBN 978-3-0348-0413-4/hbk; 978-3-0348-0414-1/ebook). Operator Theory: Advances and Applications 224, 11-24 (2012).
Pusz, Wiesław (ed.) et al., Operator algebras and quantum groups. Proceedings of the conference, Warsaw, Poland, September 19–23, 2011. Warszawa: Polish Academy of Sciences, Institute of Mathematics (ISBN 978-83-86806-17-1/pbk). Banach Center Publications 98, 433-453 (2012).
Pusz, Wiesław (ed.) et al., Operator algebras and quantum groups. Proceedings of the conference, Warsaw, Poland, September 19–23, 2011. Warszawa: Polish Academy of Sciences, Institute of Mathematics (ISBN 978-83-86806-17-1/pbk). Banach Center Publications 98, 107-198 (2012).
Banach Center Publications 98. Warszawa: Polish Academy of Sciences, Institute of Mathematics (ISBN 978-83-86806-17-1/pbk). 486 p. EUR 50.00 (2012).
Bożejko, Marek (ed.) et al., Noncommutative harmonic analysis with applications to probability. III: Papers presented at the 13th workshop, Bȩdlewo, Poland, July 11–17, 2010. Warszawa: Polish Academy of Sciences, Institute of Mathematics (ISBN 978-83-86806-15-7/pbk). Banach Center Publications 96, 299-317 (2012).
Bożejko, Marek (ed.) et al., Noncommutative harmonic analysis with applications to probability. III: Papers presented at the 13th workshop, Bȩdlewo, Poland, July 11–17, 2010. Warszawa: Polish Academy of Sciences, Institute of Mathematics (ISBN 978-83-86806-15-7/pbk). Banach Center Publications 96, 159-174 (2012).
Mathematik. München: Dr. Hut; Darmstadt: Univ. Darmstadt, Fachbereich Mathematik (Diss.) (ISBN 978-3-8439-0517-6/hbk). xxi, 171 p. EUR 60.00/hbk (2012).

all top 5
all top 5
all top 5


text:
# Introduction #

In my last blog entry here, we looked at improving the performance of the SPARKNaCl library, while preserving all the proof results.

At the end of that blog, I speculated a few ideas for further improvement that might bear a bit more fruit. This blog documents how three of those ideas turned out and the results that I got.

Don’t worry… this one is much shorter than the last…

## Starting point #

At the end of the last blog, we had the following performance results for a single NaCl “Sign” operation, running on a 32-bit RISC-V board. This time, I will also include result for worst-case stack usage (WCSU, generated by the GNATstack tool) and code size using the RM32IM (not compressed) ISA (generated by the riscv32-elf-size tool) for the entire SPARKNaCl library (so not including the effect of the linker removing unused sections). As before, execution times are in Millions of CPU cycles, while WCSU and code size are in bytes:

Level      Time    WCSU     Size
-O0    181.74    6688    49162
-O1     69.96    4824    30492
-O2     57.82    5504    31180

## The big idea #

TweetNaCl stores “big” (256 bit) integers as an array of sixteen 64-bit signed integers. This type is called “GF” in the sources. A “Normal” GF has all the digits or “limbs” in the range 0 .. 65535. Signed 64-bit “limbs” are needed to store the intermediate results of “+”, “-” and “*” operations, which can overflow and/or become negative before being “normalized” modulo 2**255-19 ready for the next operation.

1. Most of the time, a “Normal_GF” value occupies 128 bytes of space, but only really contains 32 bytes of useful information, since the most significant 6 bytes of each limb are always zero. This takes up lots of space on the stack, and copying/assignment of them is inefficient. See our earlier comments on functional programming style and “Return Slot Optimization” for how this kind of copying overhead can hit performance.
2. I also realized that the “Carry and Reduce” (CAR) functions were all using 64-bit signed arithmetic all the time. This is really slow on a 32-bit machine like our RV32 target.

So…here’s the idea: what if we store a Normal_GF as an array of 32-bit values, and only “expand out” into 64-bit arithmetic for the “*” operation?

To be sure this will work, we’d have to know exact bounds on the types produced by the “*” operation at each phase of the normalization process. Good news: we already have that information embodied in the subtype declarations for those functions. For example, the first “CAR” function that follows “*” has the declaration:

function Product_To_Seminormal (X : in Product_GF) return Seminormal_GF
with Pure_Function,
Global => null;

so let’s look at those type declarations too. These fragments are taken from sparknacl.ads at the Faster_Tweet” tag so they correspond to what we had at the end of the last blog.

--  "LM"   = "Limb Modulus"
--  "LMM1" = "Limb Modulus Minus 1"
LM   : constant := 65536;
LMM1 : constant := 65535;

--  "R2256" = "Remainder of 2**256 (modulo 2**255-19)"
R2256 : constant := 38;

--  "Maximum GF Limb Coefficient"
MGFLC : constant := (R2256 * 15) + 1;

--  "Maximum GF Limb Product"
MGFLP : constant := LMM1 * LMM1;

subtype Product_GF is GF
with Dynamic_Predicate =>
(for all I in Index_16 =>
Product_GF (I) >= 0 and
Product_GF (I) <=
(MGFLC - 37 * GF_Any_Limb (I)) * MGFLP);

----------------------------------------------------------------------
--  A "Seminormal GF" is the result of applying a single
--  normalization step to a Product_GF
--
--  Least Significant Limb ("LSL") of a Seminormal GF.
--  LSL is initially normalized to 0 .. 65535, but gets
--  R2256 * Carry added to it, where Carry is (Limb 15 / 65536)
--  The upper-bound on Limb 15 is given by substituting I = 14
--  into the Dynamic_Predicate above, so
--    (MGFLC - 37 * 14) * MGFLP = 53 * MGFLP
--  See the body of Product_To_Seminormal for the full
--  proof of this upper-bound
subtype Seminormal_GF_LSL is I64
range 0 .. (LMM1 + R2256 * ((53 * MGFLP) / LM));

--  Limbs 1 though 15 are in 0 .. 65535, but the
--  Least Significant Limb 0 is in GF_Seminormal_Product_LSL
subtype Seminormal_GF is GF
with Dynamic_Predicate =>
(Seminormal_GF (0) in Seminormal_GF_LSL and
(for all I in Index_16 range 1 .. 15 =>
Seminormal_GF (I) in GF_Normal_Limb));

The first important number to check is the worst-case upper bound on all the limbs of Product_GF. This occurs for limb 0, and has upper bound (MGFLC * MGFLP) which has a value a bit below 2**42, so we definitely can’t store one of those in 32 bits.

The second important number is Seminormal_GF_LSL’Last, which evaluates to a bit less than 2**27, so we can store a Seminormal_GF using 32-bit limbs. That also means that a Nearlynormal_GF (which is the type returned by the second “Seminormal_To_Nearlynormal” normalization step) can be represented with 32-bit limbs.

The upshot is that the “*” operation needs to “expand out” its parameters into 64-bit before it does all the sums, the first normalization (Product_To_Seminormal) has to be done in 64-bit arithmeric, but can return an array of 32-bit values, then the second and third normalization functions can all be done in 32-bit maths.

Does this make much of a difference? I did a quick-and-dirty 32-bit implementation of “CAR” for comparison. On RV32 the difference is startling - 133 Instructions for a CAR function in 32-bit arithmetic, versus 321 Instructions for the 64-bit variant, so that’s a respectable improvement and doesn’t take into account the benefit of saving memory, faster assignments and so on.

## Results (Git Tag: GF32) #

I made the changes suggested above, and used GNATprove to make sure I hadn’t made a mess of it. I didn’t test the code at all until the proof was complete and automated as before. Here are the new numbers:

Level      Time    WCSU     Size
-O0    172.59    4416    47092
-O1     65.14    3595    27448
-O2     34.94    3536    27668

So we’ve saved about 2k of stack space and 2k of code size across the board. Performance shows a reasonably improvement at -O0 and -O1, and a massive jump at -O2 with an improvement of nearly 23 Million cycles.

Getting this far took about 5 hours at the keyboard, including code modifications, re-generating the proofs, a single (successful) run of the regression tests, and getting the performance numbers above.

# Going further… (Git Tag: GF16) #

So far so good. Those results gave me the confidence to take the next logical step.

The most compact representation possible for a “Normal_GF” is to store each limb (with range 0 .. 65535) as an unsigned 16-bit value, so the whole thing fits in 32 bytes with no wasted space at all.

Using that representation doesn’t change the maths in the operations at all - we still need 64-bit math for “*” and the first CAR, and 32-bit math for the second and third CAR and for “+” and “-“, but it will halve the cost of an assignment (again), reduce memory usage and (possibly) improve data cache hit rate, so probably worth a go…

Having done it once, this change was relatively simple - about 2 hours time - although some care was needed to make sure all the proofs were still OK.

Results as follows:

Level      Time    WCSU     Size
-O0    167.48    3296    49684
-O1     32.60    2560    27460
-O2     27.43    2592    27868

So…code size is about the same, and we save another 1k of stack space across the board. Performance improves a bit at -O0, gets a massive boost at -O1, and another noticeable improvement at -O2.

Pleasantly enough, these numbers are also well ahead of the best results that we’ve seen for TweetNaCl.

I could transfer these changes to the TweetNaCl C sources I suppose, but I really don’t have the confidence to do such a subtle bunch of changes without GNATprove backing me up. Perhaps the maintainers of WireGuard will decide to pick up these ideas?

# What about -O3? (Git Tag: O3_No_Unroll) #

At the end of the last blog, I promised I’d have a look at at results with optimization level 3.

The main differences between -O2 and -O3 are that -O3 automatically inlines more subprogram calls, and automatically unrolls loops (whether you want it to or not…) These can be a double-edged sword - possibly faster through fewer branches, but possibly slower through increased code size and reduced hit rate for the instruction cache.

Using exactly the same code as we left off with, we can add a fourth row to the table:

Level      Time    WCSU     Size
-O0    167.48    3296    49684
-O1     32.60    2560    27460
-O2     27.43    2592    27868
-O3     49.87    2624    48656

Having manually unrolled the most time-critical loops already, one thing we could do is to use GNAT’s pragma Loop_Optimize (No_Unroll) to explicitly tell the compiler not to unroll all the other loops in the library. With those pragmas in place, we get:

Level      Time    WCSU     Size
-O0    167.48    3296    49684
-O1     32.60    2560    27460
-O2     27.43    2592    27868
-O3     25.80    2736    34404

so a slight improvement in performance for a small increase in stack usage, and a few more kilobytes of code size over -O2 (owing to inlining no doubt).

# Feeling the Squeeze (Part 2…) (Git Tag: RV32IMC_A4) #

And finally… I promised I’d return to looking at the performance of SPARKNaCl using the RISV-V Compressed Instruction Set. Previously, we’d seen a significant performance improvement by switching from the compressed instructions (known as “RV32IMC”) to uncompressed (known as “RV32IM”) owing to a performance penalty in the E31 core for branching to a mis-aligned instruction.

Eric Botcazou at AdaCore suggested having a look at the compiler’s options to force instruction alignment to see if we could put things right.

Section 3.11 of the GCC manual says there are four options to control the alignment of function entry points, branch labels, jump labels and loop entry points. With the RV32IMC instruction set, these all default to “2” (so all instructions appear at an address which is an exact multiple of 2), giving maximum code density, but (with our E31 core) incurring a performance penalty when we branch to a 32-bit instruction that is not aligned on a 4-byte boundary.

To see how bad the penalty is, let’s look at performance, worst-case stack usage, and code size using RV32IMC and the default alignment, based on the code as above. For fun, I have also added a build using the option “-Os” which means “optimize for smallest code size” just to see what happens.

Level      Time    WCSU     Size
-O0    170.53    3296    40430
-O1     32.51    2592    20230
-O2     62.72    2608    20972
-O3     25.73    2736    26000
-Os     30.54    2512    18280

So… we’re a bit slower at -O0, a bit faster at -O1 and -O3, but -O2 suffers some sort of terrible penalty. Worst-case stack usage is mostly unchanged, while code size comes down by about 7k across the board, with -Os winning with the whole library fitting in slightly over 18k.

Now we’ll see what happens when I set all four alignment options to “4”. The Makefile in perf/tsign.mk selects this by setting the SPARKNACL_TARGET_ARCH project variable to “rv32imc_a4”. With that, we hope to keep the code size of compressed instructions, but avoid the performance penalty of mis-aligned branches. The results are:

Level      Time    WCSU     Size
-O0    174.47    3296    41360
-O1     32.54    2592    20908
-O2     41.10    2608    21672
-O3     25.69    2736    26828
-Os     30.52    2512    18616

In terms of performance, -O0 and -O1 are a bit slower. -O2 regains some credibility, while -O3 and -Os are a bit faster than before. Worst-case stack usage is unchanged, while code size goes up by about 800 bytes at all levels owing to the need to “pad” basic blocks with “no-op” instructions to regain instruction alignment.

# Conclusions #

In summary:

Fastest performance: 25.69 Million cycles for a “Sign” at -O3 and RV32IMC_A4

Smallest worst-case stack usage: 2512 bytes at -Os (and regardless of instruction set)

Smallest code size: 18280 bytes at -Os and RV32IMC

From the initial starting point we’ve seen the performance of a “Sign” operation come down from 93.86 to 27.43 million cycles at -O2 - an improvement of about 3.4 times, without breaking the proof. I think this last point cannot be overstated… years of programming (and PSP training) have taught me what a brilliant programmer I’m not. I just don’t know if I could make subtle changes without a tool like GNATprove backing me up (and saving me from acute embarrassment).

Of course, there’s more to come. Firstly, we could attempt to reduce the number of “CAR” calls from 3 to 2 following “*”. Benoit Viguier’s work (see his recent paper on the subject) shows that this is possible (and provable!) but relies on some very deep work with the Coq theorem prover that I can’t really get my head around at the moment…

Secondly, we’ll report more performance figures as and when GNAT adopts more advanced Return Slot Optimization for SPARK, as we mentioned last time.

Then there’s GNAT-LLVM, and … I guess it all depends on when lockdown ends and if I get any “proper” work to do… we’ll see…