Mercurial > hg > ddm

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_der2_lin_sys_fft.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,34 @@
+function [A_sys,b_sys,sig_fft] = ddm_der2_lin_sys_fft(Q, R, win, win_der, mf_ders, sig, N, N_fft, fs)
+
+  assert(size(mf_ders,1) == N );
+  assert(size(mf_ders, 2) == R);
+  assert(size(win,1) == N );
+  assert(size(win,2) == 1 );
+  assert(size(win_der,1) == N );
+  assert(size(win_der,2) == 1 );
+  assert(size(sig,1) == N );
+  assert(size(sig,2) == 1 );
+  assert(N <= N_fft);
+
+  win_mat  = repmat(win, 1, R);
+  frqs_fft = [0:N_fft-1]'*fs*2*pi/N_fft;
+  [str_idx end_idx] = zpzh_idxs(N);
+  sig_mat  = repmat(sig, 1, R);
+  t = disc_time_axis(N)/fs;
+  fft_sig_win_bffr     = win .* sig * 2.* t';
+  fft_sig_win          = fft([fft_sig_win_bffr(str_idx,:);zeros(N_fft-N,1);fft_sig_win_bffr(end_idx,:)],N_fft,1);
+  fft_sig_mat_bffr     = win_mat .* sig_mat .* mf_ders; % S_t^2w, S_t^3w
+  fft_sig_win_der_bffr = win_der .* sig .* mf_ders(:,1); % S_t^2w'
+  fft_sig_win_der      = fft([fft_sig_win_der_bffr(str_idx,:);zeros(N_fft-N,1);fft_sig_win_der_bffr(end_idx,:)],N_fft,1);
+  A =   fft([fft_sig_mat_bffr(str_idx,:);zeros(N_fft-N,R);fft_sig_mat_bffr(end_idx,:)], N_fft, 1); %zero padded sheezl
+  %b = -(-1j*fft_sig_win_der - 1j* fft_sig_win - frqs_fft.* A(:,1));
+  b = -(-fft_sig_win_der - 2 * fft_sig_win + 1j*frqs_fft.* A(:,1));
+% create the actual matrixes
+  A_sys = zeros(Q,R,N_fft-Q+1);
+  b_sys = zeros(Q,1,N_fft-Q+1);
+  for k=1:Q
+    A_sys(k,:,:) = -1j*shiftdim(A(k:N_fft-Q+k,:).',-1);
+    b_sys(k,:,:) =     shiftdim(b(k:N_fft-Q+k).'  ,-1);
+  end
+  sig_fft = A(:,1);
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_der2_lin_sys_sol.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,38 @@
+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
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_der_lin_sys_fft.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,32 @@
+function [A_sys,b_sys,sig_fft] = ddm_der_lin_sys_fft(Q, R, win, win_der, mf_ders, sig, N, N_fft, fs)
+
+  assert(size(mf_ders,1) == N );
+  assert(size(mf_ders, 2) == R);
+  assert(size(win,1) == N );
+  assert(size(win,2) == 1 );
+  assert(size(win_der,1) == N );
+  assert(size(win_der,2) == 1 );
+  assert(size(sig,1) == N );
+  assert(size(sig,2) == 1 );
+  assert(N <= N_fft);
+
+  win_mat  = repmat(win, 1, R);
+  frqs_fft = [0:N_fft-1]'*fs*2*pi/N_fft;
+  [str_idx end_idx] = zpzh_idxs(N);
+  sig_mat  = repmat(sig, 1, R);
+  fft_sig_win_bffr     = win .* sig;
+  fft_sig_win          = fft([fft_sig_win_bffr(str_idx,:);zeros(N_fft-N,1);fft_sig_win_bffr(end_idx,:)],N_fft,1);
+  fft_sig_mat_bffr     = win_mat .* sig_mat .* mf_ders; % S_tw, S_t^2w
+  fft_sig_win_der_bffr = win_der .* sig .* mf_ders(:,1); % S_tw'
+  fft_sig_win_der = fft([fft_sig_win_der_bffr(str_idx,:);zeros(N_fft-N,1);fft_sig_win_der_bffr(end_idx,:)],N_fft,1);
+  A =   fft([fft_sig_mat_bffr(str_idx,:);zeros(N_fft-N,R);fft_sig_mat_bffr(end_idx,:)], N_fft, 1); %zero padded sheezl
+  b = -(-1j*fft_sig_win_der - 1j* fft_sig_win - frqs_fft.* A(:,1));
+% create the actual matrixes
+  A_sys = zeros(Q,R,N_fft-Q+1);
+  b_sys = zeros(Q,1,N_fft-Q+1);
+  for k=1:Q
+    A_sys(k,:,:) = -1j*shiftdim(A(k:N_fft-Q+k,:).',-1);
+    b_sys(k,:,:) =     shiftdim(b(k:N_fft-Q+k).',-1);
+  end
+  sig_fft = A(:,1);
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_der_lin_sys_sol.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,23 @@
+function [ddlss dlss A_plus AA_inv A_der_plus] = ddm_der_lin_sys_sol(A_sys,b_sys,A_der_sys,b_der_sys,N_fft,R,Q)
+dlss = zeros(Q+1, N_fft-R+1);
+ddlss = 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];
+
+end
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_fft.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,45 @@
+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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_fft_2.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,12 @@
+function [df2,A_sys,b_sys] = ddm_fft_2(win, win_der, mf_ders, sig, N, N_fft, fs)
+  R = 2;
+  Q = R; % nr of unknows is the same as the number of model funcs, fact!
+  [A_sys b_sys] =  ddm_lin_sys_fft(Q,R,win, win_der, mf_ders, sig, N, N_fft, fs);
+%    A_sys1 = cat(2,shiftdim(A(1:N_fft-1,1),-2),shiftdim(A(1:N_fft-1,2),-2));
+%    A_sys2 = cat(2,shiftdim(A(2:N_fft  ,1),-2),shiftdim(A(2:N_fft  ,2),-2));
+%    A_sys  = cat(1,A_sys1,A_sys2);
+%    b_sys  = cat(1,shiftdim(b(1:N_fft-1).',-1),shiftdim(b(2:N_fft).',  -1));
+%    A_sys = A_sys_;
+%    b_sys = b_sys_;
+  df2    = lin_solve_dgr_2(A_sys,b_sys,N_fft-1); %hardcoded degree 2 solver (fast)
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_fft_3.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,12 @@
+function [df3,A_sys,b_sys] = ddm_fft_3(win, win_der, mf_ders, sig, N, N_fft, fs)
+  R = 3;
+  Q = R; % nr of unknows is the same as the number of model funcs
+  [A_sys b_sys] =  ddm_lin_sys_fft(Q, R, win, win_der, mf_ders, sig, N, N_fft, fs);
+%   A_sys1 = cat(2,shiftdim(A(1:N_fft-2,1),-2),shiftdim(A(1:N_fft-2,2),-2),shiftdim(A(1:N_fft-2,3),-2));
+%   A_sys2 = cat(2,shiftdim(A(2:N_fft-1,1),-2),shiftdim(A(2:N_fft-1,2),-2),shiftdim(A(2:N_fft-1,3),-2));
+%   A_sys3 = cat(2,shiftdim(A(3:N_fft  ,1),-2),shiftdim(A(3:N_fft  ,2),-2),shiftdim(A(3:N_fft  ,3),-2));
+%   A_sys  = cat(1,A_sys1,A_sys2,A_sys3);
+%   b_sys  = cat(1,shiftdim(b(1:N_fft-2).',-1),shiftdim(b(2:N_fft-1).',  -1),shiftdim(b(3:N_fft  ).',  -1));
+  df3    =  lin_solve_dgr_3(A_sys,b_sys,N_fft-2); %hardcoded degree 3 solver (fast)
+
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_gen.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,40 @@
+% a high-level function for computing a DTFT based DDM
+% (constructs the linear system of equations using distribution derivative rule and solves it via
+% pseudo-inverse (pinv). Use this function if non-Fourier kernels are needed or
+% DTFT at specific frequency rather than FFT is required for whatever reason...
+% no zero-padding is used internally
+function [gdi,r, A,b, gdi_lns, r_lns] = ddm_gen(krnls, krlns_ders, mf_ders, sig, N, tol)
+%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...
+
+[A,b] = ddm_lin_sys(krnls, krlns_ders, mf_ders, sig, N);
+
+% solving via pinv with provided tolerance
+gdi = pinv(A,tol) * b;
+% r = rcond(A);
+% [gdi_lns r_lns]  = linsolve(A,b);
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_lin_sys.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,47 @@
+% 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
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_lin_sys_fft.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,34 @@
+% create
+function [A_sys,b_sys,sig_fft] = ddm_lin_sys_fft(Q, R, win, win_der, mf_ders, sig, N, N_fft, fs)
+  %% START: comment this out for better performance
+  assert(size(mf_ders,1) == N );
+  assert(size(mf_ders, 2) == R);
+  assert(size(win,1) == N );
+  assert(size(win,2) == 1 );
+  assert(size(win_der,1) == N );
+  assert(size(win_der,2) == 1 );
+  assert(size(sig,1) == N );
+  assert(size(sig,2) == 1 );
+  assert(N <= N_fft);
+  %% START: comment this out for better performance
+
+  win_mat  = repmat(win, 1, R);
+  frqs_fft = [0:N_fft-1]'*fs*2*pi/N_fft;
+  [str_idx end_idx] = zpzh_idxs(N);
+  sig_mat  = repmat(sig, 1, R);
+
+
+  fft_sig_mat_bffr     = win_mat .* sig_mat .* mf_ders;
+  fft_sig_win_der_bffr = win_der .* sig; % usefull hack to easily compute derivative of the matrix
+  fft_sig_win_der = fft([fft_sig_win_der_bffr(str_idx,:);zeros(N_fft-N,1);fft_sig_win_der_bffr(end_idx,:)],N_fft,1);
+  A =   fft([fft_sig_mat_bffr(str_idx,:);zeros(N_fft-N,R);fft_sig_mat_bffr(end_idx,:)], N_fft, 1); %zero padded sheezl
+  b = -(fft_sig_win_der - 1j*frqs_fft.* A(:,1));
+% create the actual matrixes
+  A_sys = zeros(Q,R,N_fft-Q+1);
+  b_sys = zeros(Q,1,N_fft-Q+1);
+  for k=1:Q
+    A_sys(k,:,:) = shiftdim(A(k:N_fft-Q+k,:).',-1);
+    b_sys(k,:,:) = shiftdim(b(k:N_fft-Q+k).',-1);
+  end
+  sig_fft = A(:,1);
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_lin_sys_fft_frq_hop.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,35 @@
+function [A_sys,b_sys] = ddm_lin_sys_fft_frq_hop(R, Q, frq_hop, win, win_der, mf_ders, sig, N, N_fft, fs)
+
+  assert(size(mf_ders,1) == N );
+  assert(size(mf_ders, 2) == Q);
+  assert(size(win,1) == N );
+  assert(size(win,2) == 1 );
+  assert(size(win_der,1) == N );
+  assert(size(win_der,2) == 1 );
+  assert(size(sig,1) == N );
+  assert(size(sig,2) == 1 );
+  assert(N <= N_fft);
+
+  win_mat  = repmat(win, 1, Q);
+  frqs_fft = [0:N_fft-1]'*fs*2*pi/N_fft;
+  [str_idx end_idx] = zpzh_idxs(N);
+  sig_mat  = repmat(sig, 1, Q);
+  fft_sig_mat_bffr     = win_mat .* sig_mat .* mf_ders;
+  fft_sig_win_der_bffr = win_der .* sig;
+  fft_sig_win_der = fft([fft_sig_win_der_bffr(str_idx,:);zeros(N_fft-N,1);fft_sig_win_der_bffr(end_idx,:)],N_fft,1);
+  A =   fft([fft_sig_mat_bffr(str_idx,:);zeros(N_fft-N,Q);fft_sig_mat_bffr(end_idx,:)], N_fft, 1); %zero padded sheeza
+  b = -(fft_sig_win_der - 1j*frqs_fft.* A(:,1));
+% create the actual matrixes
+
+  frq_bins = [1:frq_hop:N_fft-R+1];
+  frq_lims = [frq_bins;frq_bins+R];
+  nr_frqs = length(frq_bins);
+  A_sys = zeros(R,Q,nr_frqs);
+  b_sys = zeros(R,1,nr_frqs);
+
+  for k=1:frq_hop:R
+    A_sys(k,:,:) = shiftdim(A(k:N_fft-R+k,:).',-1);
+    b_sys(k,:,:) = shiftdim(b(k:N_fft-R+k).',-1);
+  end
+
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_lin_sys_fft_single.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,26 @@
+function [A_sys,b_sys] = ddm_lin_sys_fft_single(ddm_bins, R, win, win_der, mf_ders, sig, N, N_fft, fs)
+
+  assert(size(mf_ders,1) == N );
+  assert(size(mf_ders, 2) == R);
+  assert(size(win,1) == N );
+  assert(size(win,2) == 1 );
+  assert(size(win_der,1) == N );
+  assert(size(win_der,2) == 1 );
+  assert(size(sig,1) == N );
+  assert(size(sig,2) == 1 );
+  assert(N <= N_fft);
+
+  Q = length(ddm_bins);
+
+  win_mat  = repmat(win, 1, R);
+  frqs_fft = [0:N_fft-1]'*fs*2*pi/N_fft;
+  [str_idx end_idx] = zpzh_idxs(N);
+  sig_mat  = repmat(sig, 1, R);
+  fft_sig_mat_bffr     = win_mat .* sig_mat .* mf_ders;
+  fft_sig_win_der_bffr = win_der .* sig;
+  fft_sig_win_der = fft([fft_sig_win_der_bffr(str_idx,:);zeros(N_fft-N,1);fft_sig_win_der_bffr(end_idx,:)],N_fft,1);
+  A =   fft([fft_sig_mat_bffr(str_idx,:);zeros(N_fft-N,R);fft_sig_mat_bffr(end_idx,:)], N_fft, 1); %zero padded sheeza
+  b = -(fft_sig_win_der - 1j*frqs_fft.* A(:,1));
+  A_sys = A(ddm_bins,:);
+  b_sys = b(ddm_bins);
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_lin_sys_fft_zp_fact.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,25 @@
+function [A b A_sys,b_sys] = ddm_lin_sys_fft_zp_fact(Q, R, win, win_der, mf_ders, sig, N, zp_fct, fs,all_zp_bins)
+
+  assert(R == size(mf_ders, 2) );
+  N_fft = N*zp_fct;
+  win_mat  = repmat(win, 1, R);
+  frqs_fft = [0:N_fft-1]'*fs*2*pi/N_fft;
+  [str_idx end_idx] = zpzh_idxs(N);
+  sig_mat  = repmat(sig, 1, R);
+  fft_sig_mat_bffr     = win_mat .* sig_mat .* mf_ders;
+  fft_sig_win_der_bffr = win_der .* sig;
+  fft_sig_win_der = fft([fft_sig_win_der_bffr(str_idx,:);zeros(N_fft-N,1);fft_sig_win_der_bffr(end_idx,:)],N_fft,1);
+  A =   fft([fft_sig_mat_bffr(str_idx,:);zeros(N_fft-N,R);fft_sig_mat_bffr(end_idx,:)], N_fft, 1); %zero padded sheeza
+  b = -(fft_sig_win_der - 1j*frqs_fft.* A(:,1));
+% create the actual matrixes
+  A_sys = zeros(Q,Q,N_fft-Q+1);
+  b_sys = zeros(Q,1,N_fft-Q+1);
+  stp = zp_fct;
+  if ~exists('all_cp_bins','var') || all_cp_bins == 1
+    stp = 1;
+  end
+  for k=1:Q
+    A_sys(k,:,:) = shiftdim(A(k:stp:N_fft-Q+k,:).',-1);
+    b_sys(k,:,:) = shiftdim(b(k:stp:N_fft-Q+k).',-1);
+  end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_lin_sys_sol.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,9 @@
+function dlss = ddm_lin_sys_sol(A_sys,b_sys,N_fft,R)
+dlss = zeros(ddm_dgr+1, N_fft-R+1);
+
+for l=1:N_fft-R+1
+  A = A_sys(:,:,l);
+  b = b_sys(:,:,l);
+  dlss(:,l)     = [fliplr((inv(A'*A)*A' * b).') 0];
+end
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fourier_kernel.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,9 @@
+function [fk fk_ders] = fourier_kernel(frqs, N,fs)
+R = size(frqs,2);
+t = [-(N-1)/2:(N-1)/2]'/fs;
+
+t_mat = repmat(t,[1,R]);
+frq_mat = repmat(frqs,[N,1]);
+fk = exp(1j .* t_mat .* frq_mat);
+fk_ders = 1j * fk .* frq_mat;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fourier_kernel_win.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,11 @@
+function [fk fk_ders] = fourier_kernel_win(frqs, N,fs, win, win_der)
+R = size(frqs,2);
+t = [- (N - 1) / 2:(N - 1) / 2]' / fs;
+win_mat     = repmat(win,     [1,R]);
+win_der_mat = repmat(win_der, [1,R]);
+t_mat = repmat(t, [1, R]);
+frq_mat = repmat(frqs, [N, 1]);
+krnls = exp(1j .* t_mat .* frq_mat);
+fk = win_mat .* krnls;
+fk_ders = (1j *  win_mat  .* frq_mat + win_der_mat) .* krnls;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gen_ddm_2.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,42 @@
+% distribution derivative method:
+% 2nd degree system solver is used, make sure the data provided corresponds to that
+function dg2 = gen_ddm_2(krnls, krlns_ders, mf_ders, sig, N)
+% 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...
+
+
+R = size(krnls,2);
+assert(R == 2);
+assert(R == size(krlns_ders,2) );
+assert(R == size(mf_ders,2) );
+
+[A, b] = ddm_lin_sys(krnls, krlns_ders, mf_ders, sig, N); % generate the linear system of eqs (slow)
+dg2 = lin_solve_dgr_2(A,b,1); %hardcoded degree 2 solver (fast)
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gen_ddm_3.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,43 @@
+% distribution derivative method:
+% 3rd degree system solver is used, make sure the data provided corresponds to that
+function dg3 = gen_ddm_3(krnls, krlns_ders, mf_ders, sig, N)
+% 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...
+
+% 3rd degree system solver is used, make sure the data provided corresponds
+% to that
+R = size(krnls,2);
+assert(R == 3);
+assert(R == size(krlns_ders,2) );
+assert(R == size(mf_ders,2) );
+
+[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)
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gen_ddm_fft.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,39 @@
+function dg3 = gen_ddm_fft(krnls, krlns_ders, mf_ders, sig, N)
+% 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...
+
+R = size(krnls,2);
+assert(R == size(krlns_ders,2) );
+assert(R == size(mf_ders,2) );
+
+[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)
+dg2 = [];
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gen_ddm_int.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,39 @@
+% a high-level function for computing a DTFT based DDM
+% (constructs the linear system of equations using distribution derivative rule and solves it via
+% pseudo-inverse (pinv). Use this function if non-Fourier kernels are needed or
+% using FFT would be unsatisfactory due to
+function [gdi,A,b, gdi_lns, r_lns] = gen_ddm_int(krnls, krlns_ders, mf_ders, sig, N, tol)
+%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...
+
+[A,b] = ddm_lin_sys(krnls, krlns_ders, mf_ders, sig, N);
+
+% solving via pinv with provided tolerance
+gdi = pinv(A,tol) * b; % more generic than the following code:
+%[gdi_lns r_lns]  = linsolve(A,b);
+%gdi = gdi_lns;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lin_solve_dgr_2.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,7 @@
+% solve multiple 2nd degree linear systems at once - no pivoting!
+function dg2 = lin_solve_dgr_2(A,b,R)
+nrm =  (-A(2,1,:) .* A(1,2,:) + A(1,1,:) .* A(2,2,:));
+  dg2 = zeros(2,R);
+  dg2(1,:) = -(b(2,1,:) .* A(1,2,:) - b(1,1,:) .* A(2,2,:)) ./ nrm;
+  dg2(2,:) = -(b(1,1,:) .* A(2,1,:) - b(2,1,:) .* A(1,1,:)) ./ nrm;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lin_solve_dgr_3.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,26 @@
+% solve multiple 3rd degree linear systems at once - no pivoting!
+function dg3 = lin_solve_dgr_3(A,b,R)
+  nrm =  ...
+    + A(1,3,:) .* A(2,2,:) .* A(3,1,:) ...
+    - A(1,3,:) .* A(2,1,:) .* A(3,2,:) ...
+    + A(1,2,:) .* A(2,1,:) .* A(3,3,:) ...
+    - A(1,2,:) .* A(2,3,:) .* A(3,1,:) ...
+    + A(1,1,:) .* A(2,3,:) .* A(3,2,:) ...
+    - A(1,1,:) .* A(2,2,:) .* A(3,3,:);
+  dg3 = zeros(3,R);
+  dg3(1,:) = -(...
+    +b(3,1,:) .* ( A(1,2,:) .* A(2,3,:) - A(1,3,:) .* A(2,2,:) ) ...
+    -b(2,1,:) .* ( A(1,2,:) .* A(3,3,:) - A(1,3,:) .* A(3,2,:) ) ...
+    +b(1,1,:) .* ( A(2,2,:) .* A(3,3,:) - A(2,3,:) .* A(3,2,:) ) ...
+    ) ./ nrm;
+  dg3(2,:) = (...
+    +b(3,1,:) .* ( A(1,1,:) .* A(2,3,:) - A(1,3,:) .* A(2,1,:) ) ...
+    -b(2,1,:) .* ( A(1,1,:) .* A(3,3,:) - A(1,3,:) .* A(3,1,:) ) ...
+    +b(1,1,:) .* ( A(2,1,:) .* A(3,3,:) - A(2,3,:) .* A(3,1,:) ) ...
+    ) ./ nrm;
+  dg3(3,:) = -(...
+    +b(3,1,:) .* ( A(1,1,:) .* A(2,2,:) - A(1,2,:) .* A(2,1,:) ) ...
+    -b(2,1,:) .* ( A(1,1,:) .* A(3,2,:) - A(1,2,:) .* A(3,1,:) ) ...
+    +b(1,1,:) .* ( A(2,1,:) .* A(3,2,:) - A(2,2,:) .* A(3,1,:) ) ...
+    ) ./ nrm;
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lse_stat_par.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,8 @@
+function lsp = lse_stat_par(sig, krnls, w)
+  nr_krnls = size(krnls,1);
+  sig_mat = repmat(sig,nr_krnls,1);
+  w_mat   = repmat(w  ,nr_krnls,1);
+  kk = sum(  krnls.*conj(krnls).*w_mat,2);
+  ks = sum(sig_mat.*conj(krnls).*w_mat,2);
+  lsp = log(ks/kk*2);
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multi_freq_ddm_int.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,26 @@
+function mfdi = multi_freq_ddm_int(frqs, mf_ders, sig, N, tol)
+%multi-frequency distribution derivative based estimator for
+% non-stationary sinusoidal analysis. Frequencies may not be at the bins of
+% FFT...
+%
+% Uses the Fourier kernel at specified frequencies. No assumptions about
+%
+% [1] Michael Betser: Sinusoidal Polynomial Estimation Using The Distribution
+% Derivative, in IEEE Transactions on Signal Processing, Vol.57, Nr. 12,
+% December 2009
+%
+%
+%
+
+
+R = size(frqs,2);
+t = [-(N-1)/2:(N-1)/2]'/fs;
+
+t_mat = repmat(t,[1,R]);
+frq_mat = repmat(frqs,[N,1]);
+krnls = exp(1j .* t_mat .* frq_mat);
+krlns_ders = 1j * krnls .* frq_mat;
+
+mfdi = gen_ddm_int(krnls, krlns_ders, mf_ders, sig, N, tol);
+
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/multi_freq_ddm_int_test_1.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,117 @@
+% multi_freq_ddm_int tests..
+N = 2^10-1;
+N_fft = 2^12-1;
+
+fs = 44100;
+
+dsc_idxs = ([-(N-1)/2:(N-1)/2]);
+
+frqs_fft = [0:N_fft-1]*fs/N_fft;
+
+the_frq = 416;
+the_fm1 = 7000;
+the_ph = 0.2;
+the_amp = 0.1;
+pk_offst_rt = 0.01;
+%frqs = [438,442.5] * 2 *pi; % as many frequencies as desired...
+fft_bin_hz = fs/N_fft;
+t = [-(N-1)/2:(N-1)/2]/fs;
+% monomials matlab style...
+ply_len = 2;
+mdl_poly = flipud(diag(ones(ply_len+1,1))); % vector of polynomials
+mdl_poly_der = zeros(size(mdl_poly));% vector of time derivated polynomials
+
+% the actual time domain values of time derivated polynomials
+mdl_poly_der_vals = zeros(N,size(mdl_poly,1));
+
+% window and time derivated window...
+w  = 0.5+0.5*cos(2*pi*dsc_idxs/(N-1));
+w_d = -pi*fs/(N-1)*sin(2*pi*dsc_idxs/(N-1));
+w = w(:);
+w_d = w_d(:);
+
+% calculate time domain plys...
+for k = 1:size(mdl_poly,1)
+  pp = polyder(mdl_poly(k,:));
+  mdl_poly_der(k,:) = [zeros(1,ply_len+1-size(pp,2)) pp ];
+  mdl_poly_der_vals(:,k) = polyval( mdl_poly_der(k,:),t);
+end
+
+mdl_poly_der_vals = mdl_poly_der_vals(:,2:end);
+
+sig_r = real(exp(the_amp+1j*(the_fm1*t.^2 + t*the_frq*2*pi + the_ph*pi))).';
+%sig = hilbert(sig_r);
+sig = sig_r;
+
+[id_s id_e ] = zpzh_idxs(N);
+sig_zp = sig .* w;
+sig_zp = [sig_zp(id_s);zeros(N_fft-N,1);sig_zp(id_e)];
+
+sig_fft = fft(sig_zp, N_fft);
+[pk_vl,pk_idx] = max(abs(sig_fft));
+pk_frq_hz = (pk_idx - 1) * fs / N_fft;
+%pk_offst_rt = (fft_bin_hz/the_frq);
+%frqs_test = [pk_frq_hz*(1 - pk_offst_rt),pk_frq_hz, pk_frq_hz*(1 + pk_offst_rt)] * 2 *pi; % as many frequencies as desired...
+frqs_test = [pk_idx - 2, pk_idx - 1, pk_idx ] * 2 * pi * fs / N_fft;
+%frqs_test = [pk_idx - 0.5, pk_idx - 1.5, pk_idx ] * 2 * pi * fs / N_fft;
+
+
+[fks fks_drs] = fourier_kernel_win(frqs_test, N, fs, w, w_d);
+[gdi, A, b,gdi_lns,r_lns] = ...
+  gen_ddm_int(fks, fks_drs, mdl_poly_der_vals, sig, N, 0);
+
+errs = A(:,1) - sig_fft(pk_idx-1:pk_idx+1);
+assert( sum(errs < (1 + 1j) * 10^(-11)) == 3);
+
+%gdi = gdi_lns;
+
+err_rel_frq_log = log10(abs((imag(gdi(1)) / 2 / pi - the_frq)/(imag(gdi(1)) / 2 / pi) ));
+err_rel_fm1_log = log10(abs((imag(gdi(2)) - the_fm1) / imag(gdi(2)) ));
+assert(err_rel_frq_log < -3);
+assert(err_rel_fm1_log < -2);
+real(gdi(1))
+real(gdi(2))
+ply_est = fliplr(gdi.');
+frq_est = imag(ply_est(end));
+the_fm1
+imag(ply_est(end-1))
+frq_est/2/pi
+%ply_est(end) = real(ply_est(end)) + 1j*(imag(ply_est(end)) - the_frq);
+ply_est(end) = real(ply_est(end)); % hm....
+wmp = win_mod_poly(N, t, [ply_est 0], w);
+
+% this creates a bias but it's fast...
+
+krnl_dmd = exp(1j * frq_est * t');
+dtft_dmd = sum(w .* sig .* conj(krnl_dmd),1);
+r0 = dtft_dmd / wmp;
+err_rel_amp_log = log10(abs((abs(r0)*2 - exp(the_amp)) / exp(the_amp)) )
+err_rel_phs_log = log10(abs((angle(r0) / pi - the_ph)/the_ph));
+
+assert(err_rel_amp_log  < -2);
+sprintf('%15.15f', (angle(r0) / pi-the_ph) )
+
+
+
+
+
+%% plots
+figure(3);clf;
+subplot(211);
+hold on;
+plot(frqs_fft,abs(sig_fft),'.k-');
+plot(frqs_test'/2/pi,abs(A(:,1)),'r*');
+subplot(212);
+hold on
+plot(frqs_fft,angle(sig_fft),'.k-');
+plot(frqs_test'/2/pi,angle(A(:,1)),'r*');
+
+display('tests completed successfully!');
+display('accuracy:');
+display([ '     frq: ' num2str(err_rel_frq_log)]);
+display([ '     amp: ' num2str(err_rel_amp_log)]);
+display([ '     fm1: ' num2str(err_rel_fm1_log)]);
+display([ '     phs: ' num2str(err_rel_phs_log)]);
+display([ '     rcond: ' num2str(r_lns)]);
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/overlap_log_p_1.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,23 @@
+fs = 441;;
+N = 2^10-1;
+
+ply = [1000, -65, -1.2];
+
+t = [-(N-1)/2:(N-1)/2];
+
+p_val = polyval(ply,t);
+f = log(p_val);
+f_fit = polyfit(t,f,20);
+real(f_fit(end - 3:end))
+
+ply(end-1)/ply(end)
+log(ply(end))
+
+figure(1);clf;
+subplot(311);
+plot(p_val);
+subplot(312);hold on
+plot(real(f));
+plot(real(polyval(f_fit,t)),'r--');
+subplot(313);
+plot(imag(f)/pi*-2+1);
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stat_pars_ddm.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,11 @@
+function [ests r0] = stat_pars_ddm(sig,gdi,t,N,w)
+frq_est = imag(gdi(1));
+ply_est = fliplr(gdi.');
+ply_est(end) = real(ply_est(end)); % hm....
+wmp = win_mod_poly(N, t, [ply_est 0], w);
+krnl_dmd = exp(1j * frq_est * t');
+dtft_dmd = sum(w .* sig .* conj(krnl_dmd),1);
+r0 = dtft_dmd / wmp * 2;
+%ply_est(end) = ply_est(end) + 1j*frq_est;
+ests = [fliplr(gdi.') abs(r0) + 1j*angle(r0)];
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/win_mod_poly.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,9 @@
+function wmp = win_mod_poly(N, t, poly, w)
+R = size(poly,1);
+ply_val = zeros(N,R);
+w_mat = repmat (w, 1, R);
+for k = 1:R
+  ply_val(:,k) = polyval(poly(k,:),t).';
+end
+wmp = sum (w_mat.*exp(ply_val),1);
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/zpzh_idxs.m	Thu May 09 20:04:15 2013 +0100
@@ -0,0 +1,11 @@
+function [is, ie] = zpzh_idxs(N)
+md2 = mod(N,2);
+if md2 == 1
+  is =  (N + md2)/2:N;
+  ie = 1:(N + md2 - 1)/2;
+else
+  is =  N/2+1:N;
+  ie = [1:N /2];
+end
+
+end
\ No newline at end of file