4
$\begingroup$

I am currently implementing an acoustic communication system in Matlab. For synchronization purposes, every transmission contains an upchirp preamble. The receiver uses cross-correlation to find the preamble in a captured audio file and decodes the transmitted data afterwards. This works well in case the recording actually contains a signal. However, if the audio file only contains noise, cross-correlation locates the preamble somewhere in the noise where no signal exists. Decoding would be performed based on random samples in this case.

To avoid this, I first want to determine whether the recording contains a signal at all and only in the second step determine the exact offset of the signal in the recording. I believe that the solution I am looking for is similar to the mechanism Amazon Echo uses to detect the trigger word.

I am not sure how to implement this, however, I've been thinking about two approaches. For the first one, let's consider the following code snippet:

% Sample rate and time vector
Fs = 48000;
t = 0:1/Fs:1-1/Fs;

% Signal of interest: Chirp
c = chirp(t, 18000, 1, 20000);

% Generate noisy signal that contains the SOI
s = [zeros(1, 5000) c zeros(1, 5000)];
s = awgn(s, 8);

% Generate pure noise sequence
n = wgn(1, 10000 + numel(c), 0);

% Cross-correlation of the SOI with the generated signals
[r, lag] = xcorr(s, c);
[r2, lag2] = xcorr(n, c);

% Estimate offset of the SOI in the generated signals
[~, I] = max(abs(r));
chirp_offset = abs(lag(I)) + 1;

[~, I2] = max(abs(r2));
chirp_offset2 = abs(lag2(I2)) + 1;

% Plot cross-correlation
subplot(211)
plot(lag, r);
title("Cross-correlation of chirp with signal that actually contains the chirp");
subplot(212)
plot(lag2, r2);
title("Cross-correlation of chirp with pure noise");

% Output estimated offset of the SOI
disp("Estimated SOI offset in samples in the first signal (expected: 5001): " + chirp_offset);
disp("Estimated SOI offset in samples in pure noise: " + chirp_offset2);

The code generates two noisy signals, one of which contains a pre-generated upchirp. Using cross-correlation, the offset of the upchirp in both signals is estimated. The code generates the following plots: cross-correlation plots In the upper plot, we see a single high peak indicating the start of the upchirp in the signal. As the second signal only contains noise, the lower plot does not have a clear peak.

Is it possible to decide based on the result of the cross-correlation whether a recording contains the signal of interest? My approach to this would have been to search the output of the cross-correlation for a peak whose magnitude is above a certain threshold. If such a peak exists, I would assume that the signal is present.

My second idea was to measure the signal power in the frequency band where I expect the signal of interest, and compare that to the power in a frequency band where I expect only noise. If the power in the "interesting" band exceeds the noise frequency band by a certain amount, I assume that the signal is present. Again, one would have to set a decision threshold in this case.

Is any of the above suggestions a valid approach to solving my problem? If yes, how do I determine an appropriate threshold in each case? If no, is there a better solution? In general, I would prefer to have a slightly higher risk of false positives than to miss a transmission.

Any help would be highly appreciated. Thanks.

$\endgroup$
8
  • $\begingroup$ I like both these approaches, have you ran tests to determine the percentage of false-positives / false-negatives? $\endgroup$
    – Jdip
    Aug 25, 2022 at 0:39
  • $\begingroup$ The simplest solution would be to threshold the maximum value of the cross correlation (in your example 30000 vs. 600). However, this works assuming your signal is somehow normalized. $\endgroup$
    – Florian
    Aug 25, 2022 at 7:25
  • $\begingroup$ @Jdip: I haven't performed any major experiments yet, as I'm still researching how signal detection is usually implemented. $\endgroup$
    – F105
    Aug 25, 2022 at 15:33
  • $\begingroup$ @Florian: Assuming the signal is normalized, would you just hardcode the threshold? $\endgroup$
    – F105
    Aug 25, 2022 at 15:35
  • $\begingroup$ Isn't this equivalent to a matched filter and the usual tools of receiver operating characteristic (ROC), etc apply? $\endgroup$
    – D Duck
    Aug 26, 2022 at 7:27

3 Answers 3

6
$\begingroup$

Is it possible to decide based on the result of the cross-correlation whether a recording contains the signal of interest?

Certainly. The easiest way would be look at the crest factor (peak to RMS ratio) of the cross correlation and set a threshold.

Another must have is a checksum on the pre-amble and the actual data. This ensure that what you decode is actually the intended data.

Another standard practice is the number the packets. This way you can tell if you have dropped one or receiving data out of sequence.

One final note: acoustic modems are very tricky since the channel and the environment are highly variable. This includes distance between microphone and speaker, directiv, non-linear distortions in transducers and/or amplifiers, room acoustics and reverb time, strong individual reflections and dozens of types of acoustic background noises (traffic, people, air conditioner, etc)

It's important to nail down the requirements up front: the achievable data rates can easily vary over a factor of 1000 or more between a good and bad (but realistic) scenario and different environments may require different algorithms altogether.

$\endgroup$
3
  • $\begingroup$ I imagine you could also search for just the first half of the preamble, and if it's not followed by the second half, it's spurious so ignore it. $\endgroup$
    – user253751
    Aug 25, 2022 at 13:47
  • $\begingroup$ @Hilmar: Thanks for pointing out the crest factor, I didn't know about that. Would you determine an appropriate threshold through experimentation or is there a standard practice on how to set it (if certain parameters are known)? At the moment, my implementation is intended to be a proof of concept, assuming mainly white Gaussian background noise. $\endgroup$
    – F105
    Aug 25, 2022 at 15:49
  • $\begingroup$ this will depend a lot on the quality (or lack thereof) of your channel. If you put a check sum on the pre-amble, you can set it quite low. "False positives" are not a problem: If you find a peak, decode the preamble and verify the checksum. If it doesn't check out just try the next peak. I don't know how "realistic" your proof of concept is supposed to be, but "gaussian white noise" doesn't really happen in the acoustic world. $\endgroup$
    – Hilmar
    Aug 26, 2022 at 13:55
2
$\begingroup$

Have you tried setting a threshold based on the signal statistics?
For your first approach (could be applied to your second as well, only in the frequency domain), that would look something like:

threshold = mean(corr) + alpha * std(corr)

with some appropriately chosen alpha (I don't have access to your data and don't know what type of noise we're talking about here, but you would want to choose alpha based on the noise statistics).

$\endgroup$
2
  • $\begingroup$ I actually thought about something similar. But do you mean mean(corr) + alpha * std(corr) or mean(abs(corr)) + alpha * std(corr)? Without taking the absolute value, the mean of the correlation is very small. Regarding noise statistics: My implementation is mainly intended to be a proof of concept. Currently, I assume white Gaussian background noise. I assume that I have to determine the value for alpha heuristically, right? Based on the code above, I would have to set alpha to at least 7, but I suppose this depends on the actual power of the background noise. $\endgroup$
    – F105
    Aug 25, 2022 at 15:29
  • 1
    $\begingroup$ Yep, you need to experiment and see what works best for your application ;) Not sure about abs(). Again, try and see which makes more sense! $\endgroup$
    – Jdip
    Aug 25, 2022 at 15:58
0
$\begingroup$

1.- Example low frequency up-down chirp pulse

Fs=1e3; % [Hz] sampling frequency
dt=1/Fs % [s] base time step

% example chirp(t1,f1,t2,f2) 
% t1: start pulse t2: stop pulse, both t1 t2 scalars
t1 = 0; % [s] time reference for single pulse
t2= 1;

t=[t1:dt:t2];

f1=0 % [Hz] start freq
f2=25 % [Hz] stop freq
p0 = chirp(t,f1,t2,f2);  % up-chirp puse

figure(1)
plot(t,p0)
xlabel('t')
grid on
title(['chirp single pulse f1=' num2str(f1) 'Hz f2=' num2str(f2) 'Hz'])

enter image description here

A chirp complete pulse of interest for your system, a single pulse, may well be double duration than the trace generated with command chirp, because to avoid the abrupt change in frequency / amplitude / phase that most probably causes one way freq chirp pulses (up-chirp down-chirp), on could for instance do 1st half of the pulse pull up from f1 to f2, and on the 2nd half of the pulse bring down from f2 to f1.

up-down-chirp

p01 = chirp(t,f1,t2,f2);  
p10 = chirp(t,f2,t2,f1);

if to avoid abrupt 180 amplitude change right in the middle of the pulse

if p01(end-1)==p10(1)
    p010=[p01([1:end-1]) p10];
else
    p010=[p01([1:end-1]) -p10];
end
t2=[t([1:end-1]) t(end-1)+t];

figure(2)
plot(t2,p010)
xlabel('t')
grid on
title(['chirp single pulse f1=' num2str(f1) 'Hz f2=' num2str(f2) 'Hz'])

enter image description here

One could also truncate the output of chirp, and use as actual chirp pulse just a portion of the trace produced by command chirp.

With f1=100Hz f2=250Hz, just audio frequencies, the plot window already gets all-blue, with the curve, so when running chirp with the parameters supplied in the question, evenmore there's little point attempting to visualize the generated pulse

2.- Pulse as in the question, back-to-back

Fs = 48e3; % [Hz] sampling
dt=1/Fs % [s]

t1=0; % [s] time reference for single pulse
t2=1;

t=[t1:dt:t2];

% PULSE
p0 = chirp(t,f1,t2,f2);  % up-chirp puse

% Pulse generation : chirp
f1=18e3 %18e3 % [Hz] chirp start freq
f2=2e4 % [Hz] chirp stop freq

p01 = chirp(t,f1,t2,f2); % pulse to detect
p10 = chirp(t,f2,t2,f1);

% if to avoid abrupt 180 amplitude change right in the middle of the pulse
if p01(end-1)==p10(1)
    p010=[p01([1:end-1]) p10];
else
    p010=[p01([1:end-1]) -p10];
end

t2=[t([1:end-1]) t(end-1)+t];

% figure  % meaningless all-blue plot
% plot(t2,p010)
% grid on

p010(end)=[];t2(end)=[]; % forcing cs1=2

np=numel(p010); % amount samples pulse

Now a clean pulse for your system is contained in p010 with time reference t2 .

3.- Signals

q1=randi([1e5 5e5],1,1)
q2=randi([1e4 5e5],1,1)
p1 = [zeros(1,q1) p010 zeros(1,q2)]; % pulse placed in loose t

if mod(numel(p1),2)~=2 p1(end)=[]; end  % forcing cs1=1

% q1+q2+np  % check

t=dt*[0:1:numel(p1)-1]; % building time reference

figure(3)
ax3=gca
stem([t(1) t(end)],[1 1],'Color','c'); % time stamps t(1) t(end)
grid on
xlabel('t');
hold(ax3,'on')
stem([t(q1+1) t(q1+np)],[1 1],'LineWidth',2,'Color','b'); % time stamps pulse start stop
stem(t(q1+floor(.5*(np))),1,'Marker','^','LineWidth',2,'Color','r'); % time stamp pulse centre
grid on
xlabel('t');
title('signal, no noise')

% NOISE
varn1=.1;   % small noise variance: 10% signal amplitude
varn2=1.2   % large noise variance : 110% signal amplitude
n1 = randn(1,numel(t))*varn1;
n2 = randn(1,numel(t))*varn2;

% SIGNAL + NOISE(var=.1)
s1=p1+n1;

plot(ax3,t,s1)
grid on
xlabel('t');
title(ax3,['Signal + Noise(var = ' num2str(varn1) ')'])

enter image description here

while with small noise one can still appreciate pulse presence

% SIGNAL + NOISE(var=1.2)
s2=p1+n2;

figure(4)
ax4=gca
plot(ax4,t,s2)
hold(ax4,'on')
grid on
xlabel('t');
stem(ax4,[t(1) t(end)],[1 1],'Color','c'); % time stamps t(1) t(end)
stem(ax4,[t(q1+1) t(q1+np)],[1 1],'LineWidth',2,'Color','b'); % time stamps pulse start stop
stem(ax4,t(q1+floor(.5*(np))),1,'Marker','^','LineWidth',2,'Color','r'); % time stamp pulse centre

title(ax4,['Signal + Noise(var = ' num2str(varn2) ')'])

enter image description here

with large noise the pulse goes unseen, now it is below noise level, yet still detectable.

The overlapped arrow shows where the pulse centre position in time is, where the correlation/convolution should point at.

4.- Detection

% Correlation 
% [r1, nlag1] = xcorr(s1,p1);
% [r2, nlag2] = xcorr(s2,p1);

I can repeat it with xcorr but I had it already with conv so showing solution with convolution

Regardless of any even/odd signal length the following applies :

length(conv(x1,x2))=length(x1)+length(x2)-1

% Convolution low noise signal
r1 = conv(s1,p010);

% Convolution high noise signal
r2 = conv(s2,p010);

xcorr and conv both increase resulting signal length.

Matching correlation/convolution length to initial signal length can be done in many ways.

The following if-clause and switch work ok as long as not trivial cases.

r1 r2 can be shortened to s1 s2 lengths because pulses neither start half way nor end up halved.

Directly removing pulse half length from beginning of r1 r2, and another half length from end of r1 r2

if mod(numel(r1),2)==0 && mod(np,2)==0 cs1=1; end
if (mod(numel(r1),2)==0 && mod(np,2)~=0) ||...
   (mod(numel(r1),2)~=0 && mod(np,2)==0) cs1=2; end
if mod(numel(r1),2)~=0 && mod(np,2)~=0 cs1=3; end

switch cs1
    case 1  % even-even
        r1([1:np/2-1 end-np/2+1:end])=[];
    case 2  % even-odd or odd-even
         r1([1:floor(np/2)-1 end-floor(np/2)+1:end])=[];
    case 3  % odd-odd
         r1([1:floor(np/2)-1 end-floor(np/2):end])=[];
    otherwise
end

if mod(numel(r2),2)==0 && mod(np,2)==0 cs2=1; end
if (mod(numel(r2),2)==0 && mod(np,2)~=0) ||...
   (mod(numel(r2),2)~=0 && mod(np,2)==0) cs2=2; end
if mod(numel(r2),2)~=0 && mod(np,2)~=0 cs2=3; end

switch cs2
    case 1  % even-even
        r2([1:np/2-1 end-np/2+1:end])=[];
    case 2  % even-odd or odd-even
         r2([1:np/2-1 end-np/2+1:end])=[];
    case 3  % odd-odd
         r2([1:np/2-1 end-np/2:end])=[];
    otherwise
end

max(r2) is expected time stamp to look for sync

[maxval,maxind]=max(r2)
td=t(maxind)

this detection timestamp is already aligned to s1 s2 time references, not to the correlation/convolution wider time references.

stem(ax3,t(maxind),-1.5,'Marker','v','LineWidth',2,'Color','g')

stem(ax4,t(maxind),-6.5,'Marker','v','LineWidth',2,'Color','g')

enter image description here

enter image description here

5.- Overlapping clean signal

plot(ax3,t,p1,'y')
stem(ax3,t(maxind),-1.5,'Marker','v','LineWidth',2,'Color','g')

enter image description here

plot(ax4,t,p1,'y')
stem(ax4,t(maxind),-6.5,'Marker','v','LineWidth',2,'Color','g')

enter image description here

Detection works this way ok.

however, MATLAB should be translated to a low level language so that detection of for instance voice or/and video and RF/uw signals could be performed in real time.

Note that this pulse could stand even larger noise.

Such strong and reliable correlation regardless of strong noise presence is the base for code-based detection like CDMA and Spread-Spectrum modulations.

You may be able to implement the system you are after with a couple Raspberri-Pi boards.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.