diff options
Diffstat (limited to '_signal_processing')
| -rwxr-xr-x | _signal_processing/krakenSDR_signal_processor.py | 713 |
1 files changed, 713 insertions, 0 deletions
diff --git a/_signal_processing/krakenSDR_signal_processor.py b/_signal_processing/krakenSDR_signal_processor.py new file mode 100755 index 0000000..31f1836 --- /dev/null +++ b/_signal_processing/krakenSDR_signal_processor.py @@ -0,0 +1,713 @@ +# KrakenSDR Signal Processor +# +# Copyright (C) 2018-2021 Carl Laufer, Tamás Pető +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <https://www.gnu.org/licenses/>. +# +# +# - coding: utf-8 -*- + +# Import built-in modules +import sys +import os +import time +import logging +import threading +import queue +import math + +# Import optimization modules +import numba as nb +from numba import jit, njit +from functools import lru_cache + +# Math support +import numpy as np +import numpy.linalg as lin +#from numba import jit +import pyfftw + +# Signal processing support +import scipy +from scipy import fft +from scipy import signal +from scipy.signal import correlate +from scipy.signal import convolve + +from pyapril import channelPreparation as cp +from pyapril import clutterCancellation as cc +from pyapril import detector as det + +c_dtype = np.complex64 + +#import socket +# UDP is useless to us because it cannot work over mobile internet + +# Init UDP +#server = socket.socket(socket.AF_INET, socket.SOCK_DGRAM, socket.IPPROTO_UDP) +#server.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEPORT, 1) +# Enable broadcasting mode +#server.setsockopt(socket.SOL_SOCKET, socket.SO_BROADCAST, 1) +# Set a timeout so the socket does not block +# indefinitely when trying to receive data. +#server.settimeout(0.2) + +class SignalProcessor(threading.Thread): + + def __init__(self, data_que, module_receiver, logging_level=10): + + """ + Parameters: + ----------- + :param: data_que: Que to communicate with the UI (web iface/Qt GUI) + :param: module_receiver: Kraken SDR DoA DSP receiver modules + """ + super(SignalProcessor, self).__init__() + self.logger = logging.getLogger(__name__) + self.logger.setLevel(logging_level) + + root_path = os.path.dirname(os.path.dirname(os.path.realpath(__file__))) + doa_res_file_path = os.path.join(os.path.join(root_path,"_android_web","DOA_value.html")) + self.DOA_res_fd = open(doa_res_file_path,"w+") + + self.module_receiver = module_receiver + self.data_que = data_que + self.en_spectrum = False + self.en_record = False + self.en_DOA_estimation = True + self.first_frame = 1 # Used to configure local variables from the header fields + self.processed_signal = np.empty(0) + + # Squelch feature + self.data_ready = False + self.en_squelch = False + self.squelch_threshold = 0.1 + self.squelch_trigger_channel = 0 + self.raw_signal_amplitude = np.empty(0) + self.filt_signal = np.empty(0) + self.squelch_mask = np.empty(0) + + # DOA processing options + self.en_DOA_Bartlett = False + self.en_DOA_Capon = False + self.en_DOA_MEM = False + self.en_DOA_MUSIC = False + self.en_DOA_FB_avg = False + self.DOA_offset = 0 + self.DOA_inter_elem_space = 0.5 + self.DOA_ant_alignment = "ULA" + self.DOA_theta = np.linspace(0,359,360) + + # PR processing options + self.PR_clutter_cancellation = "Wiener MRE" + self.max_bistatic_range = 128 + self.max_doppler = 256 + self.en_PR = True + + + # Processing parameters + self.spectrum_window_size = 2048 #1024 + self.spectrum_window = "hann" + self.run_processing = False + self.is_running = False + + + self.channel_number = 4 # Update from header + + # Result vectors + self.DOA_Bartlett_res = np.ones(181) + self.DOA_Capon_res = np.ones(181) + self.DOA_MEM_res = np.ones(181) + self.DOA_MUSIC_res = np.ones(181) + self.DOA_theta = np.arange(0,181,1) + + self.max_index = 0 + self.max_frequency = 0 + self.fft_signal_width = 0 + + self.DOA_theta = np.linspace(0,359,360) + + self.spectrum = None #np.ones((self.channel_number+2,N), dtype=np.float32) + self.spectrum_upd_counter = 0 + + + def run(self): + """ + Main processing thread + """ + + pyfftw.config.NUM_THREADS = 4 + scipy.fft.set_backend(pyfftw.interfaces.scipy_fft) + pyfftw.interfaces.cache.enable() + + while True: + self.is_running = False + time.sleep(1) + while self.run_processing: + self.is_running = True + + que_data_packet = [] + start_time = time.time() + + #-----> ACQUIRE NEW DATA FRAME <----- + self.module_receiver.get_iq_online() + + # Check frame type for processing + en_proc = (self.module_receiver.iq_header.frame_type == self.module_receiver.iq_header.FRAME_TYPE_DATA)# or \ + #(self.module_receiver.iq_header.frame_type == self.module_receiver.iq_header.FRAME_TYPE_CAL)# For debug purposes + """ + You can enable here to process other frame types (such as call type frames) + """ + + que_data_packet.append(['iq_header',self.module_receiver.iq_header]) + self.logger.debug("IQ header has been put into the data que entity") + + # Configure processing parameteres based on the settings of the DAQ chain + if self.first_frame: + self.channel_number = self.module_receiver.iq_header.active_ant_chs + self.spectrum_upd_counter = 0 + self.spectrum = np.ones((self.channel_number+2, self.spectrum_window_size), dtype=np.float32) + self.first_frame = 0 + + decimation_factor = 1 + + self.data_ready = False + + if en_proc: + self.processed_signal = self.module_receiver.iq_samples + self.data_ready = True + + first_decimation_factor = 1 #480 + + # TESTING: DSP side main decimation - significantly slower than NE10 but it works ok-ish + #decimated_signal = signal.decimate(self.processed_signal, first_decimation_factor, n = 584, ftype='fir', zero_phase=True) #first_decimation_factor * 2, ftype='fir') + #self.processed_signal = decimated_signal #.copy() + #spectrum_signal = decimated_signal.copy() + + max_amplitude = -100 + + #max_ch = np.argmax(np.max(self.spectrum[1:self.module_receiver.iq_header.active_ant_chs+1,:], axis=1)) # Find the channel that had the max amplitude + max_amplitude = 0 #np.max(self.spectrum[1+max_ch, :]) #Max amplitude out of all 5 channels + #max_spectrum = self.spectrum[1+max_ch, :] #Send max ch to channel centering + + que_data_packet.append(['max_amplitude',max_amplitude]) + + #-----> SQUELCH PROCESSING <----- + + if self.en_squelch: + self.data_ready = False + + self.processed_signal, decimation_factor, self.fft_signal_width, self.max_index = \ + center_max_signal(self.processed_signal, self.spectrum[0,:], max_spectrum, self.module_receiver.daq_squelch_th_dB, self.module_receiver.iq_header.sampling_freq) + + #decimated_signal = [] + #if(decimation_factor > 1): + # decimated_signal = signal.decimate(self.processed_signal, decimation_factor, n = decimation_factor * 2, ftype='fir') + # self.processed_signal = decimated_signal #.copy() + + + #Only update if we're above the threshold + if max_amplitude > self.module_receiver.daq_squelch_th_dB: + self.data_ready = True + + + #-----> SPECTRUM PROCESSING <----- + + if self.en_spectrum and self.data_ready: + + spectrum_samples = self.module_receiver.iq_samples #spectrum_signal #self.processed_signal #self.module_receiver.iq_samples #self.processed_signal + + + N = self.spectrum_window_size + + N_perseg = 0 + N_perseg = min(N, len(self.processed_signal[0,:])//25) + N_perseg = N_perseg // 1 + + # Get power spectrum + f, Pxx_den = signal.welch(self.processed_signal, self.module_receiver.iq_header.sampling_freq//first_decimation_factor, + nperseg=N_perseg, + nfft=N, + noverlap=int(N_perseg*0.25), + detrend=False, + return_onesided=False, + window= ('tukey', 0.25), #tukey window gives better time resolution for squelching #self.spectrum_window, #('tukey', 0.25), #self.spectrum_window, + #window=self.spectrum_window, + scaling="spectrum") + + self.spectrum[1:self.module_receiver.iq_header.active_ant_chs+1,:] = np.fft.fftshift(10*np.log10(Pxx_den)) + + self.spectrum[0,:] = np.fft.fftshift(f) + + + # Create signal window for plot +# signal_window = np.ones(len(self.spectrum[1,:])) * -100 + # signal_window[max(self.max_index - self.fft_signal_width//2, 0) : min(self.max_index + self.fft_signal_width//2, len(self.spectrum[1,:]))] = max(self.spectrum[1,:]) + #signal_window = np.ones(len(max_spectrum)) * -100 + #signal_window[max(self.max_index - self.fft_signal_width//2, 0) : min(self.max_index + self.fft_signal_width//2, len(max_spectrum))] = max(max_spectrum) + + #self.spectrum[self.channel_number+1, :] = signal_window #np.ones(len(spectrum[1,:])) * self.module_receiver.daq_squelch_th_dB # Plot threshold line + que_data_packet.append(['spectrum', self.spectrum]) + + #-----> Passive Radar <----- + conf_val = 0 + theta_0 = 0 + if self.en_PR and self.data_ready and self.channel_number > 1: + + ref_ch = self.module_receiver.iq_samples[0,:] + surv_ch = self.module_receiver.iq_samples[1,:] + + td_filter_dimension = self.max_bistatic_range + + start = time.time() + + if self.PR_clutter_cancellation == "Wiener MRE": + surv_ch, w = Wiener_SMI_MRE(ref_ch, surv_ch, td_filter_dimension) + #surv_ch, w = cc.Wiener_SMI_MRE(ref_ch, surv_ch, td_filter_dimension) + + surv_ch = det.windowing(surv_ch, "Hamming") #surv_ch * signal.tukey(surv_ch.size, alpha=0.25) #det.windowing(surv_ch, "hamming") + + max_Doppler = self.max_doppler #256 + max_range = self.max_bistatic_range + + #RD_matrix = det.cc_detector_ons(ref_ch, surv_ch, self.module_receiver.iq_header.sampling_freq, max_Doppler, max_range, verbose=0, Qt_obj=None) + RD_matrix = cc_detector_ons(ref_ch, surv_ch, self.module_receiver.iq_header.sampling_freq, max_Doppler, max_range) + + end = time.time() + print("Time: " + str((end-start) * 1000)) + + que_data_packet.append(['RD_matrix', RD_matrix]) + + # Record IQ samples + if self.en_record: + # TODO: Implement IQ frame recording + self.logger.error("Saving IQ samples to npy is obsolete, IQ Frame saving is currently not implemented") + + stop_time = time.time() + que_data_packet.append(['update_rate', stop_time-start_time]) + que_data_packet.append(['latency', int(stop_time*10**3)-self.module_receiver.iq_header.time_stamp]) + + # If the que is full, and data is ready (from squelching), clear the buffer immediately so that useful data has the priority + if self.data_que.full() and self.data_ready: + try: + #self.logger.info("BUFFER WAS NOT EMPTY, EMPTYING NOW") + self.data_que.get(False) #empty que if not taken yet so fresh data is put in + except queue.Empty: + #self.logger.info("DIDNT EMPTY") + pass + + # Put data into buffer, but if there is no data because its a cal/trig wait frame etc, then only write if the buffer is empty + # Otherwise just discard the data so that we don't overwrite good DATA frames. + try: + self.data_que.put(que_data_packet, False) # Must be non-blocking so DOA can update when dash browser window is closed + except: + # Discard data, UI couldn't consume fast enough + pass + + """ + start = time.time() + end = time.time() + thetime = ((end - start) * 1000) + print ("Time elapsed: ", thetime) + """ +@jit(fastmath=True) +def Wiener_SMI_MRE(ref_ch, surv_ch, K): + """ + Description: + ------------ + Performs Wiener filtering with applying the Minimum Redundance Estimation (MRE) technique. + When using MRE, the autocorrelation matrix is not fully estimated, but only the first column. + With this modification the required calculations can be reduced from KxK to K element. + + Parameters: + ----------- + :param K : Filter tap number + :param ref_ch : Reference signal array + :param surv_ch: Surveillance signal array + + :type K : int + :type ref_ch : 1 x N complex numpy array + :type surv_ch: 1 x N complex numpy array + Return values: + -------------- + :return filt: Filtered surveillance channel + :rtype filt: 1 x N complex numpy array + + :return None: Input parameters are not consistent + """ + + N = ref_ch.shape[0] # Number of time samples + R, r = pruned_correlation(ref_ch, surv_ch, K, N) + R_mult = R_eye_memoize(K) + w = fast_w(R, r, K, R_mult) + + #return surv_ch - np.convolve(ref_ch, w)[0:N], w # subtract the zero doppler clutter + return surv_ch - signal.oaconvolve(ref_ch, w)[0:N], w # subtract the zero doppler clutter #oaconvolve saves us about 100-200 ms + +@njit(fastmath=True, parallel=True, cache=True) +def fast_w(R, r, K, R_mult): + # Complete the R matrix based on its Hermitian and Toeplitz property + + for k in range(1, K): + R[:, k] = shift(R[:, 0], k) + #R[:, K] = shift(R[:,0], K) + + R += np.transpose(np.conjugate(R)) + R *= R_mult #(np.ones(K) - np.eye(K) * 0.5) + + #w = np.dot(lin.inv(R), r) # weight vector + w = lin.inv(R) @ r #np.dot(lin.inv(R), r) # weight vector #matmul (@) may be slightly faster that np.dot for 1D, 2D arrays. + # inverse and dot product run time : 1.1s for 2048*2048 matrix + + return w + +#Memoize ~50ms speedup? +@lru_cache(maxsize=2) +def R_eye_memoize(K): + return (np.ones(K) - np.eye(K) * 0.5) + +#Modified pruned correlation, returns R and r directly and saves one FFT +@jit(fastmath=True, cache=True) +def pruned_correlation(ref_ch, surv_ch, clen, N): + """ + Description: + ----------- + Calculates the part of the correlation function of arrays with same size + The total length of the cross-correlation function is 2*N-1, but this + function calculates the values of the cross-correlation between [N-1 : N+clen-1] + Parameters: + ----------- + :param x : input array + :param y : input array + :param clen: correlation length + + :type x: 1 x N complex numpy array + :type y: 1 x N complex numpy array + :type clen: int + Return values: + -------------- + :return corr : part of the cross-correlation function + :rtype corr : 1 x clen complex numpy array + + :return None : inconsistent array size + """ + R = np.zeros((clen, clen), dtype=c_dtype) # Autocorrelation mtx. + + # --calculation-- + # set up input matrices pad zeros if not multiply of the correlation length + cols = clen - 1 #(clen = Filter drowsimension) + rows = np.int32(N / (cols)) + 1 + + zeropads = cols * rows - N + x = np.pad(ref_ch, (0, zeropads)) + + # shaping inputs into matrices + xp = np.reshape(x, (rows, cols)) + + # padding matrices for FFT + ypp = np.vstack([xp[1:, :], np.zeros(cols, dtype=c_dtype)]) #vstack appears to be faster than pad + yp = np.concatenate([xp, ypp], axis=1) + + # execute FFT on the matrices + xpw = fft.fft(xp, n = 2*cols, axis=1, workers=4, overwrite_x=True) + bpw = fft.fft(yp, axis=1, workers=4, overwrite_x=True) + + # magic formula which describes the unified equation of the universe + # corr_batches = np.fliplr(fft.fftshift(fft.ifft(corr_mult(xpw, bpw), axis=1, workers=4, overwrite_x=True)).conj()[:, 0:clen]) + corr_batches = fft.fftshift(fft.ifft(corr_mult(xpw, bpw), axis=1, workers=4, overwrite_x=True)).conj()[:, 0:clen] + + # sum each value in a column of the batched correlation matrix + R[:,0] = np.fliplr([np.sum(corr_batches, axis=0)])[0] + + #calc r + y = np.pad(surv_ch, (0, zeropads)) + yp = np.reshape(y, (rows, cols)) + ypp = np.vstack([yp[1:, :], np.zeros(cols, dtype=c_dtype)]) #vstack appears to be faster than pad + yp = np.concatenate([yp, ypp], axis=1) + bpw = fft.fft(yp, axis=1, workers=4, overwrite_x=True) + #corr_batches = np.fliplr(fft.fftshift(fft.ifft(corr_mult(xpw, bpw), axis=1, workers=4, overwrite_x=True)).conj()[:, 0:clen]) + corr_batches = fft.fftshift(fft.ifft(corr_mult(xpw, bpw), axis=1, workers=4, overwrite_x=True)).conj()[:, 0:clen] + #r = np.sum(corr_batches, axis=0) + r = np.fliplr([np.sum(corr_batches, axis=0)])[0] + + return R, r + +@njit(fastmath=True, cache=True) +def shift(x, i): + """ + Description: + ----------- + Similar to np.roll function, but not circularly shift values + Example: + x = |x0|x1|...|xN-1| + y = shift(x,2) + x --> y: |0|0|x0|x1|...|xN-3| + Parameters: + ----------- + :param:x : input array on which the roll will be performed + :param i : delay value [sample] + + :type i :int + :type x: N x 1 complex numpy array + Return values: + -------------- + :return shifted : shifted version of x + :rtype shifted: N x 1 complex numpy array + """ + + N = x.shape[0] + if np.abs(i) >= N: + return np.zeros(N, dtype=c_dtype) + if i == 0: + return x + shifted = np.roll(x, i) + if i < 0: + shifted[np.mod(N + i, N):] = np.zeros(np.abs(i), dtype=c_dtype) + if i > 0: + shifted[0:i] = np.zeros(np.abs(i), dtype=c_dtype) + return shifted + + +@njit(fastmath=True, parallel=True, cache=True) +def resize_and_align(no_sub_tasks, ref_ch, surv_ch, fs, fD_max, r_max): + surv_ch_align = np.reshape(surv_ch,(no_sub_tasks, r_max)) # shaping surveillance signal array into a matrix + pad_zeros = np.expand_dims(np.zeros(r_max, dtype=c_dtype), axis=0) + surv_ch_align = np.vstack((surv_ch_align, pad_zeros)) # padding one row of zeros into the surv matrix + surv_ch_align = np.concatenate((surv_ch_align[0 : no_sub_tasks,:], surv_ch_align[1 : no_sub_tasks +1, :]), axis = 1) + + ref_ch_align = np.reshape(ref_ch, (no_sub_tasks, r_max)) # shaping reference signal array into a matrix + pad_zeros = np.zeros((no_sub_tasks, r_max),dtype = c_dtype) + ref_ch_align = np.concatenate((ref_ch_align, pad_zeros),axis = 1) # shaping + + return ref_ch_align, surv_ch_align + +@njit(fastmath=True, cache=True) +def corr_mult(surv_fft, ref_fft): + return np.multiply(surv_fft, ref_fft.conj()) + +@jit(fastmath=True, cache=True) +def cc_detector_ons(ref_ch, surv_ch, fs, fD_max, r_max): + """ + Parameters: + ----------- + :param N: Range resolution - N must be a divisor of the input length + :param F: Doppler resolution, F has a theoretical limit. If you break the limit, the output may repeat + itself and get wrong results. F should be less than length/N otherwise use other method! + Return values: + -------------- + :return None: Improper input parameters + + """ + N = ref_ch.size + + # --> Set processing parameters + fD_step = fs / (2 * N) # Doppler frequency step size (with zero padding) + Doppler_freqs_size = int(fD_max / fD_step) + no_sub_tasks = N // r_max + + # Allocate range-Doppler maxtrix + mx = np.zeros((2*Doppler_freqs_size+1, r_max),dtype = c_dtype) #memoize_zeros((2*Doppler_freqs_size+1, r_max), c_dtype) #np.zeros((2*Doppler_freqs_size+1, r_max),dtype = nb.c8) + + ref_ch_align, surv_ch_align = resize_and_align(no_sub_tasks, ref_ch, surv_ch, fs, fD_max, r_max) + + # row wise fft on both channels + ref_fft = fft.fft(ref_ch_align, axis = 1, overwrite_x=True, workers=4) #pyfftw.interfaces.numpy_fft.fft(ref_ch_align_a, axis = 1, overwrite_input=True, threads=4) #fft.fft(ref_ch_align_a, axis = 1, overwrite_x=True, workers=4) + surv_fft = fft.fft(surv_ch_align, axis = 1, overwrite_x=True, workers=4) #pyfftw.interfaces.numpy_fft.fft(surv_ch_align_a, axis = 1, overwrite_input=True, threads=4) #fft.fft(surv_ch_align_a, axis = 1, overwrite_x=True, workers=4) + + corr = corr_mult(surv_fft, ref_fft) #np.multiply(surv_fft, ref_fft.conj()) + + corr = fft.ifft(corr,axis = 1, workers=4, overwrite_x=True) + + corr_a = pyfftw.empty_aligned(np.shape(corr), dtype=c_dtype) + corr_a[:] = corr #.copy() + + # This is the most computationally intensive part ~120ms, overwrite_x=True gives a big speedup, not sure if it changes the result though... + corr = fft.fft(corr_a, n=2* no_sub_tasks, axis = 0, workers=4, overwrite_x=True) # Setting the output size with "n=.." is faster than doing a concat first. + + # crop and fft shift + mx[ 0 : Doppler_freqs_size, 0 : r_max] = corr[2*no_sub_tasks - Doppler_freqs_size : 2*no_sub_tasks, 0 : r_max] + mx[Doppler_freqs_size : 2 * Doppler_freqs_size+1, 0 : r_max] = corr[ 0 : Doppler_freqs_size+1 , 0 : r_max] + + return mx + + + + + + +#NUMBA optimized center tracking. Gives a mild speed boost ~25% faster. +@njit(fastmath=True, cache=True, parallel=True) +def center_max_signal(processed_signal, frequency, fft_spectrum, threshold, sample_freq): + + # Where is the max frequency? e.g. where is the signal? + max_index = np.argmax(fft_spectrum) + max_frequency = frequency[max_index] + + # Auto decimate down to exactly the max signal width + fft_signal_width = np.sum(fft_spectrum > threshold) + 25 + decimation_factor = max((sample_freq // fft_signal_width) // 2, 1) + + # Auto shift peak frequency center of spectrum, this frequency will be decimated: + # https://pysdr.org/content/filters.html + f0 = -max_frequency #+10 + Ts = 1.0/sample_freq + t = np.arange(0.0, Ts*len(processed_signal[0, :]), Ts) + exponential = np.exp(2j*np.pi*f0*t) # this is essentially a complex sine wave + + return processed_signal * exponential, decimation_factor, fft_signal_width, max_index + + +# NUMBA optimized MUSIC function. About 100x faster on the Pi 4 +@njit(fastmath=True, cache=True) +def DOA_MUSIC(R, scanning_vectors, signal_dimension, angle_resolution=1): + # --> Input check + if R[:,0].size != R[0,:].size: + print("ERROR: Correlation matrix is not quadratic") + return np.ones(1, dtype=nb.c16)*-1 #[(-1, -1j)] + + if R[:,0].size != scanning_vectors[:,0].size: + print("ERROR: Correlation matrix dimension does not match with the antenna array dimension") + return np.ones(1, dtype=nb.c16)*-2 + + #ADORT = np.zeros(scanning_vectors[0,:].size, dtype=np.complex) #CHANGE TO nb.c16 for NUMBA + ADORT = np.zeros(scanning_vectors[0,:].size, dtype=nb.c16) + M = R[:,0].size #np.size(R, 0) + + # --- Calculation --- + # Determine eigenvectors and eigenvalues + sigmai, vi = lin.eig(R) + sigmai = np.abs(sigmai) + + idx = sigmai.argsort()[::1] # Sort eigenvectors by eigenvalues, smallest to largest + #sigmai = sigmai[idx] # Eigenvalues not used again + vi = vi[:,idx] + + # Generate noise subspace matrix + noise_dimension = M - signal_dimension + #E = np.zeros((M, noise_dimension),dtype=np.complex) + E = np.zeros((M, noise_dimension),dtype=nb.c16) + for i in range(noise_dimension): + E[:,i] = vi[:,i] + + theta_index=0 + for i in range(scanning_vectors[0,:].size): + S_theta_ = scanning_vectors[:, i] + S_theta_ = S_theta_.T + ADORT[theta_index] = 1/np.abs(S_theta_.conj().T @ (E @ E.conj().T) @ S_theta_) + theta_index += 1 + + return ADORT + +# Numba optimized version of pyArgus corr_matrix_estimate with "fast". About 2x faster on Pi4 +@njit(fastmath=True, cache=True) #(nb.c8[:,:](nb.c16[:,:])) +def corr_matrix(X): + M = X[:,0].size + N = X[0,:].size + #R = np.zeros((M, M), dtype=nb.c8) + R = np.dot(X, X.conj().T) + R = np.divide(R, N) + return R + +# Numba optimized scanning vectors generation for UCA arrays. About 10x faster on Pi4 +# LRU cache memoize about 1000x faster. +@lru_cache(maxsize=8) +def uca_scanning_vectors(M, DOA_inter_elem_space): + + thetas = np.linspace(0,359,360) # Remember to change self.DOA_thetas too, we didn't include that in this function due to memoization cannot work with arrays + + x = DOA_inter_elem_space * np.cos(2*np.pi/M * np.arange(M)) + y = -DOA_inter_elem_space * np.sin(2*np.pi/M * np.arange(M)) # For this specific array only + + scanning_vectors = np.zeros((M, thetas.size), dtype=np.complex) + for i in range(thetas.size): + scanning_vectors[:,i] = np.exp(1j*2*np.pi* (x*np.cos(np.deg2rad(thetas[i])) + y*np.sin(np.deg2rad(thetas[i])))) + + return scanning_vectors + # scanning_vectors = de.gen_scanning_vectors(M, x, y, self.DOA_theta) + +@njit(fastmath=True, cache=True) +def DOA_plot_util(DOA_data, log_scale_min=-100): + """ + This function prepares the calulcated DoA estimation results for plotting. + + - Noramlize DoA estimation results + - Changes to log scale + """ + + DOA_data = np.divide(np.abs(DOA_data), np.max(np.abs(DOA_data))) # Normalization + DOA_data = 10*np.log10(DOA_data) # Change to logscale + + for i in range(len(DOA_data)): # Remove extremely low values + if DOA_data[i] < log_scale_min: + DOA_data[i] = log_scale_min + + return DOA_data + +@njit(fastmath=True, cache=True) +def calculate_doa_papr(DOA_data): + return 10*np.log10(np.max(np.abs(DOA_data))/np.mean(np.abs(DOA_data))) + +# Old time-domain squelch algorithm (Unused as freq domain FFT with overlaps gives significantly better sensitivity with acceptable time resolution expense +""" + K = 10 + self.filtered_signal = self.raw_signal_amplitude #convolve(np.abs(self.raw_signal_amplitude),np.ones(K), mode = 'same')/K + + # Burst is always started at the begining of the processed block, ensured by the squelch module in the DAQ FW + burst_stop_index = len(self.filtered_signal) # CARL FIX: Initialize this to the length of the signal, incase the signal is active the entire time + self.logger.info("Original burst stop index: {:d}".format(burst_stop_index)) + + min_burst_size = K + burst_stop_amp_val = 0 + for n in np.arange(K, len(self.filtered_signal), 1): + if self.filtered_signal[n] < self.squelch_threshold: + burst_stop_amp_val = self.filtered_signal[n] + burst_stop_index = n + burst_stop_index-=K # Correction with the length of filter + break + + #burst_stop_index-=K # Correction with the length of filter + + + self.logger.info("Burst stop index: {:d}".format(burst_stop_index)) + self.logger.info("Burst stop ampl val: {:f}".format(burst_stop_amp_val)) + self.logger.info("Processed signal length: {:d}".format(len(self.processed_signal[0,:]))) + + # If sign + if burst_stop_index < min_burst_size: + self.logger.debug("The length of the captured burst size is under the minimum: {:d}".format(burst_stop_index)) + burst_stop_index = 0 + + if burst_stop_index !=0: + self.logger.info("INSIDE burst_stop_index != 0") + + self.logger.debug("Burst stop index: {:d}".format(burst_stop_index)) + self.logger.debug("Burst stop ampl val: {:f}".format(burst_stop_amp_val)) + self.squelch_mask = np.zeros(len(self.filtered_signal)) + self.squelch_mask[0 : burst_stop_index] = np.ones(burst_stop_index)*self.squelch_threshold + # Next line removes the end parts of the samples after where the signal ended, truncating the array + self.processed_signal = self.module_receiver.iq_samples[: burst_stop_index, self.squelch_mask == self.squelch_threshold] + self.logger.info("Raw signal length when burst_stop_index!=0: {:d}".format(len(self.module_receiver.iq_samples[0,:]))) + self.logger.info("Processed signal length when burst_stop_index!=0: {:d}".format(len(self.processed_signal[0,:]))) + + #self.logger.info(' '.join(map(str, self.processed_signal))) + + self.data_ready=True + else: + self.logger.info("Signal burst is not found, try to adjust the threshold levels") + #self.data_ready=True + self.squelch_mask = np.ones(len(self.filtered_signal))*self.squelch_threshold + self.processed_signal = np.zeros([self.channel_number, len(self.filtered_signal)]) +""" + + |