Re: SNR calculations
Hi Fons I am detecting beeps on CW that happen about 1/second with a pulse length of 18ms. At a decimated sampling rate of 10240 this equates to 189 samples. I found that after the smoothing filter (orange line) I could calculate the orange beep length and divide by 2 to get close to blue beep length. The blue line as you guessed is after a filter and abs. Al On Mon, 1 Apr 2024, 22:21 Fons Adriaensen, wrote: > On Mon, Apr 01, 2024 at 01:50:32PM +1300, Al Grant wrote: > > > I have a block of code in my wildlife tracker that detects high/low beeps > > in a frequency. > > It's not clear what the signal (blue trace in your png) actually > represents and what exactly you are trying to do. > > - Detect any signal that exceeds the noise level by some margin ? > > In that case S/N ratio makes sense. > > - Detect a signal at a particular frequency ? > > In that case S/N isn't relevant, S/No (signal to noise density > ratio) is the relevant metric. > > So until you provide a bit more detail, just a few comments. > > I assume the signal (blue trace) is the absolute value or > square of something - it must be nonnegative for what you > do to make sense. > > > samples = signal.convolve(samples, [1]*189, 'same')/189 > > So you are convolving with a rectangular pulse in the time domain, > which becomes a sinc in the frequency domain. This is not a very > good choice. If you are convolving anyway (there may be simpler > solutions), better use a shape that transforms to a good lowpass > filter. This would be a truncated or windowed sinc in the time > domain. > > > 189 is the number of expected samples > > What do you mean by 'expected' ? Why 189 ? > > > 1. Beep length calculations are off because of the extended > > length of "high samples" > > Look at the points in your plot where the orange line intersects > the blue pulse. This happens at half the peak value of the orange > trace. A threshold at that level will give you a good estimate of > the pulse lenght. > > > 2. SNR calculations are considerably higher (and I am not sure > > which SNR is "correct") > > Which may be because S/N isn't the relevant parameter (see above). > You need to provide more details about what comes before (how > you get the blue signal). > > Ciao, > > -- > FA > > > >
SNR calculations
Hi, I have a block of code in my wildlife tracker that detects high/low beeps in a frequency. When the signal is weak the smoothing is particularly helpful at detecting real beeps from the background - but it causing me issues with calculating beep length and SNR calculations - please let me explain: This is my smoothing: # smoothing samples = signal.convolve(samples, [1]*189, 'same')/189 189 is the number of expected samples and after smoothing looks like this ( triangular orange line / blue line before smoothing ) [image: 90b66410-2d68-4d5f-bc29-0eb568eb077a.png] The high/low samples detection looks like: # Get a boolean array for all samples higher or lower than the threshold self.threshold = np.median(samples) * 1.5 # go just above noise floor low_samples = samples < self.threshold high_samples = samples >= self.threshold # Compute the rising edge and falling edges by comparing the current value to the next with # the boolean operator & (if both are true the result is true) and converting this to an index # into the current array rising_edge_idx = np.nonzero(low_samples[:-1] & np.roll(high_samples, -1)[:-1])[0] falling_edge_idx = np.nonzero(high_samples[:-1] & np.roll(low_samples, -1)[:-1])[0] if len(rising_edge_idx) > 0: self.rising_edge = rising_edge_idx[0] if len(falling_edge_idx) > 0: self.falling_edge = falling_edge_idx[0] While the smoothing (orange line) improves beep detection when the signal is weak, it introduces two new issues: 1. Beep length calculations are off because of the extended length of "high samples" and; 2. SNR calculations are considerably higher (and I am not sure which SNR is "correct") Any suggestions on how to work with this issue would be appreciated greatly. Thanks AL
Clipping - Gain too high
Good Morning, I am trying to work out when my beeps are getting clipped due to gain too high. This is in pure python (apologies if too OT here). What I have come up with so far - which doesn't seem to work: ```python def clipping(samples, rising_edge_idx, falling_edge_idx, beep_slice): rising_edge_idx = rising_edge_idx*100 falling_edge_idx = falling_edge_idx*100 clipping = np.sqrt(np.mean(samples[rising_edge_idx:falling_edge_idx])) clipping = np.sqrt(np.square(np.real(clipping))+np.square(np.imag(clipping))) print(clipping) return clipping Rising Edge/Falling edge gives me the boundaries of a beep. I think clipping occurs when magnitude is approaching 1. The above equation is maximum at 0.5 no matter how high I ramp the gain up. Any thoughts? Thanks Al -- "Beat it punk!" - Clint Eastwood
fft beep detection with overlap
Hi All, My code is setup to detect regular beeps in a real time buffered stream of samples from a SDR. For efficiency reasons I perform beep detection using fft - and if no beeps are detected no further processing is done. This works pretty well except for one edge case - where the beep is on the edge of the fft, then it gets skipped. I think I need to overlap the fft processing. My beeps are 0.018s long with a usual sample rate of 2.048e6. So this is my current code without overlapping: def process(self, samples: SamplesT): # look for the presence of a beep within the chunk and : # (1) if beep found calculate the offset # (2) if beep not found iterate the counters and move on start_time = time.time() fft_size = self.fft_size f = np.linspace(self.sample_rate/-2, self.sample_rate/2, fft_size) num_ffts = len(samples) // fft_size # // is an integer division which rounds down fft_thresh = 0.1 beep_freqs = [] for i in range(num_ffts): fft = np.abs(np.fft.fftshift(np.fft.fft(samples[i*fft_size:(i+1)*fft_size]))) / fft_size if np.max(fft) > fft_thresh: beep_freqs.append(np.linspace(self.sample_rate/-2, self.sample_rate/2, fft_size)[np.argmax(fft)]) finish_time = time.time() # if not beeps increment and exit early if len(beep_freqs) == 0: self.stateful_index += (samples.size/100) + 14 return And this is my current attempt at overlapping: def process(self, samples: SamplesT): # look for the presence of a beep within the chunk and : # (1) if beep found calculate the offset # (2) if beep not found iterate the counters and move on fft_size = self.fft_size f = np.linspace(self.sample_rate/-2, self.sample_rate/2, fft_size) size = fft_size # 8704 # step = self.sample_rate * self.beep_duration + 1000 step = 7704 samples_to_send_to_fft = [samples[i : i + size] for i in range(0, len(samples)//fft_size, step)] num_ffts = len(samples_to_send_to_fft[0]) # // is an integer division which rounds down fft_thresh = 0.1 beep_freqs = [] for i in range(num_ffts): # fft = np.abs(np.fft.fftshift(np.fft.fft(samples[i*fft_size:(i+1)*fft_size]))) / fft_size fft = np.abs(np.fft.fftshift(np.fft.fft(samples_to_send_to_fft[0][i]))) / fft_size if np.max(fft) > fft_thresh: beep_freqs.append(np.linspace(self.sample_rate/-2, self.sample_rate/2, fft_size)[np.argmax(fft)]) # if no beeps increment and exit early if len(beep_freqs) == 0: self.stateful_index += (samples.size/100) return However it doesn't work at all. I think the way the array is constructed with this overlapping is different structure? My overlapping technique is as per here : https://stackoverflow.com/questions/38163366/split-list-into-separate-but-overlapping-chunks I am very interested if anyone can help with this please? Thank you Al
npy file into URH
Hi All, It would be super handy if I could get my npy files into a format that URH can digest. Would make it much easier to validate expected results. I hope this is not too OT here. Is there a way? Thanks Al
Python Code - Suspect dropping samples
Hi All, Some of you may recall that I am trying to sample at 2.4M and count beeps between samples (usually 80 Beeps per minute, sometimes about 47 beeps per minute). 80 beeps per second equates to a beep every 1.3sec. Anyway I have saved a 10 sec numpy array to disk, then processed it in small chunks, lets say about 0.2 seconds worth of chunks. This works perfectly. The issue is when I try to do this in realtime from the RTLSDR I am getting more like 120 beeps per minute which is wrong. I suspected the sampling rate at 2.4M might be to blame so I went as low as possible to 250k and the beeps per minute got to 86 (closer). I still need to be certain if I am dropping samples, and the code that acquires data is threaded and buffered, so quite complicated. I assume I can add some time stampss to check for missing samples against expected samples? I change the sample rate variable to 250k, but I dont think I needed to change any other calculations, since freq offset stays the same? Code is here : https://github.com/bigalnz/test_fft/blob/read_from_sdr/test_fft.py Thanks Al -- "Beat it punk!" - Clint Eastwood
Wildlife tracker - python
Good Morning, I have with some generous help from others managed to get my wildlife tracker processing signals in Python. By way of background I am trying to workout the time between CW pulses on 160Mhz where the pulses are at 80 / per min (1.33/s) 48 / min (0.8/s) or 30 /min (0.5/s) and the signal strength of each pulse and for the moment output to console in near real-time. So goals: 1. Output beeps per minute (80/48/30) 2. Output sigal strength for each beep 3. Calculate beep high state duration is 0.0176s (validation from noise) 4. Make the logic resilient to handle samples that may have no beeps in them So far it outputs: BPM: 80.0 RSSI: 1.04 BPM: 80.0 RSSI: 104. I am capturing at 2.4e6 and after decimation the samp_rate is 24e3. There is a loop in the code which finds rising edge, falling edge and then when it find another rising edge, the calculation's are done for beep interval and RSSI. A rolling index keeps the calculations stateful between chunks. Problems: 1. If the chunk happens to end in the middle of a beep the RSSI calculation is wrong since the rssi needs to mean the high state of the beep and half the beep is in the previous chunk which is not kept. 2. I probably needs some optimisation and when it runs, its a bit/start stop'ish 3. Would threading help here? The way it is at the moment I am not sure it is reading in next sample while processing current sample. The full code is here : https://github.com/bigalnz/test_fft/blob/main/test_fft.py Many thanks in advance, Al
Vector Sink & Source
Hello, Experimenting with Vector sources and sinks. Why can a Vector Source not connect to a GUI Vector Sink? https://imgur.com/a/BMxmpde I get the types are different, but why cant a vector source. well, be a source of vectors? Cheers Al
Re: Rational Resampler - block executor error
Thanks Marcus. Time to put this to the test: https://imgur.com/a/hRsYAFu I need a way to test that the values I detect in java are the same values that I get in GR. I tried a Message Debug connected to the output of "Complex to Mag" but that gave me a red arrow so not allowed. Any suggestions? Thanks A On Tue, Nov 28, 2023 at 2:07 AM Marcus Müller wrote: > > Hi Al, > > On 27.11.23 06:44, Al Grant wrote: > > Thanks for the detailed reply Marcus. > > you're most welcome! > > > *M=800* > > Where is M=800 in my rational resampler? I see looking at the link to the > > picture I posted I had Decimation=100 (is M shorthand for decimation?) > > https://imgur.com/a/B2HqCKc shows "Rational Resampler" with "Decimation: > 800", not 100. > > > *PICTURE* > > To your picture: > > > > [image: image.png] > > > > But my sample only has a repeating tone on 1 frequency, but in the picture > > above with frequency on the x axis (and power on the y?), and 3 peaks I > > have circled in red, at first glance I would have said this shows a signal > > on the 3 different frequencies? > > No, it doesn't have 3 different frequencies! You missed my point there: this > is *one* > signal. The frequency axis extends to - infinity and to + infinity; there's > infinitely > many (light green) repetitions, not 3, in a discrete-time signal. (I just > couldn't draw > infinitely wide pictures!) > We just conveniently decide that, to look at the signal, it suffices to look > at the > baseband one repetition, the dark green one. > > It does suffice to describe all the information in the signal! But: it breaks > down when > you do things like decimating or interpolating, because suddenly you get the > effects of > these repetitions you "just tried very hard to forget about" :D > > > I am assuming you didn't mean what I have written above, but I > > can't reconcile the picture with the concept you are trying to convey. > > Hope the above helps! > > > > *FILTER TYPE* > > I see you have suggested a series of resampling filters, instead of 1 big > > one, is that for computing efficiency or because it wont work the way I > > expect the way I have done it? > > > > What is the difference between just a straight decimating block and > > rational resampler? > > If you choose to use a filter that cuts out the lower 1/M of the original > bandwidth, it's > the same. (assuming the rational resampler is set to "Interpolation: 1") > > > *BASEBAND* > > Since my original post I have a slightly better understanding of a baseband > > file. Correct me if I get it wrong, but for example a RTL-SDR can capture > > baseband at 2.4Mhz (i.e. spectral width). > > Sounds right! > > > To my example I am interested in > > 160.100Mhz to 160.200 with 100 channels spaced at 160.100, 160.110Mhz etc > > etc. > > So, a bandwidth of 100 · 0.01 MHz, if I get you correctly, or 1 MHz. Add a > bit "left and > right", because the analog filters aren't perfect, so maybe 1.2 MHz, just > because we can > trivially afford to be so generous. > > > So in one baseband file I can capture all 100 channels. Cool. > > Exactly! > > > For the moment I just want to focus on 1 disaster (channel) at a time, and > > am interested in getting the file into Java and doing the processing there. > > And you seem to be doing the right thing in principle (I didn't look at the > numbers being > sensible here): you're selecting a single channel, say "OK, because that > channel is only > 1/M of the overall bandwidth, I decimate by M". > > The thing I'm not understanding about your flow graph is then why the > following low-pass > filter at all? > (it's also incorrectly parameterized, as far as I can tell, because your > input sampling > rate wasn't 32 kHz (as specified in the low-pass filter) times 800 (=25.6 > MHz) (or 100, > assuming the 800 was a typo, so 100·32 kHz = 3.2 MHz), but probably a > bandwidth that your > RTL dongle actually supports.) > > > Best regards, > Marcus > -- "Beat it punk!" - Clint Eastwood
Rational Resampler - block executor error
I am still trying to learn how GR works. Coming from Java the idea of being able to do some processing there interests me. So I am trying to use a baseband file from SDR++ as a file source, and process it in such a way that I can get the amplitude in Java. I presume this would mean reading in the bin file and decoding the bits to the I and Q values. The source file is an unmodulated pulse on about 160.7807Mhz about 2 times per second. Here is my block setup: https://imgur.com/a/B2HqCKc And a link to the github project with the baseband file : https://github.com/bigalnz/java_test The issue I am currently having : block_executor :error: sched: (1)> is requesting more input data than we can provide. ninput_items_required = 27862 max_possible_items_available = 8191 If this is a filter, consider reducing the number of taps. What is going on here and how do I fix this? Thanks Al
Wildlife Tracker (again)
Background:
I am trying to get GR to receive a signal from a wild life transmitter
on about 160.7073Mhz and calculate the beeps per minute - with an
ultimate goal of logging it to a db or csv.
The project is not for profit and I work with an endangered species
which are fitted with these transmitters. Because the birds are in a
very remote location we would like to use a Rpi4+ and RTL-SDR to
record daily the status of the birds as we only get to the remote
location occasionally ( a few times a year).
I may even invest in a SDRPlay which I understand is a better build of
receiver. I have some java programming experience, and python looks
easy enough for this...anyway to detail what I have done so far:
Signal:
The transmitted pulse is an unmodulated burst of RF energy about 0.2s
in duration and the time between pulses can vary to indicate the
animal is alive (30 beeps per minute), dead (90bpm) or incubating
(48bpm).
Every 10 minutes the are a series of other pulses which convery even
more information, but for now I just want to determine the time
between pulses and workout beeps per minute from that.
In my research I see that a very similar question has already been
asked by others both on this list and here
https://dsp.stackexchange.com/questions/50103/checking-for-vhf-pulse-with-sdr-and-python
with great answers from Andy Walls. I have read all of that and:
1. Read the beginner tutorials
2. Read
https://dsp.stackexchange.com/questions/50103/checking-for-vhf-pulse-with-sdr-and-python
3. With help from another GR user got a baseband recording of the
signal isolated in a GR project
4. Have completed the tags tutorial on the wiki :
https://wiki.gnuradio.org/index.php?title=Python_Block_Tags
5. Tried to modify the tags tutorial as it applies to my project
6. Got as far as : https://github.com/bigalnz/tracker
This is where I got stuck - I tried to create a embedded python block
using the tags wiki article to my pulse isolator.
I expected `for index in range(len(input_items[0])):` to result in a
counter from 0,1,2,3,.x but it just seemed to be 0 and 1?
The second issue I encountered was I thought I could detect peaks with
this logic:
if input_items[0][index-1] > input_items[0][index] and
input_items[0][index-1] > input_items[0][index-2]:
But that did not put the tags on the peaks - so before I can even
think about calculating time gaps I presume I need to get the peak
detection working, and to do that I need to understand why the index
is not iterating beyond 0 and 1.
I am very grateful to those that have helped so far, and appreciate
any further input I can get.
Being able to capture this data will be really helpful for the future
of this species - full disclosure this is not a project I am marked
on, I am not a student and am not making any money from this.
Much thanks in advance
Al
** CODE SNIPPET **
def work(self, input_items, output_items):
print(f"length of input_items : {len(input_items[0])}")
for index in range(len(input_items[0])):
#print(f" the index is increasing : {index}")
if (input_items[0][index] >= self.threshold and
self.readyForTag == True):
#print("inside loop *")
#print(f"index : {index} input itmes index :
{input_items[0][index]}")
#print(f"index : {index} input itmes -1 :
{input_items[0][index-1]}")
#print(f"index : {index} input itmes -2:
{input_items[0][index-2]}")
if input_items[0][index-1] > input_items[0][index] and
input_items[0][index-1] > input_items[0][index-2]:
#print(f"peak at {input_items[0][index-1]}")
key = pmt.intern("detect")
value =
pmt.from_float(np.round(float(input_items[0][index-1]),2))
writeIndex = self.nitems_written(0) + index-1
#print(f"Index:{writeIndex}")
#print(f"Value:{value}")
self.add_item_tag(0, writeIndex, key, value )
self.readyForTag = False
if (self.readyForTag == False):
self.timer = self.timer + 1
if (self.timer >= self.report_period):
self.timer = 0
self.readyForTag = True
output_items[0][:] = input_items[0]
return len(output_items[0])
