summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-x_signal_processing/krakenSDR_signal_processor.py177
1 files changed, 7 insertions, 170 deletions
diff --git a/_signal_processing/krakenSDR_signal_processor.py b/_signal_processing/krakenSDR_signal_processor.py
index 651fc39..d0b6dba 100755
--- a/_signal_processing/krakenSDR_signal_processor.py
+++ b/_signal_processing/krakenSDR_signal_processor.py
@@ -421,6 +421,9 @@ def pruned_correlation(ref_ch, surv_ch, clen, N):
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)
+ print("pruned corr xp shape: " + str(xp.shape))
+ print("pruned corr yp shape: " + str(yp.shape))
+
# 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)
@@ -523,6 +526,10 @@ def cc_detector_ons(ref_ch, surv_ch, fs, fD_max, r_max):
ref_ch_align, surv_ch_align = resize_and_align(no_sub_tasks, ref_ch, surv_ch, fs, fD_max, r_max)
+ print("ref_ch_align shape: " + str(ref_ch_align.shape))
+ print("surv_ch_align shape: " + str(surv_ch_align.shape))
+
+
# 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)
@@ -543,173 +550,3 @@ def cc_detector_ons(ref_ch, surv_ch, fs, fD_max, 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)])
-"""
-
-