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
|