20100215

Attempt cross-correlation of signals from multiple receivers

Read about xcorr at http://www.abiscus.com/resources/Signals/matlabsignal.pdf .

40 ms of samples around tag transmission shown for band-pass filtered (200KHz +- 2KHz) signals 1 and 4 from sig11_tag150.2_J20062.csv along with their cross-correlation in this window.

foo = csvread('sig11_tag150.2_J20062.csv'); 
Fs = 1000000;  % Sampling Frequency

% Construct a 200KHz BPF 
N      = 10;      % Order
Fpass1 = 198000;  % First Passband Frequency
Fpass2 = 202000;  % Second Passband Frequency
Apass  = 1;       % Passband Ripple (dB)
h  = fdesign.bandpass('N,Fp1,Fp2,Ap', N, Fpass1, Fpass2, Apass, Fs);
Hd = design(h, 'cheby1');

samples = 670000:710000;

% Filter the four channels
out1 = filter(Hd, foo(samples,2));
out2 = filter(Hd, foo(samples,3));
out3 = filter(Hd, foo(samples,4));
out4 = filter(Hd, foo(samples,5));

% Now compare the cross-correlation

subplot(3,1,1)
plot(out1)
title('Out1')
subplot(3,1,2)
plot(out4)
title('Out4')
subplot(3,1,3)
N = xcorr(out1, out4);
R = N(1:length(out1));
Rrx = fliplr(R');
plot(Rrx)
title('Cross correlation')

This doesn't look terribly useful - they're most correlated in the first sample and become less uncorrelated after about 1.5e4 samples, which is when the blips move past each other - this indicates that the noise itself is uncorrelated.

Cross-correlation of band-pass filtered samples 1 and 4 as above, showing first 200 samples.

Ah, let's zoom right in:

plot(Rrx(1:200))
[a,b] = max(Rrx)

It looks like they're most correlated and in phase at sample 3. Compare all four samples to see how phase differs:

subplot(3,3,1)
N = xcorr(out1, out2);
R = N(1:length(out1));
Rrx = fliplr(R');
[a, maxes(1,2)] = max(Rrx)
subplot(3,3,1)
plot(Rrx(1:20))
title('out1 & out2')

N = xcorr(out1, out3);
R = N(1:length(out1));
Rrx = fliplr(R');
[a, maxes(1,3)] = max(Rrx)
subplot(3,3,2)
plot(Rrx(1:20))
title('out1 & out3')

N = xcorr(out1, out4);
R = N(1:length(out1));
Rrx = fliplr(R');
[a, maxes(1,4)] = max(Rrx)
subplot(3,3,3)
plot(Rrx(1:20))
title('out1 & out4')

N = xcorr(out2, out3);
R = N(1:length(out1));
Rrx = fliplr(R');
[a, maxes(2,3)] = max(Rrx)
subplot(3,3,5)
plot(Rrx(1:20))
title('out2 & out3')

N = xcorr(out2, out4);
R = N(1:length(out1));
Rrx = fliplr(R');
[a, maxes(2,4)] = max(Rrx)
subplot(3,3,6)
plot(Rrx(1:20))
title('out2 & out4')

N = xcorr(out3, out4);
R = N(1:length(out1));
Rrx = fliplr(R');
[a, maxes(3,4)] = max(Rrx)
subplot(3,3,9)
plot(Rrx(1:20))
title('out3 & out4')

First 20 samples of the cross-correlation of all pairs of band-pass filtered samples from sig11_tag150.2_J20062.csv in a 40 ms window around transmission.

Which gives:

     1     2     3     4
     0    81    72     3   1
     0     0    93   125   2
     0     0     0    33   3

Does a 125 sample difference between signal 1 and 2 make any sense?

t = 125 * 1/Fs % 0.125 ms
t * c

= 37500 metres

Not really! However finding the first maximum (always occurred in the first five samples at Fs=1e6) should give the phase difference. In this sample, out1 and out3 seemed to be in phase (a local maximum at sample 1), out1 and out2 were one sample out of phase (to the left) and out2 and out3 were one sample out of phase (to the right).

Compare that visually to the original samples for verification:

figure
subplot(4,1,1)
plot(out1(15000:15020))
title('out1')
subplot(4,1,2)
plot(out2(15000:15020))
title('out2')
subplot(4,1,3)
plot(out3(15000:15020))
title('out3')
subplot(4,1,4)
plot(out4(15000:15020))
title('out4')

Four band-pass filtered signals from sig11_tag150.2_J20062.csv showing samples 685000-685020, which occurred during the CW transmission, and detailed view of phase differences.

This matches! out1 and out3 are in phase, while out2's first peak occurs one sample earlier. out4 is two samples ahead (or three behind) out1 and out3, which agrees with the xcorr.

This does make it look like resolution is a bit poor for accurately determining phase differences - with this signal at 200 KHz and 1 Msps, the wavelength is 5 samples. To locate accurately we would probably want bearings of better than 5 degrees. What will the relationship between signal phase differences and bearing accuracy be?

Perhaps the filter is mangling the shape of the signals? The leading and trailing edges of the 'blip' have a bit of a ramp. Try a lower order filter! Maybe a Butterworth one.

N      = 4;      % Order
Fpass1 = 198000;  % First Passband Frequency
Fpass2 = 202000;  % Second Passband Frequency
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fpass1, Fpass2, Fs);
Hd = design(h, 'butter');
butt1 = filter(Hd, foo(samples,2));
butt2 = filter(Hd, foo(samples,3));
butt3 = filter(Hd, foo(samples,4));
butt4 = filter(Hd, foo(samples,5));

Which gives

maxes =

     0    76    72     3
     0     0    93   125
     0     0     0    33

The filter shape doesn't seem to change much (visually), other than attenuation, which is worse (maximum of 50% less) in this Butterworth filter.

Curiously, only one figure is different - 76 vs 81, but this is a difference of one wavelength. I don't think think this particular data is very useful for anything.

The document I read had an example for computing phase differences programmatically:

Prx=csd(out2,out3);	% Estimate the Cross Spectrum
pha=angle(Prx);	% Get the phase
phase=unwrap(pha); % Unwrap the phase
plot(phase);

Try to work out what this means.