% constructs the linear system of equations using distribution derivative rule
function [A,b] = ddm_lin_sys(krnls, krlns_ders, mf_ders, sig, N)
%generic multi-frequency distribution derivative based estimator for
% non-stationary sinusoidal analysis
%
%
% [1] Michael Betser: Sinusoidal Polynomial Estimation Using The Distribution
% Derivative, in IEEE Transactions on Signal Processing, Vol.57, Nr. 12,
% December 2009
%
% krnls: matrix of all the kernels... N x R , where R is the number of 
%      non-static parameters to estimate and at the same time, the number 
%      of kernels
%
% krlns_ders: matrix of all the kernel time derivatives... N x R , where R 
%     is the number of non-static parameters to estimate and at the same 
%     time, the number of kernels
%
% mf_ders: matrix of all the model function time derivatives... N x Q , where Q 
%     is the number of model functions
%
%
% sig: vector - signal, N x 1 (CAUTION: MUST be column vector!!!)
%
% N: odd integer - signal buffer length, ...
%
% For any reasonable use, Q equals R, otherwise it makes little sense.
% Kernels must include the window function...
% 

R = size(krnls,2);
Q = size(mf_ders,2);
assert(R == size(krlns_ders, 2) );
assert(R >= Q);
% constructing the matrixes A and B from equation III.4 in [1]
% 1st dimension is the discrete time

sig_mat   = repmat(sig,   [1,R,Q]);
krnls_mat = repmat(krnls, [1,1,Q]);

mf_ders_mat = repmat(reshape(mf_ders,[N,1,Q]),[1,R,1]);

% inner product for left and right hand side of the eq.
A = shiftdim(sum(  conj(krnls_mat)  .* mf_ders_mat .* sig_mat, 1), 1);
b = shiftdim(sum(- conj(krlns_ders) .* sig_mat(:,:,1), 1), 1);

end