# HG changeset patch # User Dan Stowell # Date 1368126255 -3600 # Node ID a4a7e34050624c0fb53b5e742798795f072ca9ac Import DDM code by Sašo Muševič diff -r 000000000000 -r a4a7e3405062 ddm_der2_lin_sys_fft.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_der2_lin_sys_sol.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_der_lin_sys_fft.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_der_lin_sys_sol.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_fft.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_fft_2.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_fft_3.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_gen.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_lin_sys.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_lin_sys_fft.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_lin_sys_fft_frq_hop.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_lin_sys_fft_single.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_lin_sys_fft_zp_fact.m --- /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 diff -r 000000000000 -r a4a7e3405062 ddm_lin_sys_sol.m --- /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 diff -r 000000000000 -r a4a7e3405062 fourier_kernel.m --- /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 diff -r 000000000000 -r a4a7e3405062 fourier_kernel_win.m --- /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 diff -r 000000000000 -r a4a7e3405062 gen_ddm_2.m --- /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) + + diff -r 000000000000 -r a4a7e3405062 gen_ddm_3.m --- /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) + + diff -r 000000000000 -r a4a7e3405062 gen_ddm_fft.m --- /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 = []; + + diff -r 000000000000 -r a4a7e3405062 gen_ddm_int.m --- /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 diff -r 000000000000 -r a4a7e3405062 lin_solve_dgr_2.m --- /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 diff -r 000000000000 -r a4a7e3405062 lin_solve_dgr_3.m --- /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 diff -r 000000000000 -r a4a7e3405062 lse_stat_par.m --- /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 diff -r 000000000000 -r a4a7e3405062 multi_freq_ddm_int.m --- /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 diff -r 000000000000 -r a4a7e3405062 multi_freq_ddm_int_test_1.m --- /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)]); + + diff -r 000000000000 -r a4a7e3405062 overlap_log_p_1.m --- /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 diff -r 000000000000 -r a4a7e3405062 stat_pars_ddm.m --- /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 diff -r 000000000000 -r a4a7e3405062 win_mod_poly.m --- /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 diff -r 000000000000 -r a4a7e3405062 zpzh_idxs.m --- /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