Scientific Computing

draft

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].

References

[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
• https://ntrs.nasa.gov/api/citations/20180001309/downloads/20180001309.pdf

SolarIrradiance.py is a LOWTRAN model accessed through Python of Solar Irradiance, as well as transmittance ground-space.

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.

If when trying to import numpy you get message:

Intel MKL FATAL ERROR: Cannot load libmkl_avx2.so or libmkl_def.so

with Anaconda/Miniconda Python, the fix is to force reinstall Numpy. It appears to be a recent glitch with Anaconda and Numpy

conda install --force-reinstall numpy

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)

Simulation

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:

f_beat = 2 * v * f_tx/(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.

Performance

These plots demonstrate a system (CW or FMCW) that results in a true f_beat = 2.0 Hz. The estimation is performed by ESPRIT in piradar/CWsubspace.py

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.

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.

The glitches are customarily eliminated in a straightforward way by dropping those samples (or skipping them) since they are inherently periodic.

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

Property: DAQmx_Read_Offset Requested Value: 0

Task Name: _unnamedTask<0>

Status Code: -200279

Problematic nidaqmx synchronous read code

This code reads 3 binary lines at once via the Port interface

from time import time
import nidaqmx as nd
with 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=0
        while True:
            tic=time()
            data = task.read(number_of_samples_per_channel=Nsampread)
            print(time()-tic)

When testing why nidaqmx-python was taking so long to read data and write to HDF5, we benchmarked h5py boolean write speed: h5py_write_speed.py.

h5py write speed benchmark

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-only complete Python API via ctypes for NI-DAQmx.

The Python nidaqmx_examples are similar to those in LabVIEW, but with text based code.

Matlab data cursors are handy for pointing out multiple data values for sharing or publication.

1. place a data cursor
2. right click, select Create New Datatip

data point

literally the value of data at a particular point.

datatip

Text annotation pointing to the data point

To move the datatip, left-click the datatip and hold the button down while moving mouse around. Put the datatip to NW,NE,SW,SE of data point

Copy plot data value to workspace by: right-click, select Export Cursor Data to Workspace

Related: Adjust Matlab data cursor precision

A first step in calibrating a whole-slide fluorescence cytometer is using a uranyl slide proving spatial uniformity to spatially normalize LED illumination and camera vignetting. We then measure the coefficient of variation (CV) for each pixel, expecting the CV to increase away from the center of slide as in the following Blue and UV excited fluorescence images.

In this data, the estimated CV is over 50 images and 1.3 ms exposure. The plots are created by PlotCV.py.

The next step is to test the system using fluorescent beads of known diameter, where the CV will be worse as a figure of merit vs. other systems.