Search code examples
matlabnoisenoise-reductiontransformation-matrix

How to find transformation matrix from the output with Gaussian noise?


For the below given input and output, the matrix A can be found out by pseudoinverse or mrdivision in MATLAB. Similarly, I would now like to know, how to determine A, if my output signal Y matrix contains additive zero mean, uncorrelated, Gaussian noise?

x1 = [1 1 1]';
x2 = [0 1 1]';
x3 = [0 0 1]';
x4 = [1 0 1]';

y1 = [1 2 0]';
y2 = [-1 0 3]';
y3 = [3 1 1]';
y4 = [5 3 -2]';

X = [x1 x2 x3 x4];
Y = [y1 y2 y3 y4];

A = Y/X

Also, I have modelled the unknown noisy output as below:

y1_n = y1 + sqrt(var(y1))*randn(size(y1));
y2_n = y2 + sqrt(var(y2))*randn(size(y2));
y3_n = y3 + sqrt(var(y3))*randn(size(y3));
y4_n = y4 + sqrt(var(y4))*randn(size(y4));
Y = [y1_n y2_n y3_n y4_n];

Solution

  • The statement A = Y/X solves the linear system of equations A*X = Y. If the system is overdetermined, as in your case, the solution given is the least squares solution. Thus, if you have additive, zero mean, uncorrelated, Gaussian noise on Y, then A = Y/X will give you the best possible, unbiased, estimate of A.

    Note that the noise you add to your Y matrix is quite large, hence the estimate of A is far away from the ideal. If you add less noise, the estimate will be closer:

    x1 = [1 1 1]';
    x2 = [0 1 1]';
    x3 = [0 0 1]';
    x4 = [1 0 1]';
    X = [x1 x2 x3 x4];
    
    y1 = [1 2 0]';
    y2 = [-1 0 3]';
    y3 = [3 1 1]';
    y4 = [5 3 -2]';
    Y = [y1 y2 y3 y4];
    
    for n = [1,0.1,0.01,0]
       Y_n = Y + n*randn(size(Y));
       A = Y_n/X;
       fprintf('n = %f, A = \n',n)
       disp(A)
    end
    

    Output:

    n = 1.000000, A = 
        2.9728   -5.5407    2.8011
        2.6563   -1.3166    0.6596
       -3.3366    1.1349    1.5342
    
    n = 0.100000, A = 
        2.0011   -4.0256    2.9402
        1.9223   -1.0029    1.0921
       -3.1383    1.9874    1.0913
    
    n = 0.010000, A = 
        1.9903   -3.9912    2.9987
        1.9941   -1.0001    1.0108
       -3.0015    2.0001    1.0032
    
    n = 0.000000, A = 
        2.0000   -4.0000    3.0000
        2.0000   -1.0000    1.0000
       -3.0000    2.0000    1.0000
    

    Of course if you make X and Y larger by adding more vectors you'll get a better estimate too, and will be able to compensate more noisy data.