How to find delay of signal sample

I have two signals that I am trying to align. When I perform the xcorr, the index of the peak does not match with the offset I introduced.

%Generate true signal and time
trT = 10:0.1:1500;
trPos = sin(trT*pi/200).*sin(trT*pi/50)+3*sin(trT*pi/300)+((trT-200)/200).^2;

%use part of tr data for ms
msPos = trPos(2401:5901);

%correlate signals
[r,lags] = xcorr((msPos),(trPos));
[~,idx] = max(abs(r));
figure(42)
plot(abs(r),'b.')
figure(43)
plot(trPos,'.')
hold on
plot(msPos,'.')
%zero-pad to shift the msPos to the correct position
msPos = [zeros(1,2400) msPos];
plot(msPos,'.-')
hold off
legend('True', 'Meas', 'Meas offset','location', 'northwest')

When I run the code, idx is 2520 when I would expect it to be closer to 2400. What am I missing?

[Alexander]

Very good question. I'm interested of the answeres from the experts.

[Reed]

I played around with the signals I was correlating and with simpler signals the alignment was better. I'm not sure why xcorr was producing such a poor result, though. See my workaround below.

Answer 1

[Jon]

Your problems are coming from the fact that your signal is not zero mean. The peak correlation will only correspond to the best alignment in time if the signals have no offset (zero mean). Below I have modified your code, eliminating the last term in you expression for trPos so that it is a zero mean signal. Now I think the rest of your code acts as you expect.

By the way, you can use the code button on the toolbar to insert code and have it nicely formatted.

%Generate true signal and time
trT = 10:0.1:1500;
% use zero mean signal instead (eliminate last term that isn't zero mean)
trPos = sin(trT*pi/200).*sin(trT*pi/50)+3*sin(trT*pi/300);% +((trT-200)/200).^2;
%use part of tr data for ms
msPos = trPos(2401:3401); % use shorter segment so they don't overlap
% % msPos = trPos(2401:5901);
%correlate signals
[r,lags] = xcorr((msPos),(trPos));
[~,idx] = max(abs(r));
figure(42)
plot(abs(r),'b')

figure(43)
plot(trPos,'-')
hold on
plot(msPos,'-')
%zero-pad to shift the msPos to the correct position
msPos = [zeros(1,2400) msPos];
plot(msPos,'.-')
hold off
legend('True', 'Meas', 'Meas offset','location', 'northwest')

[Alexander]

Could you please explain? The shift of 2400 is known. How can I get the index of best fit of two unknown signals with xcorr (or xconv, same problem). Thank you.

[Reed]

@Alexander You can look at the code I put together below to manually check the signal alignment. Doesn't require a zero-mean signal. Not the most elegant solution, but it works pretty well.

@Jon I think you're right. I posted my nonzero-mean workaround code below

[Jon]

Can you please be more specific, what aspect of using xcorr are you having difficulty understanding? If you haven't already, a good start would be to look at the documentation for xcorr. To see this just type doc xcorr on your command line in MATLAB

[Jon]

In case it is helpful I have written the attached function to locate a subsequence, x, in a longer sequence y. If the subsequence, x, does not occur exactly in x, the location where it has the closest fit (min sum square error) is found.

Here I have used it instead of xcorr in your original code

%Generate true signal and time
trT = 10:0.1:1500;
trPos = sin(trT*pi/200).*sin(trT*pi/50)+3*sin(trT*pi/300)+((trT-200)/200).^2;
%use part of tr data for ms
msPos = trPos(2401:5901);

% find the shift
idx = locsubseq(msPos,trPos);


figure(43)
plot(trPos,'.')
hold on
plot(msPos,'.')
%zero-pad to shift the msPos to the correct position
msPos = [zeros(1,idx-1) msPos];
plot(msPos,'.-')
hold off
legend('True', 'Meas', 'Meas offset','location', 'northwest')

Here's the function (also attached)

function idx = locsubseq(x,y)
%locsubseq Locates subsequence x in longer sequence y
% Given a length m vector x, and a length n vector y, with n>=m
% iStrt = locsubseq(x,y), provides index, iStrt, in y such that
% y(iStrt:iStrt + m -1) is the best fit (Euclidean distance) to y

% Get some dimensions
m = numel(x);
n = numel(y);

% Make sure x and y are column vectors
x = x(:);
y = y(:);

% zero pad x so it has same length as y
xpad = [x;zeros(n-m,1)];

% Make a toeplitz matrix where each column shifts x by one step (index)
X = toeplitz(xpad,[xpad(1) zeros(1,n-m)]);

% Find the sum square error from x to each subsequence in y
% 
e = (X - y).*X; % .*X multiplies unwanted values by zero
sse = sum(e.^2);

% Find index that provides best match
[~,idx] = min(sse);

[Jon]

Note this function also "slides" the subsequence by the main one looking for the best fit (as you did), it just does it all with matrices instead of loops, in case that is of interest

Answer 2

[Reed]

I solved my problem by avoiding xcorr and just "sliding" the measured signal past the true signal and minimizing the error.

%Generate true signal and time

tStep = 0.1;

trT = 10:tStep:1500;

trPos = sin(trT*pi/200).*sin(trT*pi/50)+3*sin(trT*pi/200)+((trT-200)/200).^2;

%use part of tr data for ms

tShift = 20;

msPos = trPos((2101+tShift/tStep):(6501+tShift/tStep));

msT = trT(2101:6501);

% %correlate signals

% [r,lags] = xcorr((msPos),(trPos));

% [~,idx] = max(abs(r));

%

% figure(43)

% plot(abs(r),'b.')

%Manual alignment

msError = zeros(1,(length(trPos)-length(msPos)));

msMean = mean(msPos);

for ii=1:(length(trPos)-length(msPos))

msPosPadded = msMean*ones(1,length(trPos));

msPosPadded(ii:(ii+length(msPos)-1)) = msPos;

msError(ii) = sum(abs(trPos-msPosPadded));

end

[minError,minErrorIdx] = min(msError);

tOffset = msT(1)-trT(minErrorIdx)

figure(42)

plot(trT,trPos,'.-')

hold on

plot(msT,msPos,'.')

plot(msT-tOffset,msPos,'.')

hold off

legend('True', 'Meas', 'Meas corrected','location', 'northwest')

Answer 3

[Taylor]

I think it has to do with the fact that your original signal is just multiple permutations of a sine wave. If you closely inspect the original signal, you can actually see some morphological similarities in at different time points. So xcorr is determining that spot to be lcoation of maximum correlation. If you use random data as your original signal, you don't experience this issue.

% Define the signals
signal = randn(1, 1000); % Original signal
offset = 200; % Actual offset for demonstration
segmentLength = 25; % Length of the segment
segment = signal(offset:(offset+segmentLength-1)); % Extracted segment

% Compute cross-correlation to find the offset
[correlation, lags] = xcorr(signal, segment);
[~, idx] = max(correlation);
offsetEstimate = lags(idx);

figure;
plot(signal, 'b', 'DisplayName', 'Original Signal');
hold on;
plot(segment, 'r', 'DisplayName', 'Segment');
xlabel('Sample Index');
ylabel('Amplitude');
title('Original Signal and Segment');
legend;
hold off;

% Correct the segment's position
correctedPosition = (1:length(segment)) + offsetEstimate;

figure;
plot(signal, 'b', 'DisplayName', 'Original Signal');
hold on;
plot(correctedPosition, segment, 'r', 'DisplayName', 'Corrected Segment');
xlabel('Sample Index');
ylabel('Amplitude');
title('Original Signal and Corrected Segment');
legend;
hold off;

[Jon]

While it is true that the cross correlation between periodic signals will have multiple maximums, I don't think this is the key problem with the op's case. The problem is the non-zero mean signal. In your example the signals are zero mean, so that the cross correlation peak provides the best shift.

[Jon]

To be a little more precise, the problems come from having two different length sequences, that are non zero mean. Imagine a short constant positive sequence being cross correlated with a longer positive increasing sequence. The peak correlation will occur when the right edge of the short sequence is aligned with the right edge of the increasing sequence

Answer 4

[Alexander]

Maybe this would be helpful:

https://www.mathworks.com/matlabcentral/answers/2107496-how-to-find-delay-of-signal-sample

@Manikanta Aditya used xcov but at least it is a very similar problem.

Answer 5

[AH]

As an alternative solution, you may want to use the function findsignal in Signal Processing Toolbox that finds the location of the measured signal in the true signal as below

%Generate true signal and time
trT = 10:0.1:1500;
trPos = sin(trT*pi/200).*sin(trT*pi/50)+3*sin(trT*pi/300)+((trT-200)/200).^2;
%use part of tr data for ms
msPos = trPos(2401:5901);

Let's first examine how the function works

findsignal(trPos,msPos,'NormalizationLength',length(msPos))

Upon querying the beginning and end indexes as the ouput arguments, the findsignal returns

[istart,istop,dist] = findsignal(trPos,msPos,'NormalizationLength',length(msPos))

istart = 2401

istop = 5901

dist = 1.4552e-11

The dist argument indeed indicates that the found match in trPos is really close to the measured signal.