dan@0: % constructs the linear system of equations using distribution derivative rule
dan@0: function [A,b] = ddm_lin_sys(krnls, krlns_ders, mf_ders, sig, N)
dan@0: %generic multi-frequency distribution derivative based estimator for
dan@0: % non-stationary sinusoidal analysis
dan@0: %
dan@0: %
dan@0: % [1] Michael Betser: Sinusoidal Polynomial Estimation Using The Distribution
dan@0: % Derivative, in IEEE Transactions on Signal Processing, Vol.57, Nr. 12,
dan@0: % December 2009
dan@0: %
dan@0: % krnls: matrix of all the kernels... N x R , where R is the number of 
dan@0: %      non-static parameters to estimate and at the same time, the number 
dan@0: %      of kernels
dan@0: %
dan@0: % krlns_ders: matrix of all the kernel time derivatives... N x R , where R 
dan@0: %     is the number of non-static parameters to estimate and at the same 
dan@0: %     time, the number of kernels
dan@0: %
dan@0: % mf_ders: matrix of all the model function time derivatives... N x Q , where Q 
dan@0: %     is the number of model functions
dan@0: %
dan@0: %
dan@0: % sig: vector - signal, N x 1 (CAUTION: MUST be column vector!!!)
dan@0: %
dan@0: % N: odd integer - signal buffer length, ...
dan@0: %
dan@0: % For any reasonable use, Q equals R, otherwise it makes little sense.
dan@0: % Kernels must include the window function...
dan@0: % 
dan@0: 
dan@0: R = size(krnls,2);
dan@0: Q = size(mf_ders,2);
dan@0: assert(R == size(krlns_ders, 2) );
dan@0: assert(R >= Q);
dan@0: % constructing the matrixes A and B from equation III.4 in [1]
dan@0: % 1st dimension is the discrete time
dan@0: 
dan@0: sig_mat   = repmat(sig,   [1,R,Q]);
dan@0: krnls_mat = repmat(krnls, [1,1,Q]);
dan@0: 
dan@0: mf_ders_mat = repmat(reshape(mf_ders,[N,1,Q]),[1,R,1]);
dan@0: 
dan@0: % inner product for left and right hand side of the eq.
dan@0: A = shiftdim(sum(  conj(krnls_mat)  .* mf_ders_mat .* sig_mat, 1), 1);
dan@0: b = shiftdim(sum(- conj(krlns_ders) .* sig_mat(:,:,1), 1), 1);
dan@0: 
dan@0: end