view ddm_fft.m @ 4:72c011ed1977 tip

more elaborate example with non-stat. estimate explanation

function [dg rs] = ddm_fft(R,w, w_d, mf_ders,sig, N, N_fft,fs)
% multi-frequency distribution derivative based non-stationary sinusoidal estimator 
%  can estimate only 1 sinusoid at the once
%
%
% [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 x K, where R is the number of 
%      non-static parameters to estimate and at the same time, the number 
%      of kernels for each sinusoid, K
%
% krlns_ders: matrix of all the kernel time derivatives... N x R x K , 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, ...
%
% K: number of sinusoids to estimate - NOT OVERLAPPING!!!
%
% For any reasonable use, Q equals R, otherwise it makes little sense.
% Kernels must include the window function...


Q = size(mf_ders,2);

[A_,b_,As, bs] = ddm_lin_sys_fft(R, size(mf_ders, 2), w, w_d, mf_ders, sig, N, N_fft, fs);
%[A, b] = ddm_lin_sys(krnls, krlns_ders, mf_ders, sig, N); % generate the linear system of eqs (slow)
%dg2 = lin_solve_dgr_3(A,b,1); %hardcoded degree 2 solver (fast)
dg = zeros(Q,N_fft-R+1);
rs = zeros(1,N_fft-R+1);
for k = 1:N_fft-R+1
  A = As(:,:,k);
  b = bs(:,k);
  Asq = A.'*A;
  rs(k) = rcond(Asq);
  dg(:,k) = inv(Asq)*A.'*b;
end