From 7ce5b82ff0b00ecb61e29d58df632b892e099b30 Mon Sep 17 00:00:00 2001
From: Joris van Rantwijk <joris@jorisvr.nl>
Date: Sat, 16 Apr 2016 09:09:28 +0200
Subject: [PATCH] * eval_sine_quality.py: Add Python program to evaluate sine
 waveform.

---
 tools/eval_sine_quality.py | 153 +++++++++++++++++++++++++++++++++++++
 1 file changed, 153 insertions(+)
 create mode 100644 tools/eval_sine_quality.py

diff --git a/tools/eval_sine_quality.py b/tools/eval_sine_quality.py
new file mode 100644
index 0000000..a8a9bcf
--- /dev/null
+++ b/tools/eval_sine_quality.py
@@ -0,0 +1,153 @@
+#!/usr/bin/python
+
+"""
+Evaluate quality of generated sine/cosine waveform.
+
+Reads output file from VHDL testbench and reports
+key figures for the quality of the sine wave.
+
+Usage:
+  python eval_sine_quality.py datafile.dat
+"""
+
+
+from __future__ import print_function
+import sys
+import numpy
+
+
+def read_data(fname):
+
+    n = 0
+    data = numpy.zeros((1024, 2), dtype=numpy.int64)
+
+    with open(fname, 'r') as f:
+        for s in f:
+            if n == data.shape[0]:
+                data = numpy.resize(data, (2*n, 2))
+            w = s.split()
+            assert len(w) == 2
+            if len(w) != 2:
+                raise ValueError("Expecting two columns in file")
+            data[n,0] = int(w[0])
+            data[n,1] = int(w[1])
+            n += 1
+
+    return data[:n,:]
+
+
+def eval_sine_quality(data):
+
+    assert len(data.shape) == 2
+    assert data.shape[1] == 2
+    n = data.shape[0]
+    assert n >= 4 and n % 4 == 0
+
+    # Extract sine and cosine colums.
+    dsin = numpy.copy(data[:,0])
+    dcos = data[:,1]
+    del data
+
+    # Check that cosine is an exact phase-shifted version of sine.
+    dcos[0:3*n//4] -= dsin[n//4:]
+    dcos[3*n//4:]  -= dsin[0:n//4]
+    tmin = numpy.amin(dcos)
+    tmax = numpy.amax(dcos)
+    del dcos
+
+    if tmin == 0 and tmax == 0:
+        print('cos(x) == sin(x+pi/2) exactly')
+    else:
+        print('cos(x) == sin(x+pi/2) + (%d .. %d)' % (tmin, tmax))
+
+    # Check that 180 degree phase shift is equivalent to sign negation.
+    t = dsin[0:n//2] + dsin[n//2:]
+    tmin = numpy.amin(t)
+    tmax = numpy.amax(t)
+    del t
+
+    if tmin == 0 and tmax == 0:
+        print('sin(x) == - sin(x+pi) exactly')
+    else:
+        print('sin(x) == - sin(x+pi) + (%d .. %d)' % (tmin, tmax))
+    print()
+
+    # Determine offset (mean value of waveform).
+    offs = numpy.mean(dsin)
+    print('offset =        %20.12f lsb' % offs)
+
+    # Determine amplitude and phase.
+    tref = numpy.sin(2 * numpy.pi / n * numpy.arange(n))
+    asin = numpy.sum(dsin * tref) * 2.0 / n
+    del tref
+    tref = numpy.cos(2 * numpy.pi / n * numpy.arange(n))
+    acos = numpy.sum(dsin * tref) * 2.0 / n
+    del tref
+
+    ampl   = numpy.sqrt(asin**2 + acos**2)
+    phase  = numpy.arctan2(acos, asin)
+
+    print('amplitude =     %20.12f lsb' % ampl)
+    print('phase offset =  %20.12f rad' % phase)
+    print()
+
+    # Determine peak and rms deviation.
+    tref = ampl * numpy.sin(2 * numpy.pi / n * numpy.arange(n) + phase)
+    terr = dsin - tref
+    del tref
+
+    peakerr = numpy.amax(numpy.abs(terr))
+    rmserr  = numpy.std(terr)
+    del terr
+
+    print('peak error =    %20.12f lsb' % peakerr)
+    print('rms error =     %20.12f lsb rms' % rmserr)
+
+    # Calculate SNR and effective number of bits.
+    sinad = 20 * numpy.log10(ampl * numpy.sqrt(0.5) / rmserr)
+    print('SINAD =         %12.4f dB' % sinad)
+    print('ENOB =          %12.4f bits' % ((sinad - 1.76) / 6.02))
+
+    # Determine spurious-free dynamic range.
+    q = numpy.fft.rfft(dsin)
+    tampl = numpy.abs(q[1])
+    tspur = numpy.amax(numpy.abs(q[2:]))
+    del q
+    sfdr = 20 * numpy.log10(tampl / tspur)
+
+    print('SFDR =          %12.4f dB' % sfdr)
+    print()
+
+
+def main():
+
+    if len(sys.argv) != 2:
+        print(__doc__, file=sys.stderr)
+        print("ERROR: Invalid/missing command line arguments", file=sys.stderr)
+        sys.exit(1)
+
+    fname = sys.argv[1]
+
+    # Read data from file.
+    print("reading", fname, "...")
+    data = read_data(fname)
+
+    print("got array", data.shape)
+    print()
+
+    # Check array shape.
+    if len(data.shape) != 2 or data.shape[1] != 2:
+        print("ERROR: Expected array of shape (N, 2)", file=sys.stderr)
+        sys.exit(1)
+
+    # Check number of samples.
+    if data.shape[0] < 4 or (data.shape[0] & (data.shape[0] - 1)) != 0:
+        print("ERROR: Expected power-of-two record length", file=sys.stderr)
+        sys.exit(1)
+ 
+    eval_sine_quality(data)
+
+
+if __name__ == '__main__':
+    main()
+