Sensor Fusion - Test Bench

Experimental Setup
Huddle Test

Experimental Setup

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

Huddle Test

Load Data

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

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);

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 = 1e2*lsim(inv(G_geo), geo_1, t);
  geo_2 = lsim(inv(G_geo), geo_2, t);

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;

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]);

Dynamical Uncertainty

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

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');

5.8 KiB Raw Blame History