# -*- coding: utf-8 -*-
"""
Radar Class.
Contains Radar class and a RadarTypes class
Contains RadarArray class to hold timeseries data for voltage measurements,
FFT slices, and max frequency.
Author: Jason Merlo
Maintainer: Jason Merlo (merlojas@msu.edu)
"""
import numpy as np # Storing data
from pyratk.datatypes.ts_data import TimeSeries # storing data
import scipy.constants as spc # speed of light
from scipy import signal # Upsampling, bh window
from pyratk.datatypes.geometry import Point # radar location
# from pyratk.datatypes.motion import StateMatrix # track location
from pyqtgraph import QtCore
from pyratk.formatting import warning
from profilehooks import profile
from collections import namedtuple
class Receiver(object):
"""
Object to hold Radar Receiver information and signal processing.
It is assumed radars will have both I and Q channels.
NOTE: Radar class should not contain state history (TimeSeries).
Work on removing these and placing them in Tracker or DataLogger.
"""
# === INITIALIZATION METHODS ==============================================
def __init__(self, daq, daq_index, transmitter, loc=Point(),
fast_fft_size=2**12, fast_fft_window_type=None,
slow_fft_size=2**12, slow_fft_window_type=None,
slow_fft_len=20):
super().__init__()
"""
Initialize radar object.
Note: Currently only FMCW or CW radars are supported.
Args:
daq
[daq] object for DAQ configurations
daq_index
Index of radar I & Q channels in DAQ output
transmitter
Transmitter object
loc
Point() indicating the location of the radar
"""
# Data stream parameters
# TODO: create DataStream object representing I or I/Q data
self.daq = daq
self.daq_index = daq_index
# Physical parameters
self.transmitter = transmitter
self.loc = loc
# Processing parameters
self.fast_fft_size = fast_fft_size
self.fast_fft_window_type = fast_fft_window_type
self.slow_fft_size = slow_fft_size
self.slow_fft_window_type = slow_fft_window_type
self.slow_fft_len = slow_fft_len
self.datacube = DataCube(self)
self.init_data()
self.connect_signals()
def init_data(self):
# Derived processing parameters
self.fast_center_bin = np.ceil(self.fast_fft_size / 2)
self.fast_bin_size = self.daq.sample_rate / self.fast_fft_size
self.slow_center_bin = np.ceil(self.slow_fft_size / 2)
self.slow_bin_size = self.daq.sample_rate * self.transmitter.pulses[0].delay / self.slow_fft_size
self.fast_fft_data = np.ones(self.fast_fft_size)
self.fft_mat = np.ones((self.fast_fft_size, self.slow_fft_size))
self.fast_fmax = 0
self.slow_fmax = 0
self.data = None
self.fast_fft_len=int(round(self.daq.sample_rate * self.transmitter.pulses[0].delay))
self.mti_window = np.transpose(np.tile(np.fft.fftshift(signal.windows.chebwin(self.slow_fft_len,at=60)),self.fast_fft_size).reshape((-1,self.slow_fft_len)))
def connect_signals(self):
# self.daq.reset_signal.connect(self.reset)
pass
# === HELPER METHODS ======================================================
def freq_to_vel(self, freq):
"""Compute the velocity for a given frequency and the radar fc."""
c = spc.speed_of_light
# velocity = (c * self.fc / (freq + self.fc)) - c
velocity = (freq * c) / (2 * self.fc)
return velocity
def bin_to_freq(self, bin):
"""Compute frequency based on bin location."""
return (bin - self.fast_center_bin) * float(self.fast_bin_size)
def compute_cfft(self, complex_data, fft_size):
"""Compute fft and fft magnitude for plotting."""
# Create blackman-harris window
# window = signal.blackmanharris(complex_data.shape[0])
# window = signal.hamming(complex_data.shape[0])
window = np.ones(complex_data.shape[0])
# Create zero-padded array to be transformed
fft_array = np.zeros((fft_size,), dtype=np.complex64)
fft_array[:complex_data.size] = complex_data * window
# Normalize FFT magnitude to window size
fft_complex = np.fft.fft(fft_array, norm='ortho')
# Adjust fft so DC is at the center
fft_complex = np.fft.fftshift(fft_complex)
# Display only magnitude
fft_mag = np.linalg.norm([fft_complex.real, fft_complex.imag], axis=0)
return fft_mag
def compute_fft2(self, data, shape):
"""
Compute 2D FFT over range and Doppler.
data - data to be transformed
shape - output shape of transform (with zero-padding)
"""
# Create blackman-harris window
# window = signal.blackmanharris(complex_data.shape[0])
# window = signal.hamming(complex_data.shape[0])
window = np.ones(data.shape)
# Create zero-padded array to be transformed
fft_array = np.zeros(shape, dtype=np.complex64)
fft_array[:data.shape[0], :data.shape[1]] = data * window
# Normalize FFT magnitude to window size
fft_complex = np.fft.fft2(fft_array, norm='ortho')
# Adjust fft so DC is at the center
fft_complex = np.fft.fftshift(fft_complex)
# Display only magnitude
fft_mag = np.linalg.norm([fft_complex.real, fft_complex.imag], axis=0)
return fft_mag
# === CONTROL METHODS =====================================================
def update(self, data):
# Get data from daq
self.data = data_slice = np.array((data[self.daq_index[0]], data[self.daq_index[1]]))
#
# # self.data_buffer = np.append(self.data_buffer, iq_data_slice)
#
# # Get window of FFT data
# window_idx = self.fast_fft_size // self.daq.sample_chunk_size
# window_slice = data_slice[0, :].flatten() \
# + data_slice[ 1, :].flatten() * 1.0j
# # window_slice = window_slice_pair[:, 1, :].flatten()
#
# # Subtract any DC component
# window_slice -= np.mean(window_slice.real) + np.mean(window_slice.imag) * 1.0j
#
#
# # Calculate fast-time complex FFT
# self.fast_fft_data = self.compute_cfft(window_slice, self.fast_fft_size)
#
# # Find maximum frequency
# fmax_bin = np.argmax(self.fast_fft_data)
# self.fast_fmax = self.bin_to_freq(fmax_bin)
# Calculate slow-time complex FFT
# print('dc.shape',dc.shape)
self.fft_mat = self.compute_fft2(self.datacube[-1], (self.slow_fft_size, self.fast_fft_size))
self.fft_mat=np.multiply(self.fft_mat,self.mti_window)
# print('fft_mat.shape', self.fft_mat.shape)
# if self.datacube[-1].shape == self.datacube[-2].shape:
# if hasattr(self, 'zero_fft_mat'):
# self.fft_mat -= self.zero_fft_mat
# else:
# self.zero_fft_mat = self.fft_mat
# Power Thresholding
# if self.cfft_data[vmax_bin] < POWER_THRESHOLD:
# self.fmax = 0
# else:
# self.fmax = self.bin_to_freq(vmax_bin)
# self.drho = self.freq_to_vel(self.fmax)
def reset(self):
"""Reset all radar data."""
self.init_data()
class DataCube(object):
def __init__(self, receiver):
self.receiver = receiver
# TODO: only supports single pulse
self.samples_per_pulse = int(self.receiver.daq.sample_rate * self.receiver.transmitter.pulses[0].delay)
def get_frame(self, key):
"""
Create datacube with specified number of most recent frames.
Returns a complex datacube. Final shape will be:
slow_fft_len x samples_per_pulse
"""
start = self.receiver.slow_fft_len * key
end = self.receiver.slow_fft_len * (key + 1)
idx = self.receiver.daq_index
if end == 0:
datacube = self.receiver.daq.ts_buffer[start:, idx[0]] \
+ 1.0j * self.receiver.daq.ts_buffer[start:, idx[1]]
else:
datacube = self.receiver.daq.ts_buffer[start:end, idx[0]] \
+ 1.0j * self.receiver.daq.ts_buffer[start:end, idx[1]]
# print('datacube.shape:', datacube.shape)
# Subtract mean
mean_i = np.mean(datacube.real)
mean_q = np.mean(datacube.imag)
datacube -= mean_i + mean_q *1.0j
return datacube
def __getitem__(self, key):
if isinstance(key, slice):
# start, stop, step = key.indices(len(self))
# return Seq([self[i] for i in range(start, stop, step)])
raise NotImplementedError('Slicing multiple frames of datacube not yet implemented.')
elif isinstance(key, int):
return self.get_frame(key)
elif isinstance(key, tuple):
raise NotImplementedError('Tuple as index')
else:
raise TypeError('Invalid argument type: {}'.format(type(key)))
class Transmitter(object):
"""
Object to hold Radar Transmitter information.
Note: This class may perform transmitter control functions in the future.
"""
def __init__(self, data_mgr, pulses, loc=Point()):
super().__init__()
"""
Initialize radar object.
Note: Currently only FMCW or CW radars are supported.
Args:
data_mgr
[data_mgr] object for DAQ configurations
loc
Point() indicating the location of the radar
pulses
List of pulse namedtuple objects containing fc, bw, and delay
fc
Center frequency of transmit waveform
bw
Bandwidth of transmit chirp waveform
delay
Delay time between start of chirp
"""
self.data_mgr = data_mgr
self.pulses = pulses
# Currently only one pulse is supported due to datacube restrictions
if len(self.pulses) != 1:
raise NotImplementedError('Only one pulse segment is currently'
' supported due to datacube processing restrictions')
class Radar(QtCore.QObject):
"""
Holds radar Transmitter and Receiver objects.
Note: all iterable functions still only iterate over receiver for
backwards compatibility with old RadarArray structure.
Attributes:
tranmitters (namedtuple)
location (Point)
pulses (namedtuple)
fc
Pulse waveform center frequency
bw
Pulse waveform bandwidth
delay
Time between start of waveform modulation
receivers (namedtuple)
daq_index (tuple)
tuple containing index of I and Q channels on DAQ
location (Point)
fast_fft_size
fast_fft_window_type
slow_fft_size
slow_fft_window_type
slow_fft_len
"""
data_available_signal = QtCore.pyqtSignal()
def __init__(self,
daq,
transmitter_list,
receiver_list,
fast_fft_size=2**10,
fast_fft_window_type=None,
slow_fft_size=2**12,
slow_fft_window_type=None,
slow_fft_len=100
):
super().__init__()
self.daq = daq
# Used for iterator magic functions
self.idx = 0
# Flag to iqnore any stale data after a reset
self.last_sample_index = -1
# Create recievers and transmitters
self.transmitters = []
self.receivers = []
for transmitter in transmitter_list:
self.transmitters.append(
Transmitter(
self.daq,
transmitter.pulses
)
)
for receiver in receiver_list:
self.receivers.append(
Receiver(
self.daq,
receiver.daq_index,
self.transmitters[0],
receiver.location,
fast_fft_size,
fast_fft_window_type,
slow_fft_size,
slow_fft_window_type,
slow_fft_len
)
)
# Configure Signals
self.connect_signals()
def connect_signals(self):
"""Connect signals for event driven operation."""
self.daq.data_available_signal.connect(self.update)
self.daq.reset_signal.connect(self.reset)
def reset(self):
"""Clear all temporal data from radars in array."""
for receiver in self.receivers:
receiver.reset()
self.last_sample_index = -1
# @profile(immediate=True)
def update(self, data_tuple):
"""Update all radars in array."""
# start_time = time.time()
data, sample_index = data_tuple
# print('(radar.py) sample_num:', sample_index)
if sample_index > self.last_sample_index:
self.reset_flag = False
for receiver in self.receivers:
receiver.update(data)
# Emit data available signal for dependant tasks
self.data_available_signal.emit()
# print('Updating radar, current_index: ', sample_index)
else:
# print('(radar.py) last_sample_index: {}\tcurrent_index:{}'.format(self.last_sample_index, sample_index))
# warning('(radar.py) ignoring stale data...')
pass
self.last_sample_index = sample_index
# print('(radar.py) radar_array.update() ran in {:} (s)'
# .format(time.time() - start_time))
def __getitem__(self, i):
"""Return radar in array."""
assert(i < len(self)), "index greater than flattened array size"
return self.receivers[i]
def __len__(self):
"""Return length of radar array."""
# compute sumproduct of array shape
return len(self.receivers)
def __iter__(self):
"""Return self for iterator."""
return self
def __next__(self):
"""Return next radar object for iterator."""
if self.idx < len(self):
self.idx += 1
return self.receivers[self.idx - 1]
else:
self.idx = 0
raise StopIteration()