Sensor Fusion - Test Bench

1. Experimental Setup
2. Huddle Test

1 Experimental Setup


Accelerometer	PCB 393B05 - Vertical (link)
Geophone	Mark Product L4C - Vertical

2 Huddle Test

2.1 Load Data

load('./mat/huddle_test.mat', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 't');
dt = t(2) - t(1);

2.2 Data

acc_1 = acc_1 - mean(acc_1);
acc_2 = acc_2 - mean(acc_2);
geo_1 = geo_1 - mean(geo_1);
geo_2 = geo_2 - mean(geo_2);

2.3 Scale Data

From raw data to estimated velocity. This takes into account the sensibility of the sensor and possible integration to go from acceleration to velocity.

G0 = 1.02; % [V/(m/s2)]

G_acc = tf(G0);

T = 276;
xi = 0.5;
w = 2*pi;

G_geo = -T*s^2/(s^2 + 2*xi*w*s + w^2);

acc_1 = lsim(inv(G_acc), acc_1, t);
acc_2 = lsim(inv(G_acc), acc_2, t);
geo_1 = lsim(inv(G_geo), geo_1, t);
geo_2 = lsim(inv(G_geo), geo_2, t);

2.4 Compare Time Domain Signals

figure;
hold on;
plot(t, acc_1);
plot(t, acc_2);
plot(t, geo_1);
plot(t, geo_2);
hold off;

2.5 Compute PSD

We first define the parameters for the frequency domain analysis.

Fs = 1/dt; % [Hz]

win = hanning(ceil(1*Fs));

Then we compute the Power Spectral Density using pwelch function.

[p_acc_1, f] = pwelch(acc_1, win, [], [], Fs);
[p_acc_2, ~] = pwelch(acc_2, win, [], [], Fs);
[p_geo_1, ~] = pwelch(geo_1, win, [], [], Fs);
[p_geo_2, ~] = pwelch(geo_2, win, [], [], Fs);

figure;
hold on;
plot(f, sqrt(p_acc_1));
plot(f, sqrt(p_acc_2));
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD Accelerometers $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
xlim([1, 5000]);

figure;
hold on;
plot(f, sqrt(p_geo_1));
plot(f, sqrt(p_geo_2));
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD Geophones $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
xlim([1, 5000]);

2.6 Dynamical Uncertainty

[T_acc, ~] = tfestimate(acc_1, acc_2, win, [], [], Fs);
[T_geo, ~] = tfestimate(geo_1, geo_2, win, [], [], Fs);

2.7 Sensor Noise

[coh_acc, ~] = mscohere(acc_1, acc_2, win, [], [], Fs);
[coh_geo, ~] = mscohere(geo_1, geo_2, win, [], [], Fs);

pN_acc = p_acc_1.*(1 - coh_acc);
pN_geo = p_geo_1.*(1 - coh_geo);

figure;
hold on;
plot(f, pN_acc, '-', 'DisplayName', 'Accelerometers');
plot(f, pN_geo, '-', 'DisplayName', 'Geophones');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the Measurement Noise $\left[\frac{m/s}{\sqrt{Hz}}\right]$');
xlim([1, 5000]);
legend('location', 'northeast');