How to find delay of signal sample
조회 수: 27 (최근 30일)
이전 댓글 표시
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
채택된 답변
Jon
2024년 4월 16일
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
Jon
2024년 4월 16일
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
2024년 4월 16일
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
추가 답변 (4개)
Taylor
2024년 4월 16일
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
Jon
2024년 4월 16일
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
2024년 4월 16일
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
Alexander
2024년 4월 16일
Maybe this would be helpful:
@Manikanta Aditya used xcov but at least it is a very similar problem.
댓글 수: 0
AH
2024년 4월 18일
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))
The dist argument indeed indicates that the found match in trPos is really close to the measured signal.
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!