annotate ddm_example_1.m @ 4:72c011ed1977 tip

more elaborate example with non-stat. estimate explanation
author smusevic
date Tue, 30 Jul 2013 09:56:27 +0100
parents e2b116f3b69b
children
rev   line source
smusevic@2 1 input_file = 'path_to_file';
smusevic@2 2 %[sig,fs] = wavread(input_file);
smusevic@2 3
smusevic@2 4 %%%%% DEMO %%%%%
smusevic@2 5 fs = 44100;
smusevic@2 6 sig = rand(100000,1);
smusevic@2 7 %%%%% DEMO %%%%%
smusevic@2 8
smusevic@2 9 % frame length
smusevic@2 10 N = 2^9;
smusevic@2 11 % DDM internal FFT size, for zero-padding, must be => N
smusevic@2 12 N_fft = N;
smusevic@2 13 % quantization, can be anything, but N_q = N makes the most sense
smusevic@2 14 N_q = N;
smusevic@2 15 % degree 2 DDM (linear FM, quadratic AM)
smusevic@2 16 ddm_dgr = 2;
smusevic@2 17 % DDM bandwidth (in bins)
smusevic@2 18 nr_bins = 5;
smusevic@2 19 % AM/FM model functions - 2nd degree polynomial
smusevic@2 20 dsc_idxs = [-(N-1)/2:(N-1)/2];
smusevic@2 21 t = dsc_idxs / fs;
smusevic@2 22 poly_model = poly_mat(t,ddm_dgr,N);
smusevic@2 23
smusevic@2 24 % hann and hann' window
smusevic@2 25 w = 0.5+0.5*cos(2*pi*dsc_idxs/(N-1));
smusevic@2 26 wd = -pi/fs*(N-1)*sin(2*pi*dsc_idxs/(N-1));
smusevic@2 27 % example
smusevic@2 28 offset = 306; % some sample offset in the audio file
smusevic@2 29 bffr = sig(offset:offset+N-1); % get the buffer of lentgh N
smusevic@2 30
smusevic@2 31
smusevic@2 32 %%% execute DDM on the whole frequency range (FFT) %%%
smusevic@2 33 % create linear systems
smusevic@2 34 [A, b, sig_fft] = ddm_lin_sys_fft(nr_bins, ddm_dgr, w', wd', poly_model', bffr, N, N_fft, fs);
smusevic@2 35 % fast pseudo-inverse
smusevic@2 36 ddm_estimates = flipud(pinv_2_fast(A,b,nr_bins,N-nr_bins+1).');
smusevic@4 37
smusevic@4 38 % phase function estimates
smusevic@4 39 ddm_frequency_estimates_for_each_fft_bin = imag(ddm_estimates(2,:));
smusevic@4 40 ddm_linear_fm_estimates_for_each_fft_bin = imag(ddm_estimates(1,:));
smusevic@4 41 % amplitude funciton estimates
smusevic@4 42 ddm_linear_am_estimates_for_each_fft_bin = real(ddm_estimates(2,:));
smusevic@4 43 ddm_2nd_dgr_am_estimates_for_each_fft_bin = real(ddm_estimates(1,:));
smusevic@4 44
smusevic@2 45 % re-assign the spectral magnitude
smusevic@4 46 ddm_spectrogram = reallocate(abs(sig_fft), ddm_estimates, N_q, fs);
smusevic@4 47
smusevic@4 48 % find the peak
smusevic@4 49 [peak_val peak_bin] = max(ddm_spectrogram);
smusevic@4 50
smusevic@4 51 % frequency trajectory of the peak bin
smusevic@4 52 frequency_trajectory = ddm_linear_fm_estimates_for_each_fft_bin(peak_bin) * t + ddm_frequency_estimates_for_each_fft_bin(peak_bin);
smusevic@4 53
smusevic@4 54