How to find delay of signal sample

24 visualizaciones (últimos 30 días)
Reed
Reed el 16 de Abr. de 2024
Respondida: AH el 18 de Abr. de 2024
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?
  2 comentarios
Alexander
Alexander el 16 de Abr. de 2024
Very good question. I'm interested of the answeres from the experts.
Reed
Reed el 16 de Abr. de 2024
Editada: Reed el 16 de Abr. de 2024
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.

Iniciar sesión para comentar.

Respuesta aceptada

Jon
Jon el 16 de Abr. de 2024
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')
  5 comentarios
Jon
Jon el 16 de Abr. de 2024
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
Jon el 16 de Abr. de 2024
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

Iniciar sesión para comentar.

Más respuestas (4)

Reed
Reed el 16 de Abr. de 2024
Editada: Reed el 16 de Abr. de 2024
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')

Taylor
Taylor el 16 de Abr. de 2024
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;
  2 comentarios
Jon
Jon el 16 de Abr. de 2024
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
Jon el 16 de Abr. de 2024
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

Iniciar sesión para comentar.


Alexander
Alexander el 16 de Abr. de 2024
Maybe this would be helpful:
@Manikanta Aditya used xcov but at least it is a very similar problem.

AH
AH el 18 de Abr. de 2024
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.

Productos


Versión

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by