use costas loop for PSK demod (#807)

Author: jopohl

data/azure-pipelines.yml

@@ -19,6 +19,8 @@ jobs:
         python.version: '3.7'
       Python38:
           python.version: '3.8'
+      Python39:
+          python.version: '3.9'
 
   steps:
   - task: UsePythonVersion@0
@@ -114,7 +116,7 @@ jobs:
 
 - job: 'Windows'
   pool:
-    vmImage: 'vs2017-win2016'
+    vmImage: 'windows-2019'
   strategy:
     matrix:
       Python37-64:
@@ -240,16 +242,12 @@ jobs:
     displayName: "download and unpack SDR drivers"
 
   - script: |
-      python -m pip install --upgrade pip
+      python -m pip install --upgrade pip wheel setuptools
       python -m pip install --upgrade -r data/requirements.txt
-      python -m pip install --upgrade pytest pytest-faulthandler
-    displayName: 'Install dependencies'
-
-  - script: |
-      brew install airspy hackrf librtlsdr portaudio uhd
-      python -m pip install --upgrade wheel twine six appdirs packaging setuptools pyinstaller pyaudio
+      HOMEBREW_NO_INSTALL_CLEANUP=TRUE brew install airspy hackrf librtlsdr portaudio uhd
+      python -m pip install --upgrade pytest pytest-faulthandler twine six appdirs packaging pyinstaller pyaudio
       python -c "import tempfile, os; open(os.path.join(tempfile.gettempdir(), 'urh_releasing'), 'w').close()"
-    displayName: "Install build dependencies"
+    displayName: 'Install dependencies'
 
   - script: python src/urh/cythonext/build.py
     displayName: "Build extensions"

src/urh/ainterpretation/AutoInterpretation.py

@@ -382,11 +382,11 @@ def estimate(iq_array: IQArray, noise: float = None, modulation: str = None) ->
         message_indices = merge_message_segments_for_ook(message_indices)
 
     if modulation == "OOK" or modulation == "ASK":
-        data = signal_functions.afp_demod(iq_array.data, noise, "ASK")
+        data = signal_functions.afp_demod(iq_array.data, noise, "ASK", 2)
     elif modulation == "FSK":
-        data = signal_functions.afp_demod(iq_array.data, noise, "FSK")
+        data = signal_functions.afp_demod(iq_array.data, noise, "FSK", 2)
     elif modulation == "PSK":
-        data = signal_functions.afp_demod(iq_array.data, noise, "PSK")
+        data = signal_functions.afp_demod(iq_array.data, noise, "PSK", 2)
     else:
         raise ValueError("Unsupported Modulation")
 

src/urh/controller/dialogs/CostaOptionsDialog.py

@@ -0,0 +1,27 @@
+from PyQt5.QtCore import Qt, pyqtSlot
+from PyQt5.QtWidgets import QDialog
+
+from urh.ui.ui_costa import Ui_DialogCosta
+
+
+class CostaOptionsDialog(QDialog):
+    def __init__(self, loop_bandwidth, parent=None):
+        super().__init__(parent)
+        self.ui = Ui_DialogCosta()
+        self.ui.setupUi(self)
+        self.setAttribute(Qt.WA_DeleteOnClose)
+        self.setWindowFlags(Qt.Window)
+
+        self.costas_loop_bandwidth = loop_bandwidth
+        self.ui.doubleSpinBoxLoopBandwidth.setValue(self.costas_loop_bandwidth)
+
+        self.create_connects()
+
+    def create_connects(self):
+        self.ui.buttonBox.accepted.connect(self.accept)
+        self.ui.buttonBox.rejected.connect(self.reject)
+        self.ui.doubleSpinBoxLoopBandwidth.valueChanged.connect(self.on_spinbox_loop_bandwidth_value_changed)
+
+    @pyqtSlot(float)
+    def on_spinbox_loop_bandwidth_value_changed(self, value):
+        self.costas_loop_bandwidth = value

src/urh/controller/widgets/SignalFrame.py

@@ -10,6 +10,7 @@
 
 from urh import settings
 from urh.controller.dialogs.AdvancedModulationOptionsDialog import AdvancedModulationOptionsDialog
+from urh.controller.dialogs.CostaOptionsDialog import CostaOptionsDialog
 from urh.controller.dialogs.FilterDialog import FilterDialog
 from urh.controller.dialogs.SendDialog import SendDialog
 from urh.controller.dialogs.SignalDetailsDialog import SignalDetailsDialog
@@ -276,7 +277,7 @@ def refresh_signal_information(self, block=True):
         self.ui.spinBoxSamplesPerSymbol.setValue(self.signal.samples_per_symbol)
         self.ui.spinBoxNoiseTreshold.setValue(self.signal.noise_threshold_relative)
         self.ui.cbModulationType.setCurrentText(self.signal.modulation_type)
-        self.ui.btnAdvancedModulationSettings.setVisible(self.ui.cbModulationType.currentText() == "ASK")
+        self.ui.btnAdvancedModulationSettings.setVisible(self.ui.cbModulationType.currentText() in ("ASK", "PSK"))
         self.ui.spinBoxCenterSpacing.setValue(self.signal.center_spacing)
         self.ui.spinBoxBitsPerSymbol.setValue(self.signal.bits_per_symbol)
 
@@ -524,7 +525,10 @@ def update_protocol(self):
         self.ui.txtEdProto.setText("Demodulating...")
         qApp.processEvents()
 
-        self.proto_analyzer.get_protocol_from_signal()
+        try:
+            self.proto_analyzer.get_protocol_from_signal()
+        except Exception as e:
+            Errors.exception(e)
 
     def show_protocol(self, old_view=-1, refresh=False):
         if not self.proto_analyzer:
@@ -1093,7 +1097,7 @@ def on_combobox_modulation_type_text_changed(self, txt: str):
                 self.scene_manager.init_scene()
                 self.on_slider_y_scale_value_changed()
 
-        self.ui.btnAdvancedModulationSettings.setVisible(self.ui.cbModulationType.currentText() == "ASK")
+        self.ui.btnAdvancedModulationSettings.setVisible(self.ui.cbModulationType.currentText() in ("ASK", "PSK"))
 
     @pyqtSlot()
     def on_signal_data_changed_before_save(self):
@@ -1301,9 +1305,25 @@ def get_advanced_modulation_settings_dialog(self):
         dialog.message_length_divisor_edited.connect(self.on_message_length_divisor_edited)
         return dialog
 
+    def get_costas_dialog(self):
+        dialog = CostaOptionsDialog(self.signal.costas_loop_bandwidth, parent=self)
+        dialog.accepted.connect(self.on_costas_dialog_accepted)
+        return dialog
+
+    @pyqtSlot()
+    def on_costas_dialog_accepted(self):
+        sender = self.sender()
+        assert isinstance(sender, CostaOptionsDialog)
+        self.signal.costas_loop_bandwidth = sender.costas_loop_bandwidth
+
     @pyqtSlot()
     def on_btn_advanced_modulation_settings_clicked(self):
-        dialog = self.get_advanced_modulation_settings_dialog()
+        if self.ui.cbModulationType.currentText() == "ASK":
+            dialog = self.get_advanced_modulation_settings_dialog()
+        elif self.ui.cbModulationType.currentText() == "PSK":
+            dialog = self.get_costas_dialog()
+        else:
+            raise ValueError("No additional settings available")
         dialog.exec_()
 
     @pyqtSlot()

src/urh/cythonext/build.py

@@ -10,8 +10,9 @@ def main():
     cur_dir = os.path.realpath(__file__)
     os.chdir(os.path.realpath(os.path.join(cur_dir, "..", "..", "..", "..")))
     # call([sys.executable, "setup.py", "clean", "--all"])
-    call([sys.executable, "setup.py", "build_ext", "--inplace", "-j{}".format(os.cpu_count())])
+    rc = call([sys.executable, "setup.py", "build_ext", "--inplace", "-j{}".format(os.cpu_count())])
+    return rc
 
 
 if __name__ == "__main__":
-    main()
+    sys.exit(main())

src/urh/cythonext/signal_functions.pyx

@@ -242,17 +242,27 @@ cdef np.ndarray[np.float32_t, ndim=1] gauss_fir(float sample_rate, uint32_t samp
         -(((np.sqrt(2) * np.pi) / np.sqrt(np.log(2)) * bt * k / samples_per_symbol) ** 2))
     return h / h.sum()
 
-cdef void phase_demod(IQ samples, float[::1] result, float noise_sqrd, bool qam, long long num_samples):
-    cdef long long i = 0
-    cdef float real = 0, imag = 0, magnitude = 0
+cdef float clamp(float x) nogil:
+    if x < -1.0:
+        x = -1.0
+    elif x > 1.0:
+        x = 1.0
+    return x
+
+cdef float[::1] costa_demod(IQ samples, float noise_sqrd, int loop_order, float bandwidth=0.1, float damping=sqrt(2.0) / 2.0):
+    cdef float alpha = (4 * damping * bandwidth) / (1.0 + 2.0 * damping * bandwidth + bandwidth * bandwidth)
+    cdef float beta = (4 * bandwidth * bandwidth) / (1.0 + 2.0 * damping * bandwidth + bandwidth * bandwidth)
+
+    cdef long long i = 0, num_samples = len(samples)
+    cdef float real = 0, imag = 0
 
     cdef float scale, shift, real_float, imag_float, ref_real, ref_imag
 
-    cdef float phi = 0, current_arg = 0, f_curr = 0, f_prev = 0
+    cdef float f1, f2, costa_freq = 0, costa_error = 0, costa_phase = 1.5
 
-    cdef float complex current_sample, conj_previous_sample, current_nco
+    cdef float complex current_sample, nco_out, nco_times_sample
 
-    cdef float alpha = 0.1
+    cdef float[::1] result = np.empty(num_samples, dtype=np.float32)
 
     if str(cython.typeof(samples)) == "char[:, ::1]":
         scale = 127.5
@@ -272,62 +282,64 @@ cdef void phase_demod(IQ samples, float[::1] result, float noise_sqrd, bool qam,
     else:
         raise ValueError("Unsupported dtype")
 
+    if loop_order > 4:
+        # TODO: Adapt this when PSK demodulation with order > 4 shall be supported
+        loop_order = 4
+
     for i in range(1, num_samples):
         real = samples[i, 0]
         imag = samples[i, 1]
-        
-        magnitude = real * real + imag * imag
-        if magnitude <= noise_sqrd:
+
+        if real * real + imag * imag <= noise_sqrd:
             result[i] = NOISE_FSK_PSK
             continue
 
         real_float = (real + shift) / scale
         imag_float = (imag + shift) / scale
 
         current_sample = real_float + imag_unit * imag_float
-        conj_previous_sample = (samples[i-1, 0] + shift) / scale - imag_unit * ((samples[i-1, 1] + shift) / scale)
-        f_curr = arg(current_sample * conj_previous_sample)
+        nco_out = cosf(-costa_phase) + imag_unit * sinf(-costa_phase)
+        nco_times_sample = nco_out * current_sample
 
-        if abs(f_curr) < M_PI / 4:  # TODO: For PSK with order > 4 this needs to be adapted
-            f_prev = f_curr
-            current_arg += f_curr
-        else:
-            current_arg += f_prev
+        if loop_order == 2:
+            costa_error = nco_times_sample.imag * nco_times_sample.real
+        elif loop_order == 4:
+            f1 = 1.0 if nco_times_sample.real > 0.0 else -1.0
+            f2 = 1.0 if nco_times_sample.imag > 0.0 else -1.0
+            costa_error = f1 * nco_times_sample.imag - f2 * nco_times_sample.real
 
-        # Reference oscillator cos(current_arg) + j * sin(current_arg)
-        current_nco = cosf(current_arg) + imag_unit * sinf(current_arg)
-        phi = arg(current_sample * conj(current_nco))
+        costa_error = clamp(costa_error)
+
+        # advance the loop
+        costa_freq += beta * costa_error
+        costa_phase += costa_freq + alpha * costa_error
+
+        # wrap the phase
+        while costa_phase > (2 * M_PI):
+            costa_phase -= 2 * M_PI
+        while costa_phase < (-2 * M_PI):
+            costa_phase += 2 * M_PI
+
+        costa_freq = clamp(costa_freq)
+
+        if loop_order == 2:
+            result[i] = nco_times_sample.real
+        elif loop_order == 4:
+            result[i] = 2 * nco_times_sample.real + nco_times_sample.imag
+
+    return result
 
-        if qam:
-            result[i] = phi * magnitude
-        else:
-            result[i] = phi
 
-cpdef np.ndarray[np.float32_t, ndim=1] afp_demod(IQ samples, float noise_mag, str mod_type):
+cpdef np.ndarray[np.float32_t, ndim=1] afp_demod(IQ samples, float noise_mag,
+                                                 str mod_type, int mod_order, float costas_loop_bandwidth=0.1):
     if len(samples) <= 2:
         return np.zeros(len(samples), dtype=np.float32)
 
     cdef long long i = 0, ns = len(samples)
-    cdef float current_arg = 0
-    cdef float noise_sqrd = 0
-    cdef float complex_phase = 0
-    cdef float prev_phase = 0
-    cdef float NOISE = 0
-    cdef float real = 0
-    cdef float imag = 0
-
-    cdef float[::1] result = np.zeros(ns, dtype=np.float32, order="C")
-    cdef float costa_freq = 0
-    cdef float costa_phase = 0
-    cdef complex nco_out = 0
+    cdef float NOISE = get_noise_for_mod_type(mod_type)
+    cdef float noise_sqrd = noise_mag * noise_mag, real = 0, imag = 0, magnitude = 0, max_magnitude
     cdef float complex tmp
-    cdef float phase_error = 0
-    cdef float costa_alpha = 0
-    cdef float costa_beta = 0
-    cdef complex nco_times_sample = 0
-    cdef float magnitude = 0
 
-    cdef float max_magnitude   # ensure all magnitudes of ASK demod between 0 and 1
     if str(cython.typeof(samples)) == "char[:, ::1]":
         max_magnitude = sqrt(127*127 + 128*128)
     elif str(cython.typeof(samples)) == "unsigned char[:, ::1]":
@@ -341,35 +353,27 @@ cpdef np.ndarray[np.float32_t, ndim=1] afp_demod(IQ samples, float noise_mag, st
     else:
         raise ValueError("Unsupported dtype")
 
-    # Atan2 yields values from -Pi to Pi
-    # We use the Magic Constant NOISE_FSK_PSK to cut off noise
-    noise_sqrd = noise_mag * noise_mag
-    NOISE = get_noise_for_mod_type(mod_type)
-    result[0] = NOISE
-
-    cdef bool qam = False
 
-    if mod_type in ("PSK", "QAM", "OQPSK"):
-        if mod_type == "QAM":
-            qam = True
+    if mod_type == "PSK":
+        return np.asarray(costa_demod(samples, noise_sqrd, mod_order, bandwidth=costas_loop_bandwidth))
 
-        phase_demod(samples, result, noise_sqrd, qam, ns)
+    cdef float[::1] result = np.zeros(ns, dtype=np.float32, order="C")
+    result[0] = NOISE
 
-    else:
-        for i in prange(1, ns, nogil=True, schedule="static"):
-            real = samples[i, 0]
-            imag = samples[i, 1]
-            magnitude = real * real + imag * imag
-            if magnitude <= noise_sqrd:  # |c| <= mag_treshold
-                result[i] = NOISE
-                continue
+    for i in prange(1, ns, nogil=True, schedule="static"):
+        real = samples[i, 0]
+        imag = samples[i, 1]
+        magnitude = real * real + imag * imag
+        if magnitude <= noise_sqrd:  # |c| <= mag_treshold
+            result[i] = NOISE
+            continue
 
-            if mod_type == "ASK":
-                result[i] = sqrt(magnitude) / max_magnitude
-            elif mod_type == "FSK":
-                #tmp = samples[i - 1].conjugate() * c
-                tmp = (samples[i-1, 0] - imag_unit * samples[i-1, 1]) * (real + imag_unit * imag)
-                result[i] = atan2(tmp.imag, tmp.real)  # Freq
+        if mod_type == "ASK":
+            result[i] = sqrt(magnitude) / max_magnitude
+        elif mod_type == "FSK":
+            #tmp = samples[i - 1].conjugate() * c
+            tmp = (samples[i-1, 0] - imag_unit * samples[i-1, 1]) * (real + imag_unit * imag)
+            result[i] = atan2(tmp.imag, tmp.real)  # Freq
 
     return np.asarray(result)
 

src/urh/signalprocessing/Signal.py

@@ -43,6 +43,7 @@ def __init__(self, filename: str, name="Signal", modulation: str = None, sample_
         self.__samples_per_symbol = 100
         self.__pause_threshold = 8
         self.__message_length_divisor = 1
+        self.__costas_loop_bandwidth = 0.1
         self._qad = None
         self.__center = 0
         self._noise_threshold = 0
@@ -255,6 +256,18 @@ def pause_threshold(self, value: int):
             if not self.block_protocol_update:
                 self.protocol_needs_update.emit()
 
+    @property
+    def costas_loop_bandwidth(self):
+        return self.__costas_loop_bandwidth
+
+    @costas_loop_bandwidth.setter
+    def costas_loop_bandwidth(self, value: float):
+        if self.__costas_loop_bandwidth != value:
+            self.__costas_loop_bandwidth = value
+            self._qad = None
+            if not self.block_protocol_update:
+                self.protocol_needs_update.emit()
+
     @property
     def message_length_divisor(self) -> int:
         return self.__message_length_divisor
@@ -362,7 +375,9 @@ def save_as(self, filename: str):
         QApplication.instance().restoreOverrideCursor()
 
     def quad_demod(self):
-        return signal_functions.afp_demod(self.iq_array.data, self.noise_threshold, self.modulation_type)
+        return signal_functions.afp_demod(self.iq_array.data, self.noise_threshold,
+                                          self.modulation_type, self.modulation_order,
+                                          self.costas_loop_bandwidth)
 
     def calc_relative_noise_threshold_from_range(self, noise_start: int, noise_end: int):
         num_digits = 4
@@ -501,7 +516,10 @@ def crop_to_range(self, start: int, end: int):
     def filter_range(self, start: int, end: int, fir_filter: Filter):
         self.iq_array[start:end] = fir_filter.work(self.iq_array[start:end])
         self._qad[start:end] = signal_functions.afp_demod(self.iq_array[start:end],
-                                                          self.noise_threshold, self.modulation_type)
+                                                          self.noise_threshold,
+                                                          self.modulation_type,
+                                                          self.modulation_order,
+                                                          self.costas_loop_bandwidth)
         self.__invalidate_after_edit()
 
     def __invalidate_after_edit(self):