xcorr and finddelay giving erroneous results for offset sine waves

19 views (last 30 days)
Hi all.
I am using the finddelay function to verify a fractional delay FIR. The signal is then upsampled before going through a finddelay function to return a time delay which I then verify against.
Unfortunately, finddelay appears to return an incorrect result. I ran an implementation of the xcorr function in parallel which returns the same result.
In this example, my testbench is set up to delay the signal by 156.25 samples; both verification methods return a 126 sample delay. I have verified that my FIR delays the signal correctly visually in Figure 2.
Both functions window the signals in their entirety by default.
Relevant function call below, followed by plots of the arguments used for finddelay.
delay_out_samp = finddelay(real(sig_chopped(1,:)),real(sig_chopped(2,:)));
Figure 1 - Plots of sig_chopped(1,:)) and real(sig_chopped(2,:))
Figure 2 - Zoomed in Figure 1 visually verifying correct time delay (i.e., 11648-11491 = 157 samples)
I cannot see any reason that these methods would fail. Both tones seem to be suitably conditioned (e.g., equal gain, no zero-padding, etc.). Frustrating for sure. Any help/guidance is much appreciated. I can provide data sets publicly, script privately. Thanks in advance.
Some additional context and observations
  • Tone in this example is 20MHz sampled at 0.3215Gsps for 100 samples
  • Upsample rate is 500
  • Requested time delay enacted by the FIR equates to 1ns
  • Verification method claims a delay of 0.8064ns, an offset of 0.1936ns
  • This offset decreases as the requested time delay decreases, though remains in the same order of magnitude.
  • Using a higher tone frequency or sampling rate reduces this systematic error
  • Changing the upsample rate only changes the precision of the verification method. The systematic error still remains
  4 Comments
James Draper
James Draper on 27 Sep 2022
Edited: James Draper on 27 Sep 2022
Hi Sharmin. Thank you for the info. I brough xcorr out of finddelay hoping the lag data would hold some clues.
Maxlag in the above case will have been 33999 when calling finddelay with no maxlags argument. I've tried using 10, 500, 10000, but these all give poor results.
The signal has been upsampled 500 times and finddelay indeed returns an integer: 126 (this should've been 156 or 157).
Thanks for the pointer to alignsignals. Unfortunately, this gave the same erroneous result.

Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 28 Sep 2022
Edited: David Goodmanson on 28 Sep 2022
Hi James,
This effect is due to a peculiarity in finding delays when your sample is cut off sharply at each end. The following code approximates your signals:
t = 0:3.4e4;
ncyc = 4.35 % number of cycles
delay = 156
tmax = t(end);
w0 = 2*pi*ncyc/tmax;
y1 = sin(w0*(t+delay));
y2 = sin(w0*t);
figure(1)
plot(t,y1,t,y2)
grid on
Here the time is in one second intervals, so time and array index are basically the same thing.
eps = finddelay(y1,y2)
eps = 126 % the problem
If the signals were not sharply cut off but had many cycles contained in a slowly varying gaussian envelope for example, I don't think there would be a problem. As it is, it looks like finddelay is finding the max absolute value of xcorr and returning that lag value. xcorr is looking at the discrete version of
Integral sin(w0*(t+delay-eps))*sin(w0*t) dt
and finding the max as a function of eps. Defining
offset = delay-eps
then it's a case of finding the maximum of
Integral sin(w0*(t+offset))*sin(w0*t) dt
as a function of offset. In your case,
offset = 30
It turns out that if the window has an even number of quarter-cycles, the offset is 0 and there is a good result, i.e. finddelay outputs the actual delay. If you have an odd number of quarter-cycles, there is a maximum possible value for the offset and a maximally bad result. Unfortunately, your signal has 4.35 cycles which is not far from 17 quarter-cycles. With some algebra you can show that
offset1 = (1/w0)*atan((1-cos(2*w0*tmax))/(2*w0*tmax-sin(2*w0*tmax)));
offset1 = 29.2742
so it agrees. (there may be off-by-one issues that do not seem worth chasing down). Consequently if w0 is known, you can determine the actual delta with
delta = finddelay(y1,y2) + offset1
It's interesting to change ncyc from 4 to 4.1 to 4.25 to ... 4.5 and watch the offset change from good to bad and back again.
  2 Comments
David Goodmanson
David Goodmanson on 30 Sep 2022
Hi James,
Thanks for pointing out the file exchange entry. Because of its very limited applicability I don't consider my answer to be a practical method, but rather an explanation of an odd and unexpected (to me anyway) effect.

Sign in to comment.

More Answers (1)

Mathieu NOE
Mathieu NOE on 28 Sep 2022
see code demo below for zero crossing (with linear interpolation for best time accuracy)
you have access to both positive slope zc points (red diamonds) as well as negative slop zc points (green diamonds)
also my demo uses coarse sine wave just to demo that the zc points are indeed obtained by interpolation and not simply picked from the nearest point;
% dummy data
% create two sine waves with 0.15 s delay
n = 20;
x= 11*(0:n-1)/n;
y1 = sin(x-0.1);
y2 = sin(x-0.25);
threshold = 0; % "zero" crossing or any other value here
[t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(y1,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
% ind => time index (samples)
% t0 => corresponding time (x) values
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
[t0_pos2,s0_pos2,t0_neg2,s0_neg2]= crossing_V7(y2,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
figure(1)
plot(x,y1,'b-*',t0_pos1,s0_pos1,'dr',t0_neg1,s0_neg1,'dg','linewidth',2,'markersize',12);grid on
hold on
plot(x,y2,'k-*',t0_pos2,s0_pos2,'dr',t0_neg2,s0_neg2,'dg','linewidth',2,'markersize',12);grid on
hold off
% time delta on positive slope zero crossing points
dt_pos = t0_pos2 - t0_pos1
% time delta on negative slope zero crossing points
dt_neg = t0_neg2 - t0_neg1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% CROSSING find the crossings of a given level of a signal
% ind = CROSSING(S) returns an index vector ind, the signal
% S crosses zero at ind or at between ind and ind+1
% [ind,t0] = CROSSING(S,t) additionally returns a time
% vector t0 of the zero crossings of the signal S. The crossing
% times are linearly interpolated between the given times t
% [ind,t0] = CROSSING(S,t,level) returns the crossings of the
% given level instead of the zero crossings
% ind = CROSSING(S,[],level) as above but without time interpolation
% [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
% par = {'none'|'linear'}.
% With interpolation turned off (par = 'none') this function always
% returns the value left of the zero (the data point thats nearest
% to the zero AND smaller than the zero crossing).
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
% if no time vector is given, use the index vector as time
t = 1:length(S);
elseif length(t) ~= length(S)
% if S and t are not of the same length, throw an error
error('t and S must be of identical length!');
end
% check the level input
if nargin < 3
% set standard value 0, if level is not given
level = 0;
end
% check interpolation method input
if nargin < 4
imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S = S - level;
% first look for exact zeros
ind0 = find( S == 0 );
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind);
s0 = S(ind);
if ~isempty(ind)
if strcmp(imeth,'linear')
% linear interpolation of crossing
for ii=1:length(t0)
%if abs(S(ind(ii))) >= eps(S(ind(ii))) % MATLAB V7 et +
if abs(S(ind(ii))) >= eps*abs(S(ind(ii))) % MATLAB V6 et - EPS * ABS(X)
% interpolate only when data point is not already zero
NUM = (t(ind(ii)+1) - t(ind(ii)));
DEN = (S(ind(ii)+1) - S(ind(ii)));
slope = NUM / DEN;
slope_sign(ii) = sign(slope);
t0(ii) = t0(ii) - S(ind(ii)) * slope;
s0(ii) = level;
end
end
end
% extract the positive slope crossing points
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);
% extract the negative slope crossing points
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);
else
% empty output
ind_pos = [];
t0_pos = [];
s0_pos = [];
ind_neg = [];
t0_neg = [];
s0_neg = [];
end
end
  2 Comments

Sign in to comment.

Categories

Find more on Get Started with Signal Processing Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!