% script for computing "active noise canceller"

%% locations
source = [1 5];
speaker1 = [1 0];
speaker2 = [-1 0];
mic1 = [0.25 0];
mic2 = [-0.25 0];

%% signal parameters
A_source = 10;
f_source = 100;
phi_source = 0;

v_sound = 340;

%% source signal as observed by microphone
d1 = sqrt(sum((source-mic1).^2));
d2 = norm(source-mic2); 

tau1 = d1/v_sound;
tau2 = d2/v_sound;

phis1 = -2*pi*f_source*tau1;
phis2 = -2*pi*f_source*tau2;

%% phase shifts due to delay from speakers to mics
d11 = norm(speaker1-mic1);
d12 = norm(speaker1-mic2);
d21 = norm(speaker2-mic1);
d22 = norm(speaker2-mic2);

phi11 = -2*pi*f_source*d11/v_sound;
phi12 = -2*pi*f_source*d12/v_sound;
phi21 = -2*pi*f_source*d21/v_sound;
phi22 = -2*pi*f_source*d22/v_sound;

%% system of equations for cancellation:
% 
D = [ exp(j*phi11) exp(j*phi21)
      exp(j*phi12) exp(j*phi22)];
y = -A_source*exp(j*phi_source)*...
    [exp(j*phis1)
     exp(j*phis2)];

 x = D\y; % instead of inv(D)*y
 
 %% convert to polar coordinates
 disp('Phasors for transmitted signals:')
 A = abs(x)
 phi = angle(x) %*180/pi - don't convert to degrees, otherwise checks don't work
 
%% checks:
% both of the following expressions should be zero
c1 = A_source*exp(j*phi_source)*exp(j*phis1) + ...
    A(1)*exp(j*phi(1))*exp(j*phi11) + ...
    A(2)*exp(j*phi(2))*exp(j*phi21);
c2 = A_source*exp(j*phi_source)*exp(j*phis2) + ...
    A(1)*exp(j*phi(1))*exp(j*phi12) + ...
    A(2)*exp(j*phi(2))*exp(j*phi22);


%% phasor at origin
% figure out phases from delays to origin (0,0)
phis0 = -2*pi*f_source*norm(source)/v_sound;
phi10 = -2*pi*f_source*norm(speaker1)/v_sound;
phi20 = -2*pi*f_source*norm(speaker2)/v_sound;

% Compute phaasors at origin
phasor0 = A_source*exp(j*phi_source)*exp(j*phis0) + ...
    A(1)*exp(j*phi(1))*exp(j*phi10) + ...
    A(2)*exp(j*phi(2))*exp(j*phi20);

% convert to polar representation
disp('Phasor at origin:')
A0 = abs(phasor0)
phi0 = angle(phasor0)
