Remove slow parallels (#39)

* init

* tests

* versionup

extra/setup.py

@@ -33,7 +33,9 @@
     extraCompileArgs = ['/DNPY_NO_DEPRECATED_API', '/O2', '/Oiy', '/Ob3', '/favor:INTEL64', '/options:strict',# '/openmp',
                         '/fp:fast', '/fp:except-', '/GL', '/Gw', '/jumptablerdata', '/MP', '/Qpar']
 
-setup(name="dsp",
+setup(name="sdrterm",
+    version="0.4.1",
+    python_requires=">=3.10",
     ext_modules=cythonize([
     Extension('dsp.fast.iq_correction',
               sources=['src/iq_correction.pyx'],

extra/src/iq_correction.pyx

@@ -37,24 +37,19 @@ from numpy cimport ndarray
 cdef class IQCorrection:
     cdef unsigned long _fs
     cdef readonly double inductance
-    cdef double complex _off
 
     def __cinit__(self, const unsigned long fs, const unsigned short impedance = 50):
         self._fs = fs
         self.inductance = impedance / fs
-        self._off = 0j
 
     @cython.cdivision(True)
-    cpdef void correctIq(self, ndarray[np.complex128_t] data):
+    cpdef void correctIq(self, ndarray[np.complex128_t] data, ndarray[np.complex128_t] off):
         cdef Py_ssize_t i
-        cdef Py_ssize_t size = data.size
-        # for i in prange(size, nogil=True):
+        cdef Py_ssize_t size = data.shape[0]
         with nogil:
             for i in range(size):
-                # *SIGH* VS bitches about the arithmetic assignment
-                # operator, bc ofc it does
-                data[i] = data[i] - self._off
-                self._off += data[i] * self.inductance
+                data[i] = data[i] - off[0]
+                off[0] += data[i] * self.inductance
 
     @property
     def fs(self):
@@ -73,4 +68,3 @@ cdef class IQCorrection:
     @impedance.setter
     def impedance(self, const unsigned short impedance):
         self.inductance = impedance / self._fs
-        self._off = 0j

pyproject.toml

@@ -8,7 +8,7 @@ sources = ["src"]
 
 [project]
 name = "sdrterm"
-version = "0.4"
+version = "0.4.1"
 dependencies = [
     'scipy',
     'numpy>=1.26',

src/dsp/demodulation.py

@@ -17,50 +17,62 @@
 # You should have received a copy of the GNU General Public License
 # along with this program. If not, see <http://www.gnu.org/licenses/>.
 #
-from numba import njit, guvectorize, complex128, float64
-from numpy import angle, ndarray, conj, abs, real, imag, dtype, complexfloating, floating, reshape
+from numba import guvectorize, complex128, float64
+from numpy import angle, ndarray, conj, abs, real, imag, dtype, complexfloating, floating, square
 from scipy.signal import resample
 
 
 @guvectorize([(complex128[:], float64[:])], '(n)->(n)',
              nopython=True,
              cache=True,
              boundscheck=False,
-             target='parallel')
+             fastmath=True)
 def _fmDemod(data: ndarray[any, dtype[complexfloating]], res: ndarray[any, dtype[floating]]):
     for i in range(0, data.shape[0], 2):
         res[i >> 1] = angle(data[i] * conj(data[i + 1]))
 
 
 def fmDemod(data: ndarray[any, dtype[complexfloating]], tmp: ndarray[any, dtype[floating]]):
-    # tmp1 = resample(angle(data[::2] * conj(data[1::2])), data.size)
-    if data.ndim < 2:
-        data[:] = reshape(data, (1, data.shape[0]))
     _fmDemod(data, tmp)
     for i in range(tmp.shape[0]):
         tmp[i] = resample(tmp[i][:tmp.shape[1] >> 1], tmp.shape[1])
 
 
-@njit(cache=True, nogil=True, error_model='numpy', boundscheck=False)
+@guvectorize([(complex128[:], float64[:])], '(n)->(n)',
+             nopython=True,
+             cache=True,
+             boundscheck=False,
+             fastmath=True)
 def amDemod(data: ndarray[any, dtype[complexfloating]], res: ndarray[any, dtype[floating]]):
-    res[:] = abs(data)
+    for i in range(data.shape[0]):
+        res[i] = abs(square(data[i]))
 
 
-@njit(cache=True, nogil=True, error_model='numpy', boundscheck=False)
+@guvectorize([(complex128[:], float64[:])], '(n)->(n)',
+             nopython=True,
+             cache=True,
+             boundscheck=False,
+             fastmath=True)
 def realOutput(data: ndarray[any, dtype[complexfloating]], res: ndarray[any, dtype[floating]]):
-    res[:] = real(data)
+    for i in range(data.shape[0]):
+        res[i] = real(data[i])
 
 
-@njit(cache=True, nogil=True, error_model='numpy', boundscheck=False)
+@guvectorize([(complex128[:], float64[:])], '(n)->(n)',
+             nopython=True,
+             cache=True,
+             boundscheck=False,
+             fastmath=True)
 def imagOutput(data: ndarray[any, dtype[complexfloating]], res: ndarray[any, dtype[floating]]):
-    res[:] = imag(data)
+    for i in range(data.shape[0]):
+        res[i] = imag(data[i])
 
 
 @guvectorize([(complex128[:], complex128[:, :], complex128[:, :])], '(n),(m,n)->(m,n)',
              nopython=True,
              cache=True,
              boundscheck=False,
-             target='parallel')
+             fastmath=True)
 def shiftFreq(y: ndarray[any, dtype[complexfloating]],
               shift: ndarray[any, dtype[complexfloating]],
               res: ndarray[any, dtype[complexfloating]]) -> None:

src/dsp/dsp_processor.py

@@ -21,8 +21,7 @@
 from sys import stdout
 from typing import Callable, Iterable, Any
 
-from numpy import ndarray, dtype, complexfloating, floating, exp, arange, pi, empty, array, reshape, \
-    broadcast_to
+from numpy import ndarray, dtype, complexfloating, floating, exp, arange, pi, empty, array
 from scipy.signal import decimate, dlti, savgol_filter, sosfilt, ellip
 
 from dsp.data_processor import DataProcessor
@@ -141,7 +140,7 @@ def selectOutputImag(self):
     def _processChunk(self,
                       x: ndarray[any, dtype[complexfloating]],
                       y: ndarray[any, dtype[complexfloating]],
-                      z: ndarray[any, dtype[floating]])->None:
+                      z: ndarray[any, dtype[floating]]) -> None:
         if self._shift is not None:
             shiftFreq(x[0], self._shift, x)
             # y = y * self._shift

src/misc/read_file.py

@@ -22,8 +22,8 @@
 from sys import stdin
 from typing import Iterable
 
-from numba import njit
-from numpy import frombuffer, ndarray, complexfloating, dtype, empty, uint8
+from numba import guvectorize, complex128 as nbComplex128
+from numpy import frombuffer, ndarray, complexfloating, dtype, empty, uint8, complex128, array
 
 from misc.general_util import vprint, eprint, tprint, applyIgnoreException
 
@@ -40,6 +40,7 @@ def readFile(bitsPerSample: dtype = None,
              correctIq: bool = False,
              normalize: bool = False,
              isSocket: bool = False,
+             impedance: int = 50,
              **_) -> None:
     if fs is None:
         raise ValueError('fs is not specified')
@@ -49,42 +50,56 @@ def readFile(bitsPerSample: dtype = None,
 
     dataType = dtype([('re', bitsPerSample), ('im', bitsPerSample)])
     buffer = empty(readSize, dtype=uint8)
+    offset = array([complex128(0j)])
 
-    def _correctIq(_: ndarray[any, dtype[complexfloating]]) -> None:
+    def _correctIq(*_) -> None:
         pass
 
     if correctIq:
         try:
             from dsp.fast.iq_correction import IQCorrection
             tprint('Imported pre-compiled IQCorrection class')
+            _correctIq = IQCorrection(fs).correctIq
         except ImportError:
-            from dsp.iq_correction import IQCorrection
-            tprint('Falling-back to jit IQCorrection class')
-        _correctIq = IQCorrection(fs).correctIq
+            tprint('Falling back to local IQCorrection')
+            inductance: float = impedance / fs
+
+            @guvectorize([(nbComplex128[:], nbComplex128[:])], '(n)->()',
+                         nopython=True,
+                         cache=True,
+                         boundscheck=False,
+                         fastmath=True)
+            def _correctIq(z, res) -> None:
+                off = res[0]
+                for i in range(z.shape[0]):
+                    z[i] = z[i] - off
+                    off += z[i] * inductance
+                res[0] = off
 
-    def noNormalize(*_) -> None:
+    def _normalize(*_) -> None:
         pass
 
-    if not normalize:
-        _normalize = noNormalize
-    else:
+    if normalize:
         ret = generateDomain(bitsPerSample.char)
-        if ret is None:
-            _normalize = noNormalize
-        else:
+        if ret is not None:
             xmin, xMaxMinDiff = ret
 
-            @njit(cache=True, nogil=True, error_model='numpy', boundscheck=False)
-            def _normalize(z: ndarray[any, dtype[complexfloating]]) -> None:
+            @guvectorize([(nbComplex128[:], nbComplex128[:])], '(n)->(n)',
+                         nopython=True,
+                         cache=True,
+                         boundscheck=False,
+                         fastmath=True)
+            def _normalize(z: ndarray[any, dtype[complexfloating]],
+                           res: ndarray[any, dtype[complexfloating]]) -> None:
                 for i in range(z.shape[0]):
-                    z[i] = 1.6 * (z[i] - xmin) * xMaxMinDiff - 0.8
+                    res[i] = 1.6 * (z[i] - xmin) * xMaxMinDiff - 0.8
     procs = list(processes)
     clients = list(buffers)
 
     def feedBuffers(y: ndarray) -> None:
         z = y['re'] + 1j * y['im']
-        _normalize(z)
-        _correctIq(z)
+        _normalize(z, z)
+        _correctIq(z, offset)
         for proc, client in zip(procs, clients):
             if proc.exitcode is not None:
                 tprint(f'Process : {proc.name} ended; removing {client} from queue')

src/sdrterm.py

@@ -172,6 +172,12 @@ def main(fs: Annotated[int, Option('--fs', '-r',
 def __setStartMethod():
     from multiprocessing import get_all_start_methods
     from misc.general_util import printException
+    from os import environ
+    from sys import warnoptions
+    from warnings import simplefilter
+    if not warnoptions:
+        simplefilter('ignore')
+        environ["PYTHONWARNINGS"] = 'ignore'
 
     if 'spawn' in get_all_start_methods():
         from multiprocessing import set_start_method

test/dsp/demodulation_test.py

@@ -31,7 +31,7 @@ def test_fm(data):
 def test_am(data):
     inp, outp = data
     dsp.amDemod(np.array(inp), outp)
-    testOutp = [math.sqrt(math.pow(z.real, 2) + math.pow(z.imag, 2)) for z in inp[0]]
+    testOutp = [math.pow(z.real, 2) + math.pow(z.imag, 2) for z in inp[0]]
     for x, y in zip(outp[0], testOutp):
         assert math.fabs(x - y) < EPSILON