dan@0: % multi_freq_ddm_int tests..
dan@0: N = 2^10-1;
dan@0: N_fft = 2^12-1;
dan@0: 
dan@0: fs = 44100;
dan@0: 
dan@0: dsc_idxs = ([-(N-1)/2:(N-1)/2]);
dan@0: 
dan@0: frqs_fft = [0:N_fft-1]*fs/N_fft;
dan@0: 
dan@0: the_frq = 416;
dan@0: the_fm1 = 7000;
dan@0: the_ph = 0.2;
dan@0: the_amp = 0.1;
dan@0: pk_offst_rt = 0.01;
dan@0: %frqs = [438,442.5] * 2 *pi; % as many frequencies as desired...
dan@0: fft_bin_hz = fs/N_fft;
dan@0: t = [-(N-1)/2:(N-1)/2]/fs;
dan@0: % monomials matlab style...
dan@0: ply_len = 2;
dan@0: mdl_poly = flipud(diag(ones(ply_len+1,1))); % vector of polynomials
dan@0: mdl_poly_der = zeros(size(mdl_poly));% vector of time derivated polynomials
dan@0: 
dan@0: % the actual time domain values of time derivated polynomials
dan@0: mdl_poly_der_vals = zeros(N,size(mdl_poly,1)); 
dan@0: 
dan@0: % window and time derivated window...
dan@0: w  = 0.5+0.5*cos(2*pi*dsc_idxs/(N-1));
dan@0: w_d = -pi*fs/(N-1)*sin(2*pi*dsc_idxs/(N-1));
dan@0: w = w(:);
dan@0: w_d = w_d(:);
dan@0: 
dan@0: % calculate time domain plys...
dan@0: for k = 1:size(mdl_poly,1)
dan@0:   pp = polyder(mdl_poly(k,:));
dan@0:   mdl_poly_der(k,:) = [zeros(1,ply_len+1-size(pp,2)) pp ];
dan@0:   mdl_poly_der_vals(:,k) = polyval( mdl_poly_der(k,:),t);
dan@0: end
dan@0: 
dan@0: mdl_poly_der_vals = mdl_poly_der_vals(:,2:end);
dan@0: 
dan@0: sig_r = real(exp(the_amp+1j*(the_fm1*t.^2 + t*the_frq*2*pi + the_ph*pi))).';
dan@0: %sig = hilbert(sig_r);
dan@0: sig = sig_r;
dan@0: 
dan@0: [id_s id_e ] = zpzh_idxs(N);
dan@0: sig_zp = sig .* w;
dan@0: sig_zp = [sig_zp(id_s);zeros(N_fft-N,1);sig_zp(id_e)];
dan@0: 
dan@0: sig_fft = fft(sig_zp, N_fft);
dan@0: [pk_vl,pk_idx] = max(abs(sig_fft));
dan@0: pk_frq_hz = (pk_idx - 1) * fs / N_fft;
dan@0: %pk_offst_rt = (fft_bin_hz/the_frq);
dan@0: %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...
dan@0: frqs_test = [pk_idx - 2, pk_idx - 1, pk_idx ] * 2 * pi * fs / N_fft; 
dan@0: %frqs_test = [pk_idx - 0.5, pk_idx - 1.5, pk_idx ] * 2 * pi * fs / N_fft; 
dan@0: 
dan@0: 
dan@0: [fks fks_drs] = fourier_kernel_win(frqs_test, N, fs, w, w_d);
dan@0: [gdi, A, b,gdi_lns,r_lns] = ...
dan@0:   gen_ddm_int(fks, fks_drs, mdl_poly_der_vals, sig, N, 0);
dan@0: 
dan@0: errs = A(:,1) - sig_fft(pk_idx-1:pk_idx+1);
dan@0: assert( sum(errs < (1 + 1j) * 10^(-11)) == 3);
dan@0: 
dan@0: %gdi = gdi_lns;
dan@0: 
dan@0: err_rel_frq_log = log10(abs((imag(gdi(1)) / 2 / pi - the_frq)/(imag(gdi(1)) / 2 / pi) ));
dan@0: err_rel_fm1_log = log10(abs((imag(gdi(2)) - the_fm1) / imag(gdi(2)) ));
dan@0: assert(err_rel_frq_log < -3);
dan@0: assert(err_rel_fm1_log < -2);
dan@0: real(gdi(1)) 
dan@0: real(gdi(2))
dan@0: ply_est = fliplr(gdi.');
dan@0: frq_est = imag(ply_est(end));
dan@0: the_fm1
dan@0: imag(ply_est(end-1))
dan@0: frq_est/2/pi
dan@0: %ply_est(end) = real(ply_est(end)) + 1j*(imag(ply_est(end)) - the_frq);
dan@0: ply_est(end) = real(ply_est(end)); % hm....
dan@0: wmp = win_mod_poly(N, t, [ply_est 0], w);
dan@0: 
dan@0: % this creates a bias but it's fast...
dan@0: 
dan@0: krnl_dmd = exp(1j * frq_est * t');
dan@0: dtft_dmd = sum(w .* sig .* conj(krnl_dmd),1);
dan@0: r0 = dtft_dmd / wmp;
dan@0: err_rel_amp_log = log10(abs((abs(r0)*2 - exp(the_amp)) / exp(the_amp)) )
dan@0: err_rel_phs_log = log10(abs((angle(r0) / pi - the_ph)/the_ph));
dan@0: 
dan@0: assert(err_rel_amp_log  < -2);
dan@0: sprintf('%15.15f', (angle(r0) / pi-the_ph) )
dan@0: 
dan@0: 
dan@0: 
dan@0: 
dan@0: 
dan@0: %% plots
dan@0: figure(3);clf;
dan@0: subplot(211);
dan@0: hold on;
dan@0: plot(frqs_fft,abs(sig_fft),'.k-');
dan@0: plot(frqs_test'/2/pi,abs(A(:,1)),'r*');
dan@0: subplot(212);
dan@0: hold on
dan@0: plot(frqs_fft,angle(sig_fft),'.k-');
dan@0: plot(frqs_test'/2/pi,angle(A(:,1)),'r*');
dan@0: 
dan@0: display('tests completed successfully!');
dan@0: display('accuracy:');
dan@0: display([ '     frq: ' num2str(err_rel_frq_log)]);
dan@0: display([ '     amp: ' num2str(err_rel_amp_log)]);
dan@0: display([ '     fm1: ' num2str(err_rel_fm1_log)]);
dan@0: display([ '     phs: ' num2str(err_rel_phs_log)]);
dan@0: display([ '     rcond: ' num2str(r_lns)]);
dan@0: 
dan@0: