view ddm_der2_lin_sys_sol.m @ 4:72c011ed1977 tip

more elaborate example with non-stat. estimate explanation

function [dddlss ddlss dlss A_der2_plus A_der_plus A_plus] = ...
  ddm_der2_lin_sys_sol(A_sys,b_sys,A_der_sys,b_der_sys,A_der2_sys,b_der2_sys,N_fft,R,Q)
dlss = zeros(Q+1, N_fft-R+1);
ddlss = zeros(Q+1, N_fft-R+1);
dddlss = zeros(Q+1, N_fft-R+1);

for l=1:N_fft-R+1
  A = A_sys(:,:,l);
  b = b_sys(:,:,l);

  AA_inv = inv(A'*A);
  A_plus = AA_inv*A';
  ddm_est_  = fliplr((A_plus * b).');
  dlss(:,l) = [ddm_est_ 0];
  
  A_der = A_der_sys(:,:,l);
  b_der = b_der_sys(:,:,l);    
  AA_mix = ( A_der' * A + A' * A_der );
  A_der_plus = (-AA_inv * AA_mix * AA_inv * A' + AA_inv * A_der');
  ddm_est_der_ = fliplr((A_der_plus *b + A_plus * b_der).');
  ddlss(:,l) = [ddm_est_der_ 0];

  A_der2 = A_der2_sys(:,:,l);
  b_der2 = b_der2_sys(:,:,l);  
  
  AA_inv_der  = - AA_inv * AA_mix * AA_inv;
  AA_mix_der  = A_der2'*A + 2*A_der'*A_der + A'*A_der2;
  A_der2_plus = ...
    - AA_inv_der * AA_mix     * AA_inv     * A'...
    - AA_inv     * AA_mix_der * AA_inv     * A'...
    - AA_inv     * AA_mix     * AA_inv_der * A'...
    - AA_inv     * AA_mix     * AA_inv     * A_der'...
    + AA_inv_der * A_der' + AA_inv * A_der2';
  ddm_est_der2_ = fliplr((A_der2_plus *b +2 * A_der_plus * b_der + A_plus*b_der2).');
  dddlss(:,l) = [ddm_est_der2_ 0];

end
end