Finding inverse cross-correlation

Hello,

I am trying to compute the inverse cross-correlation of two signals ; that is, I have two signal A and B such that A should reach its positive peak when B hits its negative peak, and I want to find the delay between the two. How can I do this in Matlab?

If it helps, I am working with respiratory data, specifically co2 and 02 volume concentrations. Our respirometer starts collecting co2 data about 1-2 seconds before it starts with the o2 - hence the delay we are trying to account for.

edit: I have enclosed a graph of some sample data to be corrected - the blue signal is oxygen concentration and the orange signal is co2 concentration. As you can see, the o2 concentration is delayed (in this case) by about 1.5 seconds, which is caused entirely by our instrumentation. The goal would be to have a function which could calculate the maximal inverse cross-correlation, so that the maxima of one signal match up with the minima of the other.

Thank you!

Answer 1

[Image Analyst]

Although often true, and maybe true for your signals, the best signal overlap is not guaranteed to occur when the cross correlation is maximum. That is a very common misconception. A simple thought experiment should reveal why. That is why I never use it and use normalized cross correlation instead. There is a function for that in the Image Processing Toolbox called normxcorr2(), however it should work for 1-D signals also I believe. I'm attaching a 2D pattern recognition demo (sorry I don't have any 1-D signal demo for it all ready to go).

[Ced]

Interesting, had for some reason never heard of it. From what I can tell, the function normxcorr2 is only made for 2D signals. I take this from the doc, and in the function, they check

validateattributes(T,{...,'2d','finite'},mfilename,'T',1)

i.e. only accept 2D matrices.

However, the function documentation is very good, even the paper used for the implementation is attached. At first glance, IF the two signals have the same size (which is, I believe, your case), I don't see the difference between this normalized crosscorrelation and the crosscorrelation of the normalized signals, but maybe I am missing something crucial.

If you open the function itself, there is actually more information:

We normalize the cross correlation to get correlation 
coefficients using the definition of Haralick and Shapiro, Volume 
II (p. 317), generalized to two-dimensions.

So, if simply normalizing these signals doesn't work, maybe you could have a look there since you don't need the "generalized 2D version".

@Image Analyst: Sorry for cluttering your comment section, but it felt better suited here.

[Image Analyst]

No problem. I too once thought that signals were best aligned when the correlation was greatest until about 30 years ago in graduate scchool when my optics professor said "That's not true. What makes you think that?" I said "I don't know - it has seemed to be true for examples I've seen." He explained that correlation is simply the sum of a bunch of products, and if the elements have high values, then the products, and hence the sum, can be high regardless if the one signal has any resemblance to the other signal at all. That's when I had my hand-slapping-the-forehead moment. Of course -- now it's so obvious. Perhaps this little demo will illustrate things better.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

% Make signal 1 100 elements long.
t1 = 1 : 100;
signal1 = exp(-(t1-20).^2/25);
subplot(2, 1, 1);
plot(t1, signal1, 'b*-', 'LineWidth', 2);
grid on;
title('Signals', 'FontSize', fontSize);

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 

% Make signal 2 the same as signal 1 but shifted 10
% and only 80 elements long
t2 = 1 : 80;
signal2 = exp(-(t2-30).^2/25);
hold on;
plot(t2, signal2, 'r*-', 'LineWidth', 2);
axTop = gca;
axTop.XTick = 0 : 10 : t1(end); % Tick marks every 10.

% Create correlation with signal2.
c2 = imfilter(signal1, signal2, 'corr', 'full');
tc = -length(signal2) + 1 : length(signal1) - 1;
% Plot it.
subplot(2, 1, 2);
plot(tc, c2, 'r*-', 'LineWidth', 2);
grid on;
title('Correlations', 'FontSize', fontSize);

% Make signal 3 the same as signal 2 but with the last 20 elements having value 1
% (It's just slightly lower so we don't have complete overlap of the curves
% and we can see them both plotted at the same time.)
signal3 = 0.9*signal2;
signal3(end-19:end) = 1;
subplot(2, 1, 1);
plot(t2, signal3, 'ms-', 'LineWidth', 2);
legend('Signal1', 'Signal2', 'signal3');

% Create correlation with signal3.
c3 = imfilter(signal1, signal3, 'corr', 'full');
% Plot it.
subplot(2, 1, 2);
hold on;
plot(tc, c3, 'ms-', 'LineWidth', 2);
grid on;
axBottom = gca;
axBottom.XTick = tc(1)-1 : 10 : tc(end); % Tick marks every 10.
legend('Correlation with signal2', 'Correlation with signal3');

helpdlg('Note how max of correlation with signal 3 is not where the Gaussians overlap.');

[Ced]

Great demo, thanks! Seems quite obvious when put like that ^^. One more reason to careful filter the result.

[S.K.]

Thank you for your suggestion, Image Analyst - I will try computing normalized cross correlation and see if it gives a result appreciably different to the maximum cross-correlation. Unfortunately, I still don't know how to compute the inverse cross correlation, normalized or otherwise ...

Answer 2

[Ced]

If I interpret your description of "inverse cross-correlation" correctly, wouldn't that just be the cross-correlation with one of the signals flipped?

[ r, lag ] = xcorr(A,-B)

[S.K.]

Hi Ced, thank you for your help!

I don't think so - the values returned by finddelay(A,B) = finddelay(A,-B), for instance, where finddelay is a matlab function which estimates delay by finding maximum cross-correlation ... I think the issue is that our signals are entirely positive, and sinusoidal about a positive axis (e.g. o2 concentration varies from 13% - 25%).

Does that make sense? I have also appended a picture of the data to the original question; hopefully that makes things clearer.

[Ced]

Ah, I just saw the figure. Yes, you should remove the mean beforehand and normalize their magnitude. Have a look at ImageAnalysts response, I was not familiar with that method. Sounds interesting.

Answer 3

[Star Strider]

I’ve not worked with respirometry data with the sort of instrumentation you describe. Does it collect the O2 and CO2 simultaneously and there is a lag because of your experimental set-up (the O2 signal lags the CO2 signal), or is the lag in the instrumentation itself? If the instrumentation analyses the signals simultaneously and lags in the output, you could just shift them. Otherwise, from your description, the easiest way do assess the lags would be to negate one signal and then just use the Signal Processing Toolbox xcorr function, for example, on both the unchanged O2 signal and negated CO2 signal. Assuming that I understand your Question correctly, that’s how I would do it.

[S.K.]

Hi Star Strider, thank you for your help!

The delay is due entirely to instrumentation - the origin of the delay is that the gas from subjects nose is pulled by a motor pump through a narrow tube to the control room and reaches the CO2 monitor and then goes to O2 monitor. The 2 monitors are connected in serial, so the same parcel of air reaches co2 first, and the other later. I don't think using xcorr(o2, -co2) would work, for the reason described in Ced's response ... I've included a bit more info in my question, which hopefully clarifies matters.

[Image Analyst]

Attaching actual data - numbers in a .mat file or something - would be also useful in case someone wants to try something with it.

[Star Strider]

The problem is that ‘respiration’ exists on the cellular level so I didn’t know if you meant cellular or ventilatory ‘respiration’.

This is an instrumentation problem going back at least a half century. Going to the literature, for example ‘Time delay technique in respiratory instrumentation’ Respiration Physiology Volume 7, Issue 3, October 1969, Pages 399–402 and PubMed (with 146 related articles, such as ‘A computer-based instrumentation system for measurement of breath-by-breath oxygen consumption and carbon dioxide production’ Biomed Sci Instrum. 1994;30:1-8) can be beneficial. They’re not all free, but you should be able to get most of them from your university library.