Sunday, August 5, 2018

Detection of k=7 r=1/2 convolutional codes

We start with a random binary sequence xi which is encoded using a k=7 r=1/2 convolutional decoder into
  yI = [u0,v0,u1,v1,...]  ,
with
  ui = ajxi+j and vi = bjxi+j  ,
where "+" denotes XOR and repeated indices are summed over from 1 to 7.

The property used to find a,b is the following:
  bjuj+i = ajvj+i  .   (*)

For this we start with the full state space of 127×127 values of a,b. For each pair (ui, vi) the property (*) is evaluated and states which fail it are pruned from the search space.

The screen shot below shows the output of a simple octave script demonstrating this approach where further details of this algorithm can be found. For an encoded random bit sequence, after consuming about 28 encoded bits yI, (i.e., 14 iterations), the correct a,b values are found.

Output of signal-analysis/m/test_conv.m


This method must be already known and I would be grateful if someone showed me a reference to it.

Sunday, July 15, 2018

TDoA SuperDARN/Finland

Today signals from the SuperDARN radar in Finland were active around 12375 kHz and the most likely position matches quite well with the known CUTLASS radar position near Hankasalmi.

TDoA maps


TDoA cross correlations


There are pulse trains spaced by integer multiples of 2400 μsec (0,9,12,20,22,26,27), as descibed here:

SuperDARN pulse train

30 of these pulse trains are repeated, approximately each 95.3 msec, then there is a gap of two frames and some delay, and then the next 30 pulse trains follow:

SuperDARN pulse frames

The transmissions seem to repeat each minute with a small gap at the end of the minute.

Another way to verify the used pulse sequence is to compute the autocorrelation of each line of the plot above. All differences between all combinations of the individual pulses show up as peaks in the autocorrelation:
frame number vs. autocorrelation lag

Friday, July 13, 2018

KiwiSDR TDoA (LORAN)

Below is an example of using the new KiwiSDR TDoA extension on LORAN signals, see The bias on the most likely position is quite small. Note that despite the fact than one sample at 12kHz is 83.3μs long, the RMS values of the cross correlations are ≤ 6μs. This is because the peaks in the cross correlations are fitted with a 2nd order polynomial.

Fig 1

Fig 2


Fig 3


Fig 4


Fig 5


Fig 6


Fig 7

Wednesday, July 11, 2018

Propagation delay maps based on IRI2016 electron density profiles

In order to understand better the propagation delay maps used for this blog post, the animation below was created. References for the used method are arXiv:1104.2248 and this well-known book.

Up to 2-hop propagation is considered. In case there is more than one mode of propagation, the smallest non-ground-wave propagation delay is shown in the maps. If there is no ionospheric reflection, the ground-wave delay is shown.
  • hop=0: no ionospheric reflection -> ground-wave propagation
  • hop=1: one reflection
  • hop=2: two reflections 
Propagation delays based on IRI2016 electron density profiles

Propagation delays based on IRI2016 electron density profiles

Saturday, July 7, 2018

KiwiSDR TDoA extension

KiwiSDR v1.196 includes a KiwiSDR extension (based on https://github.com/hcab14/TDoA)

Note that ionospheric delays are not taken into account, so depending on the geometry of receivers and the target the TDoA maps will be biased.

Many thanks to John Seamons, ZL/KF6VO.

In parallel I am working on a real-time version of TDoA using python/matplotlib/cartopy.

Thursday, June 21, 2018

Interesting signal (TDMA?)

Today a signal on 5099 kHz USB (5100.8 kHz center) was picked up on several KiwiSDRs in the Northwest of Europe.

This signal consists of a sequences of sharp pulses and of tones:



Zoom around the tone:



Aligning the signal in frequency to the beginning of the tone, it becomes visible that there are two tones in three intervals. It might be a coincidence but the length of the intervals is close to 6/128 seconds which may or may not be relate to the information (here or here):

It is interesting that some residual modulation is visible, i.e., the phase is not perfectly constant or linearly rising.

The next three plots show abs(z) assuming different frame lengths:
  1. 854 samples (9.11/128 sec):



  2. 844 samples (9/128 sec):



  3. 281 samples (3/128 sec):



Monday, June 11, 2018

TDoA with propagation delays based on IRI2016 electron densities

This is a follow-up of this blog post where the method referred to there is applied to TDoA maps. The first example is a STANAG 4285 signal (FUO Toulon) centered on 8436.4 kHz (USB: 8432.6 kHz)

TDoA maps using ground-wave propagation delays


TDoA using propagation delays based on IRI2016


The maps using propagation delays based on IRI2016 electron densities match better the (assumed) true position of the transmitter. It is a small effect but noticeable, and can be seen better in the side-by-side comparison below:

TDoA comparison



Another example are the TDoA maps from the last blog post where the TDoA maps based on more realistic propagation delays are more consistent between the two measurements; both indicate a position on the South coast of Cornwall:

TDoA comparison


TDoA comparison


Sunday, June 10, 2018

TDoA code update

The TDoA code now does not explicitly use gnuplot any more. All plotting is done in octave. This will it make easier to extend the code (work in progress).

Recently some interesting signals have been independently observed by several people, see, e.g., tweet. These are frequency hopping signals.

Tuning a KiwiSDR to 7812 kHz two of the signal fit into its passband. Below are TDoA maps for these signals. Note that these maps are likely biased by the fact that some of the KiwiSDRs might receive the ground-wave signal and others a sky-wave signal.

TDoA maps


TDoA cross correlations


TDoA maps


TDoA cross corrlelations

Saturday, June 9, 2018

Towards ionospheric ray-tracing

For the work on HF TDoA a realistic method for the signal delay is needed. Up to now, the great-circle distance divided by the speed of light is used which of course is incredibly naive. Nevertheless it seems to work, at least to some degree.

Rather than attempting a full-fledged ray-tracing simulation one can use the theorems of Breit and Tuve and Maryn's theorem for obtaining the signal propagation delay from virtual vertical height profiles, see, e.g, arXix:1104.2248 and here.

In order to evaluate this technique let us use the oblique ionograms measured by the University of Twente web SDR (thanks to PA3FWM for making these available! and for him and G3PLX for explanations):

The plots below show the time delay vs. UTC for two days in June at two different frequencies. Overlaid are time delays obtained by
  • VOACAP: all propagation modes with probability>0.3 are shown,
  • simulations explained above using electron density profiles obtained from the international reference ionosphere model IRI2016
(For this comparison the oblique ionograms have been corrected for a time offset of ~0.3ms whose origin is not yet clear)

All plots are generated using octave and this octave binding of the IRI2016 model.

As the IRI model is a climatological model (as is VOACAP), the agreement between model calculation and data is not expected to be perfect. Only the general features of these backscatter ionograms are reproduced, such as the 2-hop E-layer reflection for lower frequencies during mid-day but, e.g., not sporatic E-layer reflections in the late afternoon/evenings.

Note also that this simple method does not take into account magneto-ionic effects which become important for lower frequencies (and at local night-times).







Wednesday, June 6, 2018

Blind detection of checksums (1)

This is a generalization of T207_test.m to arbitrary checksums, assuming that

  • there is a fixed frame length k
  • the last m bits in a frame are a checksum
  • the checksum is formed by a (possibly non-linear) combination of the first k-m bits in the frame
  • there are no collisions between checksums, i.e., that a given sequence of k-m data bits always map onto the same checksum

The test for a T207 checksum looks like this (thanks to Antonio for providing a bit stream):

T207_test.m

Suppose that we only know that the data consist of frames of 14 bits where the last two bits are a checksum but not know how the checksum is formed and where the frames start. Then a search for the presence of such frames results in the plot below using test_checksum.m

test_checksum(bits,14,2)
Indeed the right start position of the frames is recovered when searching for the minimum number of checksum collisions without assuming anything about the way the checksum is formed besides the frame length and the number and position of the checksum bits.

This method can be easily generalized for checksums which are formed using bits both from the present and from the previous frame, e.g., the GPS navigation message, see 20.3.2.5 in IS-GPS-200H.

Tuesday, May 29, 2018

TDoA northern Europe

Here is a Link-11 SLEW signal (see, e.g., here) picked up on 5 KiwiSDRs in Northern Europe and in Iceland which seems to come from the region around Andenes, Norway:
TDoA cross-correlations
TDoA maps
Thanks to all KiwiSDR owners who enabled GNSS reception!


The plot below shows abs(z) in frames of 1600/12000 seconds:
abs(IQ)

Thursday, May 24, 2018

STANAG 4285 survey

Perhaps not surprisingly, the frame markers described in this blog post were also found in STANAG 4285 the bit streams. The table below summarizes a 30 min scan of STANAG 4285 signals using the excellent KiwiSDR @T4FM. The bitstreams were obtained using this code.

6800 - 13500 kHz, 20180522T0945-1015Z, T4FM KiwiSDR, STANAG 4285 600L and 1200L.
format count remarks
(1) 14 7-bit frames delimited by F1, see here
(3) 14 7-bit frames delimited by F4, see here
CARB 5 ITA1 coded Channel Availability and Receipt Broadcast messages, see Antonio's blog
1536-bit ACF 3 see Antonio's blog
15-bit ACF 1 see Antonio's blog

Monday, May 21, 2018

TDoA STANAG 4285

This post is for people interested in ham band intruders. In this case it is a STANAG4285 signal on 5361.8 kHz. The frequency shown on the plots is the center frequency of the modulation (5363.6 kHz = 5361.8 kHz + 1.8 kHz).

For the TDoA analysis, IQ recordings with GNSS time stamps from four KiwiSDRs located in Europe were used. The cross correlations below show at least two different modes of propagation:
Cross-correlations

Only the points marked in red which were selected by hand are used for the time-difference of arrival analysis. The location of this signal likely is Frederikshavn, Denmark: (see also this UDXF post)
TDoA


Let's now have a closer look at the signals themselves: the modulation is STANAG 4285 with 600 bps (PSK after descrambling) using the long interleaver. In the figures below one can see the HF channel response in taps of the adaptive the channel equalization filter. The bit stream is encrypted and consists of 7-bit frames marked by a LFSR sequence (format F1 in this blog post).
STANAG 4285: top - descrambled symbols; bottom: decision-feedback channel equalizer taps
STANAG 4285: top - descrambled symbols; bottom: decision-feedback channel equalizer taps
STANAG 4285: top - descrambled symbols; bottom: decision-feedback channel equalizer taps
STANAG 4285: top - descrambled symbols; bottom: decision-feedback channel equalizer taps
The code used for this analysis can be found on GitHub: https://github.com/hcab14/TDoA, and https://github.com/hcab14/signal-analysis.


This is how the cross-correlations looked a few hours later and it serves to illustrate that the ionosphere is far from being a perfect mirror:
Cross-correlations

Sunday, May 13, 2018

Some interesting FSK signals found on HF

The following FSK modulated signals were picked up recently using KiwiSDRs located at DF0KL and at Newport,OR using kiwirecorder.py

All signals use a format similar, but not identical, to STANAG 5065 where frames are delimited by LFSR-generated pseudo-random sequences. These formats may also be related to the patent EP2220799A2.

KiwiSDR location Frequency* (kHz) Date UTC Shift (Hz) Baud Format
Newport 4905.0 4/2/18 08:49 850 50 (3)
Newport 4985.0 4/2/18 08:47 850 50 (1)
DF0KL 6376.3 4/2/18 08:32 850 50 (2)
DF0KL 13419.8 4/5/18 14:46 850 50 (2)
DF0KL 16123 4/2/18 10:35 850 50 (1)
Newport 7455.0 4/2/18 08:56 850 50 (1)
DF0KL 8518.8 4/2/18 08:35 850 50 (2)
*center frequency between mark and shift

The recorded IQ data streams were demodulated using a combination of octave and C++ code, see https://github.com/hcab14/signal-analysis

For testing if a given bit sequence is generated by a LFSR, a simple method was used which is described, e.g., in this thesis. This method uses Gaussian elimination over GF(2) for finding the the generating polynomial (which is not guaranteed to be minimal, unlike the result of the Berlekamp–Massey algorithm).

All signals shown in the table above consist of 7- or 21-bit long frames which are marked by LFSR-generated pseudo-random sequences:

Format (1): 7-bit frame delimited by F1
X X X X X X F1
1 2 3 4 5 6 7


Format (2): 21-bit frame delimited by F2, F3
X X X X X X X X X X X X X F2 X X X X X 1 F3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21


Format (3): 7-bit frame delimited by F4
X X X X X X F4
1 2 3 4 5 6 7

The pseudo-random sequences F1, F2, F3, F4 are generated the following polynomials:

sequence polynomial period recursion
F1 P1=1+x28+x31 231-1 F1(t) = F1(t-3)+F1(t-31) mod 2
F2 P2=1+x+x6 26-1 F2(t) = F2(t-5)+F2(t-6) mod 2
F3 P3=1+x+x7 27-1 F3(t) = F3(t-6)+F3(t-7) mod 2
F4=F1 P4=(1+x)P1 231-1 F4(t) = F4(t-1)+F4(t-3)+F4(t-4)+F4(t-31)+F4(t-32) mod 2


For format (2) the pseudo-random sequences could have been found by autocorrelation as they are short, whereas for the other formats the length of the pseudo-random sequence is 231-1, i.e., at 50/7 baud the sequences repeats after about 3480 days only.

Format (2) which can be easily identified as a peak at position 21 in the bit stream autocorrelation is also described, e.g., here and here.

Wednesday, January 24, 2018

Using an extended Kalman filter for TDoA HF geo-location

In this post the precision of geo-location using TDoA on shortwave is explored using the data shown here before.

The idea is to feed the time series of time differences into an extended Kalman filter (EKF) using two models for the propagation time from the transmitter to the receiver:

1) Simple ray-geometry for oblique propagation in the ionosphere:
  • based on arXiv:1104.2248 which uses the Breit and Tuve + Martyn theorems to convert vertical to oblique incidence
  • using vertical electron density profiles from the IRI2016 model as an input
  • group refractive index: 1/sqrt(1-fN2/f2), i.e., no magentoionic effects are taken into account
  • when integrating over the group refractive index in the vertical direction, analytic integration of 1/sqrt(f(h)) is performed in each height interval in order to avoid the spurious numerical singularity at the reflection point
  • spherical coordinates are used (rE=6371km)
2) Alternatively, great-circle distances on the WGS84 ellipsoid are used:

All receiver coordinates are assumed to be in WGS84, however in practice it is suspected that there are deviations from the true locations of up to 20 km, as the KiwiSDR positions are provided by the KiwiSDR owners and are not taken from the GPS position solutions used for time synchronization.


The plot below shows the output of the EKF
  • implementation in octave following ADA285972
  • state = [lat, lon] of the unknown position
  • measurements = time differences between pairs of KiwiSDRs
  • the measurements covariance matrix (constant) is build from the covariances of time differences (and is set to zero for time differences not sharing a common receiver)
  • state covariance matrix = diag(0.01, 0.01)2 deg2 (constant)
  • start covariance = diag(0.1, 0.1)2 deg2
  • there are about 60 measurements which are run 5 times sequentially through the EKF

Extended Kalman Filter (EKF) output; note that the range of the x(y) axis corresponds to approximately 17(67)km.

The simple ionosphere propagation model improves the match to the known position but there is still a remaining bias of about 30km North, 4km West. Clearly, more studies are needed in order to find out if this improvement is by chance or real.


Thursday, January 18, 2018

Cyprus OTHR

Yesterday I found the OTHR from Cyprus operating with 25 Hz pulse repetition rate and μ=500kHz/s sweep rate. The most likely geo-location is just north of the known location:

TDoA plots 17030 kHz 20180113T1126Z

Cross-correlations 173030 kHz 20180113T1126Z
I was thinking about an alternative way of geo-location using an extended Kalman filter, and found this reference: www.dtic.mil/dtic/tr/fulltext/u2/a285972.pdf. It is written very well and I was able to implement it quickly in octave using dfpdp from the optim toolbox for the Jacobian. More on this later.