I am a hobbyist programmer, and for "fun" are trying decode signals from a wildlife tracker to assist with a local conservation project.

With assistance from the GNU Radio community I have got as far as decoding frequency from a AirpyHF+ SDR in python very well using a phasor and FIR - but when I spawn threads using this method on a Rpi5 - I can only process about 5 channels at once with a Fs of 768k.

Each channel is 10khz apart, and each signal has bandwidth of about 1.5Khz.

It has been suggested that a polyphase filter might be more efficient for processing multiple channels, There are 100 channels in total, but I would be happy if I could process 10 at once on the Pi platform.

Despite a lot of reading, looking at githib libraries, and GNU radio pages I still could not even begin to write a juypter notebook to test how it might be implemented.

I am hoping that someone here might be kind enough to offer some assistance.

I can provide a link to a sample file from the SDR.

Thank you

Yes please post a link to your sample files. It's possible I may have some time tomorrow in which case I would be happy to help you generate an example of what you are doing, assuming having that in Python would be useful to you. 

That would be superhelpful. 

Here is a recording : https://geekhelp-my.sharepoint.com/:f:/g/personal/...

And I have included a info text file with the details but just in case 

Fc = 160425000 Mhz

Fs = 912000 Samp/s

The channels are 10Khz apart, but of course we often only detect a few channels and they are not always next to each other and drift slightly from published freq (generally only  a few khz max)

Frequencies in signal in Mhz:

[160242564, 160289923, 160338577, 160412115, 160439884, 160559311, 160601258]

I have made a Juypter notebook showing my original method (spawning FIR filters and smoothing on full bandwidth (not efficient)) and a fft channelized method.

The channelized method looks very efficent, but in the sample signal it has added a beep thats not there - the accurate timing between beeps is critical to the application.

A good example to show the difference is in my notebook here : 

https://github.com/bigalnz/test_fft/blob/20-scan-n... 

But I do not know why the channelize method is picking up the extra beep causing the beep rate to be 139 (see the end)?

I'm not an expert, but typically the channelizer incorporates a filter.

That could help - do you have an example of a FFT channelizer that incorporates a filter?

Perhaps a differrent appraoch to the question.

My original method for 1 signal with FIR looked like:

sample_rate = 768000

center_freq = 160425000

size = len(samples)

x = samples

t = np.arange(len(x))/912000

h = signal.firwin(501, 0.02, pass_zero=True)

p = np.exp(2j*np.pi*t*47000)

x = x * p

x = signal.convolve(x, h, 'same')

x = x[::100]

x = np.abs(x)

x = signal.convolve(x, [1]*189, 'same')/189

This produces beautful time domain signals:

And spectrogram (of unprocessed input):

from matplotlib import pyplot as plt

plt.specgram(samples, NFFT=1024, Fs=768000, Fc=160425000 )

plt.show()

# Shows 12 beeps

Now with the new FFT method I get 

import scipy.signal

#d = np.fromfile(open("test_multiple_v2.fc32"), dtype=np.complex64)

#d = np.fromfile(open("six_channels_768000.fc32"), dtype=np.complex64)

#d = np.fromfile(open("neighbouring_ch_912000_3.fc32"), dtype=np.complex64)

d = np.fromfile(open("../data/10s_160377933.fc32"), dtype=np.complex64)

f0 = 160.425e6    # Center frequency

Fs = 768e3        # Sampling rate

#Fs = 912e3        # Sampling rate

N_fft = 1024      # Number of FFT channels

# Generate array of channel frequencies

f = (np.fft.fftshift(np.fft.fftfreq(N_fft, 1/Fs)) + f0) / 1e6

# Time tag each sample

t = np.arange(len(d)) / Fs

# Reshape so we can do an FFT over an axis

d_fft = d.reshape((-1, N_fft))

D = np.fft.fftshift(np.fft.fft(d_fft, axis=1))

# Time tag each sample coming from a channel

T = np.arange(len(D)) / Fs * N_fft

# Now convert into power spectral density

# 1. Reshape to (N_timestep, N_int_per_timestep, N_fft)

# 2. Square

# 3. Sum over N_int_per_timestep axis

N_time_PSD = 1500

PSD = (np.abs(D.reshape((N_time_PSD, -1, N_fft)))**2).mean(axis=1)

# Create overall spectrum

spec = PSD.mean(axis=0)

# Find peaks (note: I hand-tuned prominence)

p = scipy.signal.find_peaks(spec, prominence=0.05)[0]

plt.figure()

plt.imshow(db(PSD), aspect='auto', extent=(f[0], f[-1], T[-1], T[0]))

plt.xlabel("Frequency [MHz]")

plt.ylabel("Elapsed time [s]")

plt.colorbar()

plt.show() 

Notice the extra beep about halfway?

Its the same input file - so what am I getting the extra beep?

https://i.imgur.com/NAiTdx8.png

I think the correct result has 13 beeps. The first beep occurs right at the beginning of the record, within the first second. Your time domain plot clearly shows it but you did not plot a red X over that first spike. Your spectrogram does not show that first spike. It may have gotten swallowed in the filter transient.

Hi Lito844,

Thanks for the reply. 

My transmitters should transmit beeps at a normalized rate of 80 beeps per/minute in the test signal. They vary slightly in time so allowed range is from 67-80 beeps per minute (bpm).

With the first FIR method the beeps counted are:

Beeps per minute : 78.62, 80.35, 70.88, 80.08, 69.43, 78.26, 79.99, 71.51, 79.93, 68.96, 78.96 

This code that gives me this result is tried and tested!

The fft channelize method gives me: 

Beeps per minute 
 ** note the unexpected 139.31 BPM output in the array ** 
 
[ 79.92895204  71.42857143  80.07117438  68.91271057  78.94736842
 139.31888545  69.01840491  78.94736842  79.92895204  71.42857143
  80.07117438  68.91271057]

(note 70.88 from above corresponds I think to the 71.42 below)

Why this 139?

Is this extra beep indicated by arrow in the middle (not worried about the total number)

My plot on left (samples loaded without modification)
FFT channelized plot on right

Might be clearer image here

https://imgur.com/a/W1NIdXG

(I note its present in your plot too - but to confirm that is not expected)

Here are my beeps from my FIR method:

Which is perfect and all within limits.

Here is the beeps from fft method (notice extra one in the middle below the 0 in 160.378)

Your bug is in the line:

D = np.fft.fftshift(np.fft.fft(d_fft, axis=1))

It should read:

D = np.fft.fftshift(np.fft.fft(d_fft, axis=1), axes=(1,))

The issue is that np.fft.fftshift() by default fftshifts both axes. Thus there was not "extra" beep, just misplaced beeps due to fftshift'ing on both frequency AND time axes.

Screenshot:

https://i.imgur.com/aDicODE.png

P.S. I had the same error in the image I posted earlier in this thread.

Yes - that was exactly it.

This FFT channelizer approach appears to be *very* efficient and testing on good strong signals so far it is working well. I am going to test what happens when I get signals on neighbouring freqs (10khz apart) and also weaker signals - will some smoothing after the FFTs help with weaker signals?

The method you prototyped is likely the simplest in both implementation logic and computational load. 

What you pay for this simplicity is lack of control over things like channel spacing, filter response, and output sample rate.

I'd stick with it unless your requirements or intellectual curiosity demand more. Bear in mind, depending on total throughput and filter requirements, running in real time on RPI may not be feasible with a more capable channelizer.

If you want to upgrade then you'd want to investigate the Polyphase Filter Bank channelizer (aka PFB or Polyphase Analysis Filter Bank) and/or the Weighted Overlap and Add Filter Bank (WOLA) channelizer. Of the two, the PFB is more popular although I can't tell you why.

PFB consists of a polyphase filter feeding an IFFT, with some shift register gymnastics for updating the polyphase state to maintain phase continuity. Harris [1] has a thorough analysis. If memory serves GNU Radio has a PFB implementation.

WOLA applies weights to the taps of a shift register followed by a stack and add time aliasing stage, followed by an FFT and finally multiplication with a rotating phasor. The best derivation I know of is in the Crochiere and Rabiner text [2]. The only open source WOLA I've seen is [3] although I cannot vouch for it.

Hope this helps.

[1] Harris, F.J. (2021). Multirate Signal Processing for Communication Systems (2nd ed.). River Publishers

[2] Crochiere, R.E. and Rabiner, L.R. (1983) Multirate Signal Processing. Prentice-Hall, Englewood Cliffs.

[3] https://github.com/TheNeuralBit/cyclo_channelizer

Yes you raise good points about the simplicity of the current approach, and I need to do more testing.

Putting aside the PFB approach for a minute, is there a case for smoothing the channels that have a signal. 
I feel like the result of the channelizing means I no longer have a true IQ signal - although I can treat it in a similar fashion.

With the FIR approach I was hitting a limit of about 5 channels in realtime.

With the FFT channelizer inital estimates are showing I am more likely to be able to do 40 channels at once.

In relation to the PFB - does this allow you select what how many channels, what size of channel, and what frequency to center each channel on? (versus the current approach which is to get 1024 channels of equal width without any control over Fc for each channel and then dump the ones not needed)

An FFT filter bank is capable of perfect reconstruction (PR) meaning you can take its outputs, run them through a synthesis filter bank that performs inverse processing steps, and reconstruct the original signal.

In terms of smoothing, applications of this type of short-time Fourier transform processing commonly windows each block prior to FFT. For PR there are conditions for overlapping successive blocks when windowing. I'll refer you to the Harris text or others on this forum to provide details.

On your second question, PFB implementations that I've worked with allow for controlling number of channels (M), decimation (D), filter bandwidth (BW), and possibly other filter parameters such as transition bandwidth (TBW).

M can be any positive integer, D can be any positive integer up to M. This flexibility allows the possibility for channels to be overlapping and/or differing bandwidths.

Thanks. 

I guess what I am asking is, with PFB, can you have for example 4 filters, at arbitary positions of differing bandwidths - like the picture below where the green boxes are the channels / filters (which is a made up example).

Arbitrary may be too strong a word but to be sure there's a lot of room for creativity, especially if you're not concerned with efficiency.

At the end of the day you'll choose an FFT size which will determine the maximum number of channels. The spacing between FFT bins determines the minimum channel spacing. You'll choose an input sampling frequency Fs which will set the bandwidth of the input. You'll choose a number of channels M and decimation D, the ratio of which (M/D) will determine how the output channels sampling rate relates to the input Fs. Perhaps you'll get creative with the FFT itself and figure out a way to stack bins for wider bandwidth channels.

Although all of these ideas (and more) are certainly in play, they may not have practical implementations. Further, channelizers with the level of flexibility you're asking about are well outside of the mainstream, meaning you're not likely to find an open source plug and play realization, though I suppose its possible.

Right - and can I still make some "measure" (dbfs?) of signal power when I am using this FFT channelizer approach?

Absolutely. In fact a common application is to center the channelizer channels on signal channels from the system you wish to process. When properly designed the stream of samples flowing out of a channelizer output is processed as effectively and in exactly the same manner as if it were produced by a discrete tuner/filter lineup. Examples of such processing include demod, signal detection and measurement. Your prototype did its own coarse measurement followed by threshold detection, proving the concept.

Try setting up your FFT channelizer so that its channel centers align with the channels you wish to process. Ensure that the output channel sample rate is sufficient to process your signal of interest.

Test performance by comparing accuracy to your FIR implementation using signals over a range of power levels with and without adjacent channel interference.