changeset 2:e2b116f3b69b

fixed reallocate and pinv_2_fast to work on Octave
added short example very quickly, might not be working properly

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ddm_example_1.m	Thu Jul 25 13:07:07 2013 +0100
@@ -0,0 +1,38 @@
+input_file = 'path_to_file';
+%[sig,fs] = wavread(input_file);
+
+%%%%% DEMO %%%%%
+fs = 44100;
+sig = rand(100000,1);
+%%%%% DEMO %%%%%
+
+% frame length
+N = 2^9;
+% DDM internal FFT size, for zero-padding, must be => N
+N_fft = N;
+% quantization, can be anything, but N_q = N makes the most sense
+N_q = N;
+% degree 2 DDM (linear FM, quadratic AM)
+ddm_dgr = 2;
+% DDM bandwidth (in bins)
+nr_bins = 5;
+% AM/FM model functions - 2nd degree polynomial
+dsc_idxs = [-(N-1)/2:(N-1)/2];
+t = dsc_idxs / fs;
+poly_model = poly_mat(t,ddm_dgr,N);
+
+% hann and hann' window
+w = 0.5+0.5*cos(2*pi*dsc_idxs/(N-1));
+wd = -pi/fs*(N-1)*sin(2*pi*dsc_idxs/(N-1));
+% example
+offset = 306; % some sample offset in the audio file
+bffr = sig(offset:offset+N-1); % get the buffer of lentgh N 
+
+
+%%% execute DDM on the whole frequency range (FFT) %%%
+% create linear systems
+[A, b, sig_fft] = ddm_lin_sys_fft(nr_bins, ddm_dgr, w', wd', poly_model', bffr, N, N_fft, fs);
+% fast pseudo-inverse
+ddm_estimates = flipud(pinv_2_fast(A,b,nr_bins,N-nr_bins+1).');
+% re-assign the spectral magnitude 
+reallocate(abs(sig_fft), ddm_estimates, N_q, fs);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pinv_2_fast.m	Thu Jul 25 13:07:07 2013 +0100
@@ -0,0 +1,19 @@
+function xx = pinv_2_fast(A,b,Q,N)
+
+b_mat = repmat(reshape(b,[Q,1,1,N]),[1 Q-1 2 1]);
+
+AA = zeros(Q,(Q-1),2,N);
+A1= zeros(Q,(Q-1),2,N);
+vct_idx = 1:Q;
+for q=1:Q
+	vct_Q = vct_idx([1:q-1,q+1:Q]);
+	A1(q,:,1,:) = shiftdim(A(vct_Q, 2,:),-1);
+	AA(q,:,1,:) = repmat(A(q,1,:),(Q-1),1) .* A(vct_Q,2,:);
+	
+	A1(q,:,2,:) = A(vct_idx([1:q-1,q+1:Q]), 1,:);
+	AA(q,:,2,:) = repmat(A(q,2,:),(Q-1),1) .* A(vct_Q,1,:);
+end
+A2 = (cat(3,conj(AA(:,:,1,:)-AA(:,:,2,:)) .* A1(:,:,1,:),conj(AA(:,:,1,:) - AA(:,:,2,:)) .* A1(:,:,2,:)));
+A3 = shiftdim(sum(sum(AA(:,:,2,:) .* conj(AA(:,:,1,:)-AA(:,:,2,:)))),3);
+A_nrm = [-A3, A3];
+xx = shiftdim(sum(sum(A2 .* b_mat)), 2).' ./ A_nrm;

--- a/reallocate.m	Thu Jul 25 10:40:55 2013 +0100
+++ b/reallocate.m	Thu Jul 25 13:07:07 2013 +0100
@@ -1,7 +1,8 @@
 function rl = reallocate(mag, ests, N_rs, fs)
-  frqs = imag(ests(end,:));
+  frqs = imag(ests(end,:))/2/pi;
   %frqs_rs = [0:N_rs-1]*fs/N_rs;
   bins_rs = round(frqs * N_rs/fs)-1;
+  bins_rs = bins_rs(bins_rs > 0 & bins_rs <= N_rs);
   rl = zeros(N_rs,1);
   % can this be done without for loop?
   for k=1:length(bins_rs)