This post documents the procedure from acquiring noisy, distorted CW radar data through drawing conclusions for the next steps in Red Pitaya radar development.
Read the data into the computer memory appropriately.
GNU Radio complex is single-precision float, marshaled as real/imag pairs of 4 byte floats.
Matlab: read GNU Radio complex data with read_complex_binary.m that reshapes and casts the data appropriately.
Python: numpy.fromfile(filename,'complex64')
An example command from the piradar program is
./PlotSpectrum.py exercise.bin 100e3
where we see the 60 Hz periodic interference, which we model in the next section.
If the sinusoids are not clearly and cleanly present, or the audio playback doesn’t sound like a pure sinusoid, something is very wrong–stop at this point and fix the problem.
Signal defects may come from software and/or hardware.
They may have subtle or visible effects.
An empirical model was created based on the appearance of ADC clipping in the time and frequency domain in the piradar program SimADCclipping.py
./SimADCclipping.py
This model gives the basis to build filters suitable for preprocessing before using subspace estimation of close-spaced CW sinusoids typical to CW radar with Red Pitaya or any other such instrument.
Without filtering, the subspace estimator homes in on the forest of 60 Hz “fake” targets.
One can create remarkable filters in the simulation or analytic space.
Whether they can be realized in a real system depends on the CPU/FPGA power available in the real system and if they require excessively precise frequency/phase control and the like.
Here’s let’s assume we can use a FIR filter with order 4096 taps to reduce signals greater than +/- 60 Hz from the CW radar center frequency by about 20 dB and much more reduction beyond that.
The center frequency chosen is for a radar baseband center frequency of 10 kHz.
The text file cwfir.asc contains FIR filter coefficients used in the next step.
FIR filters are used in perhaps the vast majority of electronic devices.
In contrast with IIR filters, which are generally realized with far fewer taps (less computational effort) for a given specification, FIR filters have linear phase in-band, while IIR filters in general do not.
The signal_subspace program has both a Python and Fortran implementation of a FIR filter.
However, there is also a filter built into scipy.signal.decimate.
This offline analysis used SciPy’s filter.
After downconversion, the signals look like the following, with and without target.
How much to downconvert is driven by the signal bandwidth.
Downconversion will act like a shift register for the analytic signal in frequency.
Downconverting the signal might take place most effectively in multiple steps if extreme ratios are involved, such as might occur in software defined radar.
This is where RFNoC or custom FPGA programming can conserve immense amount of data bandwidth and CPU/RAM resources by downconverting in the FPGA.
If the desired signal bandwidth is 50 Hz or 5 Hz, in general you don’t need the immensely higher sample rate that is the lowest sample rate the software defined radar is capable of providing.
So even on the CPU portion of the processing, further downsampling might conserve further CPU/RAM resources.
Particular analyses such as those using the Goertzel algorithm might benefit from particular sample rate choices.
{% include gallery id=“gallery”
caption=“Red Pitaya CW Radar at 47 MHz, mostly stationary human targets.” %}
{% include gallery id=“gallery2”
caption=“Red Pitaya CW Radar at 47 MHz, cyclically moving human target at 2-3 meters standoff distance.” %}
Human activity generates a range of Doppler frequencies in general for any radar.
If measuring on very short time scales relative to human motion, the Doppler frequency extent may still be finite due to arms & legs swinging vs. body trunk motion.
Simple products such as motion-detecting lights in the 2.4 and 5.8 GHz global license-free bands use analog filters.
One might also consider DFT-based methods, and look for the expected frequency bins to exceed a certain power threshold as a detection.
ESPRIT subspace estimation in signal_subspace/tests/test.py:test_esprit() is approximately O(M3), where M is “block length”.
The block length parameter M informally controls the amount of CPU power applied to find the sinusoids.
Lower SNR or more closely spaced sinusoids require a larger M, with the penalty that computation time increases O(M3).
However, where there are not pure sinusoids, such as the range of frequencies generated by human motion, alternative target detection methods may be more appropriate.
The CW mode typically measures moving targets only.
The targets manifest as sinusoids displaced in frequency–positive shift for targets approaching the radar antenna, and downward shift for targets moving away from the antenna.
When looking at a plot of the time-domain radar signal, zoomed out far enough to see several cycles of the Doppler frequency (e.g. several seconds displayed), the frequency differences show up as a sinusoidal varying of the signal envelope.
One could imagine 4 or 5 frequency bins, from low to high frequency, something like:
lowest frequency bin must be discarded, as it will pickup the noise of the transmitter, and very slowly moving objects like the window shades moving in the breeze.
someone breathing, rocking/turning in their seat, or exercising.
someone walking
someone running
fans spinning, cars driving by.
The frequency and width of the bins is chosen based on the user application.
DFT, FFT, Goertzel, etc. may be suitable analysis methods.
CFAR methods can be used in the frequency bins to mitigate false measurements.
Activity detection vs. time with CW radar at 47 MHz.
The subspace estimator outputs two vectors.
The first vector is the frequency estimate and the second vector is the corresponding singular value.
The larger the singular value, the higher the confidence in the frequency estimate.
This can be used as a qualifier to discard unlikely targets.
Also, targets outside a plausible Doppler shift (or range, for FMCW) can be discarded by inspection.
During daylight hours (and other more occasional events) D-region absorption increases immensely for VLF, LF, MW, and lower HF frequencies.
The frequency of interest passing through the D-region is very important; the daytime D-region attenuation is hundreds of dB [1].
[1] Sir Edward Appleton, W.R. Piggott, Ionospheric absorption measurements during a sunspot cycle, J. Atm. Terr. Phys., (5)1, 1954, 141-172, https://dx.doi.org/10.1016/0021-9169(54)90029-X.
Abraham, Saji, et al. “Ionospheric absorption on October 24, 1995 solar eclipse.” Geophysical research letters 25.15 (1998): 2945-2947. https://doi.org/10.1029/98GL01781
SolarIrradiance.py
is a LOWTRAN model accessed through Python of Solar Irradiance, as well as transmittance ground-space.
Lowtran solar irradiance model and transmittance.
Also added was the user data mode, where your own measured temperature, pressure, humidity, etc. can be used to compute the transmittance spectrum more accurately.
Lowtran transmittance using user-provided data for model inputs.
Lowtran solar up-radiance using user-provided data for model inputs.
A
radar based on Red Pitaya
may operate for calibration or even primary purposes on lowband VHF (25-50 MHz).
A primary reason one might use a radar on lowband VHF is for example measuring human, animal or vehicle targets using
global license-free frequency bands.
Such a radar can be handheld and battery powered, with collapsible antenna rods.
The humanitarian demand for a such a radar could be significant as a much more compact and economical radar than the
coffee can radars
that have become popular after I developed the first such “personal radar” in 2006.
There are global license-free frequencies in the 40-50 MHz range.
One can consider the
US license-free frequencies in the 40-50 MHz range.
For signal processing, using FFT-based methods breaks down for the narrowband radar because the frequency separation is too small.
A common PC runs out of RAM, and certainly a small embedded CPU would as well.
I tried to squeeze more frequency resolution out of the FFT, but a normal 8 GB laptop would run out of RAM, so it’s not workable on a 1 GB RAM Raspberry Pi 3.
By instead using subspace estimation methods to find beat frequency
CWsubspace.py
the narrowband, small beat frequency separation seems initially possible in simulation.
In initial measured data, the clipped ADC due to building 60 Hz pickup makes lots of extra fake beat frequencies.
While working on the hardware fix to the clipped ADC (1:9 or 1:16 transformer for Red Pitaya input), one can make use of existing slightly clipping data by implementing a digital bandpass filter around the known center frequency.
This is “fair” and standard practice, even without the clipping.
Signal subspace estimation techniques like ESPRIT and RootMUSIC find applications in radar, whether angle of arrival estimation or for close-spaced noisy sinusoid frequency estimation.
There are
bounds on subspace estimation performance
as with any signal estimation technique, but in a subset of cases involving short range, narrowband radar, subspace estimators may have significant advantages.
Subspace estimators give excellent performance vs. FFT for a given SNR, particularly with low SNR, with far lower RAM demands for a given accuracy vs. FFT.
Mitigate CPU usage by using optimized
signal subspace estimation FORTRAN code.
Use a priori on number of sinusoids stronger than desired signal expected.
Filter signal before estimation necessary as with FFT to constrain estimate to expected frequency range (e.g. automobiles not expected to drive 500 m/s)
To parallelize the work of modeling and designing analysis, we can model expected beat frequencies for CW and/or FMCW.
Simulate estimation with noisy sinusoids at the modeled beat frequencies.
CW only gives Doppler.
A model for beat frequency vs radial velocity is is
piradar/CW_doppler.py
as
discussed here:
fbeat = 2 * v * ftx/(c-v)
FMCW can give range and velocity via the FMCW beat frequency.
FMCW typically linearly sweeps the transmit frequency.
Non-linearities in the sweep broaden the target uncertainty (increase error).
A model for expected beat frequencies is in
CalcBeat.py.
FMCW beat frequencies in Hertz are defined by
f_beat= R * B/(T c)
where R is one-way range to target [meters], c is speed of light, T is chirp cadence [Hz].
B = f_max - f_min
chirp RF bandwidth
./CalcBeat.py -tm .1 -b 10e6
For point target ranges [1, 10, 50] meters you may expect beat tones [‘0.7’, ‘6.7’, ‘33.4’] Hz.
Using 10.0 MHz RF bandwidth, you may expect 15 meters range resolution by:
ΔR = c/(2 B)
If you sweep from 902..928 MHz, then B=26 MHz.
However, using a priori knowledge about the targets, super-resolution techniques can enable better resolution proportional to SNR.
These plots demonstrate a system (CW or FMCW) that results in a true fbeat = 2.0 Hz.
The estimation is performed by ESPRIT in piradar/CWsubspace.py
Noisy time domain received signal from CW radar. Dominant signal is carrier leakage as expected for non-nulled CW radar.
Frequency domain: receive signal is blurred together with CW carrier leakage--actual separation 2 Hz. FFT exhausts PC RAM without giving enough separation to resolve target beat frequency.
Frequency domain: receive signal is blurred together with CW carrier leakage--actual separation 2 Hz. Estimation error ~ 10%.
Estimation error with ESPRIT and block length 1000 gave beat frequency estimation error of order 10%.
This simulation from
piradar/FMCW_sim.grc
was conducted with GNU Radio and the standard built-in modules.
The parameters were chosen for convenience of display in plots, not for a particular scenario.
Parameter
Value
Units
sample rate
1e6
samples/sec
sweep freq
50
Hz
VCO modulation
ramp (sawtooth)
sweep bandwidth
100e3
Hz
The simulation was performed at baseband, so upon using DUC & DDC one could sweep from 40.7 to 40.8 MHz or whatever desired frequency.
FMCW 50 Hz ramp input voltage to VCO.
VCO output in time domain. Sweep starts at 0 Hz and goes to 100 kHz in baseband at 1 MHz sampling rate with 50 Hz cadence.
VCO output in frequency domain. Sweep starts at 0 Hz and goes to 100 kHz in baseband at 1 MHz sampling rate with 50 Hz cadence.
FMCW Demodulated output--pair of targets. Periodic glitches due to discontinuity at ramp reset.
The glitches are customarily eliminated in a straightforward way by dropping those samples (or skipping them) since they are inherently periodic.
FMCW Demodulated output spectrum--pair of targets. Periodic glitches due to discontinuity at ramp reset.
The receive spectrum shows two distinct targets–one with a beat frequency of about 7 kHz and another at about 4.7 kHz.
This corresponds to a range of about 210 km and 140 km – perhaps an F-region and E-region ionosphere radar return, respectively.
The parameters used in these plots were more suitable for HF (2-10 MHz) measurements of the ionosphere.
Using low-band VHF, we can sweep a few MHz in the license-free bands, and sweep at a faster cadence.
Let’s say the radar sweeps 5 MHz at a 1 kHz cadence, then a target 20 meters away returns a 667 Hz beat frequency.
When developing an embedded sensing system, one of the very first steps is making a method to retrieve raw data from the sensor and store it.
Prevention of later problems due to unforeseen data acquisition issue scan be easily obtained by obtaining and retaining the raw data stream an essential specification from day one.
This raw data stream is the basis for important registration and unit tests.
With experience and a forward model of what to expect, one quickly gains deep insight into what is always at least partially the black box of the embedded system.
Some of the factors thereby determined include
Sample rate verification: too fast or slow?
interference: bad grounds, environmental
signal dropouts: system can’t keep up with data rate, race conditions
unexpected signals (good or bad) mixing with desired signal
Based on the above:
Is analog and or digital filter needed?
Is a PCB or circuit redesign needed
Is more advanced signal processing or different embedded CPU or GPU needed?
It initially appears that reading from digital IO with
nidaqmx-python
may be synchronous, causing too long to loop over reads, yielding error messages as below.
I wonder if nidaqmx.streaming_readers is required to read asynchronously as the name implies?
nidaqmx.errors.DaqError: The application is not able to keep up with the hardware acquisition.
Increasing the buffer size, reading the data more frequently, or specifying a fixed number of samples to read instead of reading all available samples might correct the problem.
Property: DAQmx_Read_RelativeTo
Requested Value: DAQmx_Val_CurrReadPos
This code reads 3 binary lines at once via the Port interface
fromtimeimport time
importnidaqmxasndwith nd.Task() as task:
task.di_channels.add_di_chan('Dev1/port0/line0:2',
line_grouping=nd.constants.LineGrouping.CHAN_PER_LINE)
task.timing.cfg_samp_clk_timing(rate=Fs,
sample_mode=nd.constants.AcquisitionType.CONTINUOUS,
)
print('initializing',fn)
with h5py.File(fn, 'w') as f:
f['tstart'] = time() # TODO use GPS time over serial instead b = f.create_dataset("di",(3,Nframenight),dtype=bool)
sleep(3)
print('starting fire log loop to',fn)
i=0whileTrue:
tic=time()
data = task.read(number_of_samples_per_channel=Nsampread)
print(time()-tic)
The penalty in writing boolean arrays from lists to HDF5 via h5py is the conversion from list to Numpy array implicit in writing with h5py.
The order 10 ms h5py boolean list/array write time was not the issue, rather it was the synchronous reading of digital IO that was jamming up the system.
There is an officially-supported Windows-only Python NI-DAQmx module
and a community supported Linux + Windows Python NI-DAQmx module.
Both use Python ctypes to access the underlying C API and both require NI-DAQmx to be installed first.
Both are Python 3 compatible.
There are three NI-DAQmx download choices, pick one that suits your needs:
NI-DAQmx runtime (250 MB)
NI-DAQmx runtime + configuration (1.3 GB)
NI-DAQmx (1.9 GB)
Python NI-DAQmx: choose one of:
nidaqmx: official NI support, Windows only
pip install nidaqmx
PyDAQmx: support for Linux and Windows
pip install PyDAQmx
Python NI-DAQmx for Linux and Windows:
PyDAQmx
is the
community-supported
Python program for using NI-DAQ hardware from Python in Linux and Windows.
National Instruments provides Windows-onlycomplete Python API
via ctypes for NI-DAQmx.
The Python
nidaqmx_examples
are similar to those in LabVIEW, but with text based code.
