Help with DSP, correlation, match filtering, code, etc.

[viper]

I realized there are many chapters missing in my brain that I need to load.  I'm trying to do match filtering or similar to do air-sonar detection of a known short audio signal. 

I realized I have zero experience with digital filtering, especially in coding.  In short, I know what I need but hard for me to even paint the path to get there.  I am hoping for some good white papers or vid lectures that might help me understand?  I've seen countless code examples, but they make references like I should just know what N and corr are.  I don't! 

Part of my issue is though it seems match filtering is probably the right path for my needs, I've yet to find assistance in how to "weigh" the returns.  I expect multiple matched returns due to several factors, but only one should stand out. In digital filtering, it seems to get tricky because you don;t get a visual understanding of the return. 

the TX will be I2S to an audio amp, and the RX will be a mems with I2S to an MCU.  In short, I would send I2S audio code to amp for my target signal, then want to time that return on the mic.  the issues I very much expect are numerous returns, but I need to select the right one. 

In short, I need a bunch of help understanding DSP down to the coding levels, and this probably pushes beyond a simple forum post.  I need a mental download.  It almost appears many in the DSP world live in a "I do this, it just works" world. What happens when it doesn't? How do you even troubleshoot digital filtering?  Can't really see it....

[Marco]

Part of my issue is though it seems match filtering is probably the right path for my needs, I've yet to find assistance in how to "weigh" the returns.  I expect multiple matched returns due to several factors, but only one should stand out. In digital filtering, it seems to get tricky because you don;t get a visual understanding of the return.

I don't understand time reversal signal processing, but if you just send a chirp and convolve (ie. FIR filter) the received reflections with the time reversed chirp you can just make a 2D line plot of the output and look at it. How is that not visual?

[moffy]

For context this is in relation to a previous post: https://www.eevblog.com/forum/projects/zee-sofwher!-matlab-filter-designer-audio-filter-help/

[viper]

For context this is in relation to a previous post: https://www.eevblog.com/forum/projects/zee-sofwher!-matlab-filter-designer-audio-filter-help/

Agree.  I need to backstep a bit here.  I'm missing some critical concepts. I need a better understanding of DSP. 

[Martinn]

Agree.  I need to backstep a bit here.  I'm missing some critical concepts. I need a better understanding of DSP.

The two classics for your case are first
https://www.dspguide.com/
Meant for engineers, less math heavy. Fun to read! If you ever wondered how z and laplace transform works, read chapters 32 and 33 and be enlightened.

Second classic is Lyons https://www.dsprelated.com/books/6.php
Although there are PDFs, I have both as paperback, which I find more convenient.

- Martin

[Marco]

Those won't teach you why time reversal matching is supposedly close to optimal on an unknown channel though ... that requires a whole lot of math which I can't easily follow.

Unless you want to get a Phd or are a genius, you have to take some things on faith.

[viper]

Those won't teach you why time reversal matching is supposedly close to optimal on an unknown channel though ... that requires a whole lot of math which I can't easily follow.

Unless you want to get a Phd or are a genius, you have to take some things on faith.

Ha, I am trying to digest DSPguide now.  I did appreciate in the opening of that book that it's for engineers that want to get a task completed, not if you want a new career in DSP.....lol 

As it pertains to my matched filtering work, I really thought I had that figured out because I thought correlation math was time reversed, in which I thought I was smart in thinking that if you have a short audio sweep from say 2k-3k, the 2k leaves first, and is the first to hit the mic on return.  In short, it whipped a bich, so now is reversed as it left, so you would want to look for it as time reversed.  But now I see that is not the case per the book so......  my brain took a hit there. 

I'm still trying to figure out the digital 1s and 0s as they relate to the sample points.  I come from industrial automation so separate X and Y values, but it sort of sounds like the sample points are calculated for their amplitude vs time, and assigned a value.  I could be way off on that. 

And then the mentions of float vs fix point math in an MCU.  I have some data that indicates float will be the right play, but others saying stick with fixed.  More figuring I guess. 

[viper]

Agree.  I need to backstep a bit here.  I'm missing some critical concepts. I need a better understanding of DSP.

The two classics for your case are first
https://www.dspguide.com/
Meant for engineers, less math heavy. Fun to read! If you ever wondered how z and laplace transform works, read chapters 32 and 33 and be enlightened.

Second classic is Lyons https://www.dsprelated.com/books/6.php
Although there are PDFs, I have both as paperback, which I find more convenient.

- Martin

Excellent material!  After some short reading, I realized I needed that!  This is going to take a minute and several stupid questions to figure this stuff out, but I always feel way better once I have a firm understanding. 

[gf]

As it pertains to my matched filtering work, I really thought I had that figured out because I thought correlation math was time reversed, in which I thought I was smart in thinking that if you have a short audio sweep from say 2k-3k, the 2k leaves first, and is the first to hit the mic on return.  In short, it whipped a bich, so now is reversed as it left, so you would want to look for it as time reversed.  But now I see that is not the case per the book so......  my brain took a hit there.

The core operation of pulse compression is indeed just a cross-correlation between tx_pulse and rx_pulse in the time domain, to determine the time offsets between tx_pulse and the echoes present in rx_pulse.

Cross correlation is defined as sliding dot product, so basically there is no need to deal with any time-reversal at all. Just two nested loops and a multiplication and addition in the inner loop. But computational complexity is O(N²).

Time reversal comes into play if you want to optimize computational complexity. The key insight is that for a linear, time-invariant system, convolution with a time-reversed and conjugated signal is mathematically equivalent to cross-correlation with the original signal. This allows you to perform cross-correlation using the standard convolution operation, which can be done very efficiently with the FFT. [ And you also would not calculate time reversal explicitly, but instead calculate the conjugate complex in the frequency domain. ]

What you calculate in the end is:

cross_corr = real(ifft(fft(zero_pad(rx_pulse)).*conj(fft(zero_pad(tx_pulse)))))

Those won't teach you why time reversal matching is supposedly close to optimal on an unknown channel though

Who says it's optimal on an unknown channel? (speaker + mic + amps)
The frequency and phase response of the channel certainly does affect the resulting cross-correlation.
IMO, if there is too much distortion, it may be necessary to measure and equalize the channel response.

[viper]

Here might be a dumb question, but trying to align this in my head.  Right/wrong, cross correlation could be considered passing the input audio signal past the known recorded sample points, but is executed point by point?  I am trying to understand how the match is "scored" so to speak, as there has to be some sort of uncertainty here? 

In my mind, I imagine 4 fingers from left hand passing over 4 fingers from right.  the score as the fingers overload one by one might be 1-2-3-4-3-2-1.  4 has the best alignment. 

But when I revolve the time component considering how audio is returned wave by wave, it would seem each point has to be checked individually.  It just seems to make sense that a sample or data point is just an amplitude and time component?  But it simply does not make enough sense if the match filter is to behave like a key in a lock where it's only a hit when all data points simultaneously align. 

As for the samples or N, I'm uncertain if Nyquist principles really apply if you are not trying to generate audio, just identify and time the return? 

[radiolistener]

The core operation of pulse compression is indeed just a cross-correlation between tx_pulse and rx_pulse in the time domain, to determine the time offsets between tx_pulse and the echoes present in rx_pulse.

Cross correlation is defined as sliding dot product, so basically there is no need to deal with any time-reversal at all. Just two nested loops and a multiplication and addition in the inner loop. But computational complexity is O(N²).

Time reversal comes into play if you want to optimize computational complexity. The key insight is that for a linear, time-invariant system, convolution with a time-reversed and conjugated signal is mathematically equivalent to cross-correlation with the original signal. This allows you to perform cross-correlation using the standard convolution operation, which can be done very efficiently with the FFT. [ And you also would not calculate time reversal explicitly, but instead calculate the conjugate complex in the frequency domain. ]

What you calculate in the end is:

cross_corr = real(ifft(fft(zero_pad(rx_pulse)).*conj(fft(zero_pad(tx_pulse)))))

I've tried to implement a 1D sonar simulation in Octave, check it out:

Code: [Select]

pkg load signal;


% --- Parameters ---
Fs = 48000;                  % Sampling rate (Hz)
c = 343;                     % Speed of sound (m/s)
tx_duration = 0.1;           % Duration of transmitted chirp (s)
t = 0:1/Fs:tx_duration;      % Time vector for tx signal

% --- Transmitted signal: linear chirp ---
f1 = 1000;                   % Start frequency (Hz)
f2 = 18000;                  % End frequency (Hz)
tx_pulse = chirp(t, f1, tx_duration, f2);

% --- Simulate received signal ---
% Function to simulate received echo signal
function rx = simulate_rx_pulse(tx, Fs, c, targets)
    % tx       - transmitted signal
    % Fs       - sampling rate
    % c        - propagation speed
    % targets  - Nx2 matrix: [distance(m), amplitude]
    num_targets = size(targets, 1);
    distances = targets(:, 1);
    amplitudes = targets(:, 2);
    delays = round((2 * distances / c) * Fs);
    rx_len = length(tx) + max(delays) + Fs*0.05;
    rx = zeros(1, rx_len);
    for i = 1:num_targets
        d = delays(i);
        rx(d + (1:length(tx))) += amplitudes(i) * tx;
    end
    % Add small Gaussian noise
    rx += 0.02 * randn(size(rx));
end
targets = [3, 0.3; 7, 0.6];
rx_pulse = simulate_rx_pulse(tx_pulse, Fs, c, targets);


% --- Cross-correlation using FFT (pulse compression) ---
N = 2^nextpow2(length(rx_pulse) + length(tx_pulse));
rx_fft = fft(rx_pulse, N);
tx_fft = fft(tx_pulse, N);
cross_corr = real(ifft(rx_fft .* conj(tx_fft)));

% --- Create time axis for cross-correlation ---
tau = (0:length(cross_corr)-1) / Fs;  % seconds
range = tau * c / 2;                  % convert to meters

% --- Plot results ---
figure;
plot(range, cross_corr);
xlabel('Range (m)');
ylabel('Amplitude');
title('1D Sonar');
grid on; grid minor;
xlim([0, 10]);


[viper]

The core operation of pulse compression is indeed just a cross-correlation between tx_pulse and rx_pulse in the time domain, to determine the time offsets between tx_pulse and the echoes present in rx_pulse.

Cross correlation is defined as sliding dot product, so basically there is no need to deal with any time-reversal at all. Just two nested loops and a multiplication and addition in the inner loop. But computational complexity is O(N²).

Time reversal comes into play if you want to optimize computational complexity. The key insight is that for a linear, time-invariant system, convolution with a time-reversed and conjugated signal is mathematically equivalent to cross-correlation with the original signal. This allows you to perform cross-correlation using the standard convolution operation, which can be done very efficiently with the FFT. [ And you also would not calculate time reversal explicitly, but instead calculate the conjugate complex in the frequency domain. ]

What you calculate in the end is:

cross_corr = real(ifft(fft(zero_pad(rx_pulse)).*conj(fft(zero_pad(tx_pulse)))))

I've tried to implement a 1D sonar simulation in Octave, check it out:

Code: [Select]

pkg load signal;


% --- Parameters ---
Fs = 48000;                  % Sampling rate (Hz)
c = 343;                     % Speed of sound (m/s)
tx_duration = 0.1;           % Duration of transmitted chirp (s)
t = 0:1/Fs:tx_duration;      % Time vector for tx signal

% --- Transmitted signal: linear chirp ---
f1 = 1000;                   % Start frequency (Hz)
f2 = 18000;                  % End frequency (Hz)
tx_pulse = chirp(t, f1, tx_duration, f2);

% --- Simulate received signal ---
% Function to simulate received echo signal
function rx = simulate_rx_pulse(tx, Fs, c, targets)
    % tx       - transmitted signal
    % Fs       - sampling rate
    % c        - propagation speed
    % targets  - Nx2 matrix: [distance(m), amplitude]
    num_targets = size(targets, 1);
    distances = targets(:, 1);
    amplitudes = targets(:, 2);
    delays = round((2 * distances / c) * Fs);
    rx_len = length(tx) + max(delays) + Fs*0.05;
    rx = zeros(1, rx_len);
    for i = 1:num_targets
        d = delays(i);
        rx(d + (1:length(tx))) += amplitudes(i) * tx;
    end
    % Add small Gaussian noise
    rx += 0.02 * randn(size(rx));
end
targets = [3, 0.3; 7, 0.6];
rx_pulse = simulate_rx_pulse(tx_pulse, Fs, c, targets);


% --- Cross-correlation using FFT (pulse compression) ---
N = 2^nextpow2(length(rx_pulse) + length(tx_pulse));
rx_fft = fft(rx_pulse, N);
tx_fft = fft(tx_pulse, N);
cross_corr = real(ifft(rx_fft .* conj(tx_fft)));

% --- Create time axis for cross-correlation ---
tau = (0:length(cross_corr)-1) / Fs;  % seconds
range = tau * c / 2;                  % convert to meters

% --- Plot results ---
figure;
plot(range, cross_corr);
xlabel('Range (m)');
ylabel('Amplitude');
title('1D Sonar');
grid on; grid minor;
xlim([0, 10]);


A massive THANK YOU for that!  I am now digesting that and realizing that is days of work for me right there.  Are you hobby or pro?  Also, as I have not tried this in Octave yet, did you build the parameter headers, use a template, or?  I am, however, trying to understand your TX chirp, is that a sweep?  If so, I'm trying to find that.  I see the f1 and f2. 

Also, I see to created 4 sample targets, but only two hits? 

Also, in the real world, would you also run additional filtering to improve accuracy such as DC offset remove, windowing, etc?  Regardless, I would likely do a window and mute RX most of the time unless listing for an echo.

[radiolistener]

Also, as I have not tried this in Octave yet, did you build the parameter headers, use a template, or? 

There’s no special structure here. I just like to separate different parts of the code with comments — first come the parameters that you can change, and then the code sections with corresponding comments.

I am, however, trying to understand your TX chirp, is that a sweep?  If so, I'm trying to find that.  I see the f1 and f2. 

Yes, the TX chirp is a linear frequency sweep from f1 to f2 over the duration of tx_pulse. You can also listen to it by saving it to an audio file using:

Code: [Select]

audiowrite('tx_pulse.wav', tx_pulse, Fs);.

Also, I see to created 4 sample targets, but only two hits? 

There are only two targets in the example, each defined by a pair of values — distance and echo amplitude. You can add as many targets as you like; this is just a simple simulation. I didn’t aim to exactly reproduce real-world reflections, which could include multiple echoes or multipath effects. Here, each echo is simulated as a delayed version of the transmitted pulse for simplicity. To see how it behaves under real conditions, you would need to record actual echoes.

It’s important that the recordings contain synchronized samples of both TX and RX, so the best approach is to record TX in the left channel and RX in the right channel, ensuring proper sample alignment.

Also, in the real world, would you also run additional filtering to improve accuracy such as DC offset remove, windowing, etc?  Regardless, I would likely do a window and mute RX most of the time unless listing for an echo.

I wouldn’t rush into additional processing at first. I’d suggest starting with a simple example, and then, once it works, you can try adding filters and other processing to improve accuracy.

[gf]

Yes, the TX chirp is a linear frequency sweep from f1 to f2 over the duration of tx_pulse.

BTW, the maximum tx_duration is ~35 ms to prevent the echo from the nearest object (20 ft away) from overlapping with the speaker-to-microphone leakage of the transmitted pulse. You may want to use a slightly shorter pulse, around 30 ms, to allow time for the microphone to recover from the high-power, near-field speaker-to-microphone leakage during pulse transmission before listening for echoes. However, to maximize the time-bandwidth product, still use the maximum possible pulse duration.

[viper]

Yes, the TX chirp is a linear frequency sweep from f1 to f2 over the duration of tx_pulse.

BTW, the maximum tx_duration is ~35 ms to prevent the echo from the nearest object (20 ft away) from overlapping with the speaker-to-microphone leakage of the transmitted pulse. You may want to use a slightly shorter pulse, around 30 ms, to allow time for the microphone to recover from the high-power, near-field speaker-to-microphone leakage during pulse transmission before listening for echoes. However, to maximize the time-bandwidth product, still use the maximum possible pulse duration.

Agreed!  I was running math on that and figured I may consider working with a shorter chirp only when at the min distances, but otherwise try to run a high duration.  I had also considered having the mic active during TX, which I could use the initial pulse as a timing marker for R&D.  It would be high amplitude but probably no more accurate way. 

[gf]

Yes, the TX chirp is a linear frequency sweep from f1 to f2 over the duration of tx_pulse.

BTW, the maximum tx_duration is ~35 ms to prevent the echo from the nearest object (20 ft away) from overlapping with the speaker-to-microphone leakage of the transmitted pulse. You may want to use a slightly shorter pulse, around 30 ms, to allow time for the microphone to recover from the high-power, near-field speaker-to-microphone leakage during pulse transmission before listening for echoes. However, to maximize the time-bandwidth product, still use the maximum possible pulse duration.

Agreed!  I was running math on that and figured I may consider working with a shorter chirp only when at the min distances, but otherwise try to run a high duration.  I had also considered having the mic active during TX, which I could use the initial pulse as a timing marker for R&D.  It would be high amplitude but probably no more accurate way.

If the sound pressure is very high, the speaker-to-mic leakage could even drive the mic and/or the microphone amp into saturation.

BTW, if you could keep the leakage low enough (e.g. shielding) to transmit and receive simultaneously without compromising the dynamic range, then FMCW might be an alternative to pulse radar. For FMCW you would sweep slowly over a longer measurement interval (say 5+ seconds in your case), and you would not match pulses to measure TOF, but you would measure beat frequency between TX and RX signal.

Radiolistener's example used a high sweep bandwidth of 17 kHz and a pulse with a high T*B product of 1700. This results in highly compressible pulses and a theoretical, ideal range resolution of c/(2*B) of approximately 1 cm (despite 100ms pulse duration). If your sweep bandwidth is limited to 1 kHz, the ideal range resolution cannot be better than ~17 cm. Due to the low T*B product of 30 ms * 1 kHz = 30 of the pulses, compressibility suffers, and you likely won't reach this ideal resolution. Other non-ideal factors in the signal chain will compromise the resolution, too.

EDIT: ... where range resolution ΔR refers to the ability to discriminate two echos from similar distances, R and R+ΔR. It is not the resolution of the time axis (which is just a matter of sample rate).

[viper]

[/quote]

If the sound pressure is very high, the speaker-to-mic leakage could even drive the mic and/or the microphone amp into saturation.

BTW, if you could keep the leakage low enough (e.g. shielding) to transmit and receive simultaneously without compromising the dynamic range, then FMCW might be an alternative to pulse radar. For FMCW you would sweep slowly over a longer measurement interval (say 5+ seconds in your case), and you would not match pulses to measure TOF, but you would measure beat frequency between TX and RX signal.

Radiolistener's example used a high sweep bandwidth of 17 kHz and a pulse with a high T*B product of 1700. This results in highly compressible pulses and a theoretical, ideal range resolution of c/(2*B) of approximately 1 cm (despite 100ms pulse duration). If your sweep bandwidth is limited to 1 kHz, the ideal range resolution cannot be better than ~17 cm. Due to the low T*B product of 30 ms * 1 kHz = 30 of the pulses, compressibility suffers, and you likely won't reach this ideal resolution. Other non-ideal factors in the signal chain will compromise the resolution, too.

EDIT: ... where range resolution ΔR refers to the ability to discriminate two echos from similar distances, R and R+ΔR. It is not the resolution of the time axis (which is just a matter of sample rate).
[/quote]

You hit some solid points there!  As for FMCW, I'm not sure how that would overlay on the battery power budget.  The intent with pulse was to power gate as much time as possible with many blank seconds between samples, effectively optimizing the duty cycle. 

However, you did hit on something I am struggling to digest on the "highly compressible pulses" and theoretical distance resolution.  I had been considering playing with different frequencies to help optimize.  However, it was my thought that I could optimize the resolution with a higher sample rate?  The lower the frequency, the longer the waves, so add more sample points in there due to the short chirp window?  I realized Radiolistener actually computed FFT and conj, which I believe is way more efficient than where I was headed, but I'm still trying to digest it.  At some point this nut will crack, but I'm dense....lol 

[gf]

However, you did hit on something I am struggling to digest on the "highly compressible pulses"

The attached images show the autocorrelation of a 30 ms chirp with a 1 kHz bandwidth and a 100 ms chirp with a 17 kHz bandwidth (radiolistener's example). Both have the same time scale and sample rate of 48 kSa/s. Which auto-correlation graph is "better compressed," or more compact? Which of the two do you think can better discriminate between echos with nearly the same distance, if each echo produces a shape like the one in the graph (but at different horizontal position)?

As you can see, the 100ms chirp with high bandwidth can be compressed to a much narrower auto-correlation than the 30ms chirp with low bandwidth. That's what I mean with "compressible". Not every pulse is equally well compressible.

EDIT: Ideally, you would like a pulse whose auto-correlation is a delta impulse (i.e. a single peak), with no side lobes. But in practice, there will be more or less side lobes.

[viper]

However, you did hit on something I am struggling to digest on the "highly compressible pulses"

The attached images show the autocorrelation of a 30 ms chirp with a 1 kHz bandwidth and a 100 ms chirp with a 17 kHz bandwidth (radiolistener's example). Both have the same time scale and sample rate of 48 kSa/s. Which auto-correlation graph is "better compressed," or more compact? Which of the two do you think can better discriminate between echos with nearly the same distance, if each echo produces a shape like the one in the graph (but at different horizontal position)?

As you can see, the 100ms chirp with high bandwidth can be compressed to a much narrower auto-correlation than the 30ms chirp with low bandwidth. That's what I mean with "compressible". Not every pulse is equally well compressible.

So, tell me if I am off base here but the compression thus the max theoretical resolution is more to do with what I might call "getting more hits", or more matches higher in amplitude over time because there are just more of them to match at a higher frequency? But I guess that might not explain a sweep vs single frequency.  Possibly just how the sample points uniquely land over time to get a better fix on the target?  IE, if only a single frequency, samples just land as mere copies of each other with poor match scoring.  As with a lock and key, 7pin lock, all pins are #4.  Sure a key with all #4s will match, but a range of pin sizes gets a more unique and more secure lock? 

But the lower frequency is an issue mostly with how many cycles over such a short  T length? 

Would it be practical to build the code as you have done in Octave and horse with various sweeps and times to try to optimize? 

Honestly still trying to decode the code Radiolistener provided and how that differs from the time domain X_corr direction I was headed.  Still learning....lol

[gf]

So, tell me if I am off base here but the compression thus the max theoretical resolution is more to do with what I might call "getting more hits", or more matches higher in amplitude over time because there are just more of them to match at a higher frequency?

Look at the cross-correlation of the received signal when two targets are involved, with distances R and R+18mm. The two impulses from my previously post are are still used here. The high-bandwidth chirp can distinguish separate peaks for the two echos, but with the lower-bandwidth chirp, you cannot see that two echos are present.

[viper]

My question would be, is there a way to optimize that?  My distance range is mostly being discussed at the minimum, where using a higher frequency is no problem at all, and maybe drop lower as distance increases, if really needed!  All comes down to how much attenuation or noise there is. 

But I would say, even with your R+18mm with the higher BW, if there are truly 2 echos with equal energy only 18mm apart, there is probably no way to know what is what so just average them and go.  But I would think one would stand out.  Don't know yet. 

[gf]

I've tried to implement a 1D sonar simulation in Octave, check it out:

I noticed a mistake on my side. Traditional real cross-correlation is suboptimal for chirp signals, as it leaves a carrier oscillation in the cross-correlation result. It's better to correlate the analytic signals (which results in a complex cross correlation) and to display the magnitude. The analytic signal can be easily calculated in the frequency domain by zero-ing out negative frequencies in the spectrum. Additionally, side lobes of the peaks can be mitigated by unmatched filtering with a windowed reference signal. The trade-off is a slightly wider main lobe of the peak and 1-2dB lower SNR, but I think it's worth the trade-off.

I made the following change to your script to get cleaner peaks for the targets:

Code: [Select]

% --- Cross-correlation using FFT (pulse compression) ---
N = 2^nextpow2(length(rx_pulse) + length(tx_pulse));
% window = chebwin(length(tx_pulse), 35)';
window = taylorwin(length(tx_pulse), 5, -60)';
rx_fft = fft(rx_pulse, N);
tx_fft = fft(tx_pulse .* window, N);
rx_fft(N/2+1:end) = 0;    % make analytic
tx_fft(N/2+1:end) = 0;    % make analytic
cross_corr = abs(ifft(rx_fft .* conj(tx_fft)));

EDIT: Switched from chebwin to taylorwin, as chebwin seems to produce some artifacts (numerical problems?).
Taylor window has a quite narrow main lobe, too.