smusevic@2: input_file = 'path_to_file';
smusevic@2: %[sig,fs] = wavread(input_file);
smusevic@2: 
smusevic@2: %%%%% DEMO %%%%%
smusevic@2: fs = 44100;
smusevic@2: sig = rand(100000,1);
smusevic@2: %%%%% DEMO %%%%%
smusevic@2: 
smusevic@2: % frame length
smusevic@2: N = 2^9;
smusevic@2: % DDM internal FFT size, for zero-padding, must be => N
smusevic@2: N_fft = N;
smusevic@2: % quantization, can be anything, but N_q = N makes the most sense
smusevic@2: N_q = N;
smusevic@2: % degree 2 DDM (linear FM, quadratic AM)
smusevic@2: ddm_dgr = 2;
smusevic@2: % DDM bandwidth (in bins)
smusevic@2: nr_bins = 5;
smusevic@2: % AM/FM model functions - 2nd degree polynomial
smusevic@2: dsc_idxs = [-(N-1)/2:(N-1)/2];
smusevic@2: t = dsc_idxs / fs;
smusevic@2: poly_model = poly_mat(t,ddm_dgr,N);
smusevic@2: 
smusevic@2: % hann and hann' window
smusevic@2: w = 0.5+0.5*cos(2*pi*dsc_idxs/(N-1));
smusevic@2: wd = -pi/fs*(N-1)*sin(2*pi*dsc_idxs/(N-1));
smusevic@2: % example
smusevic@2: offset = 306; % some sample offset in the audio file
smusevic@2: bffr = sig(offset:offset+N-1); % get the buffer of lentgh N 
smusevic@2: 
smusevic@2: 
smusevic@2: %%% execute DDM on the whole frequency range (FFT) %%%
smusevic@2: % create linear systems
smusevic@2: [A, b, sig_fft] = ddm_lin_sys_fft(nr_bins, ddm_dgr, w', wd', poly_model', bffr, N, N_fft, fs);
smusevic@2: % fast pseudo-inverse
smusevic@2: ddm_estimates = flipud(pinv_2_fast(A,b,nr_bins,N-nr_bins+1).');
smusevic@4: 
smusevic@4: % phase function estimates
smusevic@4: ddm_frequency_estimates_for_each_fft_bin = imag(ddm_estimates(2,:));
smusevic@4: ddm_linear_fm_estimates_for_each_fft_bin = imag(ddm_estimates(1,:));
smusevic@4: % amplitude funciton estimates
smusevic@4: ddm_linear_am_estimates_for_each_fft_bin = real(ddm_estimates(2,:));
smusevic@4: ddm_2nd_dgr_am_estimates_for_each_fft_bin = real(ddm_estimates(1,:));
smusevic@4: 
smusevic@2: % re-assign the spectral magnitude 
smusevic@4: ddm_spectrogram = reallocate(abs(sig_fft), ddm_estimates, N_q, fs);
smusevic@4: 
smusevic@4: % find the peak
smusevic@4: [peak_val peak_bin] = max(ddm_spectrogram);
smusevic@4: 
smusevic@4: % frequency trajectory of the peak bin
smusevic@4: frequency_trajectory = ddm_linear_fm_estimates_for_each_fft_bin(peak_bin) * t + ddm_frequency_estimates_for_each_fft_bin(peak_bin);
smusevic@4: 
smusevic@4: