Estimation of the attenuation of two waves on a linear sensor array

Question

Context and objective

I am trying to estimate the attenuation coefficients ($\alpha_u$ and $\alpha_v$) of two waves ($\overrightarrow{u}$ and $\overrightarrow{v}$). These waves propagate on a linear antenna of $N$ sensors, noted hereafter $\overrightarrow{s}_{i=1, ..., N}$, equally spaced at a distance $\Delta_x$. For each sensor $M$ samples are available: $\overrightarrow{s}_{i}=(s_{i,1}, ..., s_{i,M})^T$. The two waves as well as the sensors are placed on the same line and:

• The $\overrightarrow{u}$ wave is emitted by a source located at a distance $\Delta_u$ from the first sensor $\overrightarrow{s}_{1}$
• The $\overrightarrow{v}$ wave is emitted by a source located at a distance $\Delta_v$ from the last sensor $\overrightarrow{s}_{N}$

The figure below summarizes the situation:

Hypotheses and model

• The sensors record only the $\overrightarrow{u}$ and $\overrightarrow{v}$ waves (no noise) => $\overrightarrow{s}_i = \overrightarrow{u}_i + \overrightarrow{v}_i$
• The $\overrightarrow{u}$ wave propagates with a linear attenuation coefficient $\alpha_u$ and with a constant celerity $c_u$ => $\overrightarrow{u}_i=e^{-\alpha_u d_{u,i}}\overrightarrow{u}(t-d_{u,i}/c_u) ; d_{u,i}=\Delta_u + (i-1)\Delta_x $
• The $\overrightarrow{v}$ wave propagates with a linear attenuation coefficient $\alpha_v$ and with a constant celerity $c_v$ => $\overrightarrow{v}_i=e^{-\alpha_v d_{v,i}}\overrightarrow{v}(t-d_{v,i}/c_v) ; d_{v,i}=\Delta_v + (N-i)\Delta_x $

Question

How to estimate the linear attenuation coefficients $\alpha_u$ and $\alpha_v$ from the $N$ sensor records $\overrightarrow{s}_{i=1, ..., N}$, knowing $\Delta_x$, $\Delta_u$, $\Delta_v$, $c_u$, $c_v$ but not $\overrightarrow{u}$ and $\overrightarrow{v}$?

Attempts

Energetic approach:
My first attempt was to use an energetic approach to get rid of the problems due to the propagation time of the waves. Let's start with the simple case where we only have one source, the solution is then trivial.

With a single source ($\overrightarrow{v}=\overrightarrow{0}$) and for $M$ sufficiently large we have:
$E_{s_1}\equiv \sum_{k=1}^{M}s_{i,k}^2$
$E_{s_2}\approx E_{s_1}e^{-2\alpha_u\Delta_x}$
and therefore:
$\hat{\alpha_u}=-\frac{ln(\frac{E_{s_2}}{E_{s_1}})}{2\Delta_x}$

Case with two sources:
$E_{s_i}\equiv \sum_{k=1}^{M}s_{i,k}^2=\sum_{k=1}^{M}(u_{i,k}+v_{i,k})^2=E_{u_i}+E_{v_i}$ (Assuming $\overrightarrow{u}$ and $\overrightarrow{v}$ to be independent) $E_{u_i}=E_{u}e^{-2\alpha_u d_{u,i}}$ with $E_{u}$ the energy of $\overrightarrow{u}$
$E_{v_i}=E_{v}e^{-2\alpha_v d_{v,i}}$ with $E_{v}$ the energy of $\overrightarrow{v}$

This approach leads to a system of $N$ non-linear equations with $4$ unknowns ($α_u$, $α_v$, $E_u$, $E_v$). What bothers me with this approach is:

• The difficulty to solve the system when $N$ is small
• The absence of analytical solutions
• The introduction and the determination of two unknowns $E_u$, $E_v$ of which are not useful.

Matlab code simulating the scenario

N = 10 ;
M = 100000 ;
Delta_x = 10 ;
Delta_t = 0.001 ;
Delta_u = 2 ;
Delta_v = 5 ;
c_u = 50 ;
c_v = 40 ;

alpha_u = 0.0111 ;
alpha_v = 0.0123 ;

u = randn(M,1) ;
v = randn(M,1) ;
% u = zeros(M,1) ; u(round(M/3),1) = 1 ;
% v = zeros(M,1) ; v(round(M/2),1) = 1 ;

s = zeros(M,N) ;
for ii = 1 : N
    d_u_ii = Delta_u + (ii-1)*Delta_x ;
    d_v_ii = Delta_v + (N-ii)*Delta_x ;

    tau_u = round(d_u_ii/c_u/Delta_t) ;
    tau_v = round(d_v_ii/c_v/Delta_t) ; 

    u_ii = exp(-alpha_u*d_u_ii).*[u(end-tau_u+1 : end) ; u(1:end-tau_u)] ;
    v_ii = exp(-alpha_v*d_v_ii).*[v(end-tau_v+1 : end) ; v(1:end-tau_v)] ;

    s(:,ii) = u_ii + v_ii ;
end

Answer 1

Since both waves are uncorrelated, cross correlation should do the trick here. If you calculate the cross correlation between any two antenna signals you will end up with two pronounced peaks, one at positive lags and one at negative lags. The location of the lags will give you the velocity and direction of each waves

The peak of the cross correlation will be proportional to the square of the amplitude, so you use the ratios of any two different peaks to estimate the attenuation.

The example below shows this by only using subsequent frames. There are various ways to improve the accuracy of the estimate. For example You can certainly use cross correlation between all frames and you use interactively use estimates to subtract wave 1 from wave 2 and vice versa.

EDIT: Outline for wave subtraction

Once you have an initial estimate you can do the following

1. Create a beamformer to estimate the original signals for both waves.
2. Using this signal estimate you can estimate the individual contributions of both waves at each antenna.
3. Subtract these estimates $v$ from the antenna and redo the estimation for $u$ and vice versa
4. Now you have updated estimates for attenuation and velocity. Go back to step 1 and repeat until you have reached the desired accuracy or the estimates stop moving.

Code

%% solve by cross correlation of consecutive frames
xx = zeros(2*M-1,N-1);
xmax1 = zeros(N-1,2);
xmax2 = zeros(N-1,2);
for i = 1:N-1
  xx(:,i) = xcorr(s(:,i),s(:,i+1));
  % find the peak for tau > 0 and tau < 0, find both value and location of
  % peak
  [xmax1(i,1), xmax1(i,2)] = max(xx(1:M,i));
  [xmax2(i,1), xmax2(i,2)] = max(xx(M+1:end,i));
    
end
% correct time axis 
xmax1(:,2) = xmax1(:,2) - M;
%% Estimate attenuation
% the amplitude of subsequent frames will rise/fall with e(=alpha*Delta_x)
% and hence the peak of the autocorrelation should go with the square of
% this. We can simply divide the max of subsequent frames, take the log,
% differentiate, divide by two and take the mean

est_alpha_1 = abs(0.5*mean(diff(log(xmax1(:,1))))/Delta_x);
est_alpha_2 = abs(0.5*mean(diff(log(xmax2(:,1))))/Delta_x);

% velocity is simply dX/dT
est_c_1 = Delta_x/(mean(xmax1(:,2)*Delta_t));
est_c_2 = Delta_x/(mean(xmax2(:,2)*Delta_t));

fprintf('Attenuation Estimate Wave 1 = %8.6f\n', est_alpha_1);
fprintf('Attenuation Estimate Wave 2 = %8.6f\n', est_alpha_2);
fprintf('Speed Wave 1 = %8.2f\n', est_c_1);
fprintf('Speed Wave 2 = %8.2f\n', est_c_2);