Fix CAZAC sequence

commpy/sequences.py

@@ -15,6 +15,7 @@
 """
 __all__ = ['pnsequence', 'zcsequence']
 
+import numpy as np
 from numpy import empty, exp, pi, arange, int8, fromiter, sum
 
 def pnsequence(pn_order, pn_seed, pn_mask, seq_length):
@@ -72,23 +73,40 @@ def pnsequence(pn_order, pn_seed, pn_mask, seq_length):
 
     return pnseq
 
-def zcsequence(u, seq_length):
+def zcsequence(u, seq_length, q=0):
     """
     Generate a Zadoff-Chu (ZC) sequence.
 
     Parameters
     ----------
     u : int
-        Root index of the the ZC sequence.
+        Root index of the the ZC sequence: u>0.
 
     seq_length : int
-        Length of the sequence to be generated. Usually a prime number.
+        Length of the sequence to be generated. Usually a prime number:
+        u<seq_length, greatest-common-denominator(u,seq_length)=1.
+
+    q : int
+        Cyclic shift of the sequence (default 0).
 
     Returns
     -------
     zcseq : 1D ndarray of complex floats
         ZC sequence generated.
     """
-    zcseq = exp((-1j * pi * u * arange(seq_length) * (arange(seq_length)+1)) / seq_length)
+
+    for el in [u,seq_length,q]:
+        if not float(el).is_integer():
+            raise ValueError('{} is not an integer'.format(el))
+    if u<=0:
+        raise ValueError('u is not stricly positive')
+    if u>=seq_length:
+        raise ValueError('u is not stricly smaller than seq_length')
+    if np.gcd(u,seq_length)!=1:
+        raise ValueError('the greatest common denominator of u and seq_length is not 1')
+
+    cf = seq_length%2
+    n = np.arange(seq_length)
+    zcseq = np.exp( -1j * np.pi * u * n * (n+cf+2.*q) / seq_length)
 
     return zcseq

commpy/tests/test_sequences.py

@@ -1,16 +1,17 @@
 # Authors: CommPy contributors
 # License: BSD 3-Clause
 
+import numpy as np
 from numpy import array
-from numpy.testing import run_module_suite, assert_raises, assert_equal
-
-from commpy.sequences import pnsequence
+from numpy.testing import run_module_suite, assert_raises, assert_equal, assert_almost_equal
 
+from commpy.sequences import pnsequence, zcsequence
 
 def test_pnsequence():
     # Test the raises of errors
     with assert_raises(ValueError):
         pnsequence(4, '001', '1101', 2**4 - 1)
+    with assert_raises(ValueError):
         pnsequence(4, '0011', '110', 2 ** 4 - 1)
 
     # Test output with
@@ -19,6 +20,34 @@ def test_pnsequence():
     assert_equal(pnsequence(4, (0, 0, 1, 1), array((1, 1, 0, 1)), 7), array((1, 1, 0, 0, 1, 0, 1), int),
                  err_msg='Pseudo-noise sequence is not the one expected.')
 
+def test_zcsequence():
+    # Test the raises of errors
+    with assert_raises(ValueError):
+        zcsequence(u=-1, seq_length=20, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=20, seq_length=0, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=3, seq_length=18, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=3.1, seq_length=11, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=3, seq_length=11.1, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=3, seq_length=11, q=0.1)
+
+    # Test output with
+    assert_almost_equal(zcsequence(u=1, seq_length=2, q=0), array([1.000000e+00+0.j, 6.123234e-17-1.j]),
+        err_msg='CAZAC sequence is not the expected one.')
+
+    # Test if output cross-correlation is valid
+    seqCAZAC = zcsequence(u=3, seq_length=20, q=0)
+    x = np.fft.fft(seqCAZAC) / np.sqrt(seqCAZAC.size)
+    h = (np.fft.ifft(np.conj(x) * x)*np.sqrt(seqCAZAC.size)).T
+    corr = np.absolute(h)**2/h.size
+    assert_almost_equal(corr[0], 1.,
+        err_msg='CAZAC sequence auto-correlation is not valid, first term is not 1')
+    assert_almost_equal(corr[1:], np.zeros(corr.size-1),
+        err_msg='CAZAC sequence auto-correlation is not valid, all terms except first are not 0')
 
 if __name__ == "__main__":
     run_module_suite()