matthiasm@8: function df = get_dfs(x,p)
matthiasm@8: 
matthiasm@8: % function to calculate the following onset detection functions
matthiasm@8: % df{1} = complex spectral difference
matthiasm@8: % df{2} = spectral difference
matthiasm@8: % df{3} = phase deviation
matthiasm@8: x = mean(x,2);
matthiasm@8: 
matthiasm@8: % onset analysis step increment
matthiasm@8: o_step = p.winlen*2; % should be 1024
matthiasm@8: % onset analysis winlen
matthiasm@8: o_winlen = o_step*2; % should be 2048
matthiasm@8: 
matthiasm@8: hlfwin = o_winlen/2; % will be half fft size
matthiasm@8: 
matthiasm@8: % formulate hanningz window function
matthiasm@8: win = hanningz(o_winlen);
matthiasm@8: %  win = tukeywin(o_winlen,0.5);
matthiasm@8: 
matthiasm@8: % loop parameters
matthiasm@8: N = length(x);
matthiasm@8: pin = 0;
matthiasm@8: pend = N - o_winlen;
matthiasm@8: 
matthiasm@8: % vectors to store phase and magnitude calculations
matthiasm@8: theta1 = zeros(hlfwin,1);
matthiasm@8: theta2 = zeros(hlfwin,1);
matthiasm@8: oldmag = zeros(hlfwin,1);
matthiasm@8: 
matthiasm@8: % output onset detection functions
matthiasm@8: df = {};
matthiasm@8: 
matthiasm@8: % df sample number
matthiasm@8: k = 0;
matthiasm@8: while pin<pend
matthiasm@8:     
matthiasm@8:     k=k+1;
matthiasm@8:     % calculate windowed fft frame
matthiasm@8:     segment = x(pin+1:pin+o_winlen);
matthiasm@8:     X = fft(fftshift(win.*segment));
matthiasm@8:     %X = X.*fftshift(hanningz(length(X)));
matthiasm@8:     X = X(floor(length(X)/2)+1:length(X),:);
matthiasm@8:   %  X = X.*[length(X):-1:1]';
matthiasm@8:     % separate into magnitude
matthiasm@8:     mag = (abs(X));
matthiasm@8: 	
matthiasm@8: %  	mag = [0 abs(diff(log(eps+mag)))']';
matthiasm@8:     % and phase
matthiasm@8:     theta = angle(X);
matthiasm@8:     
matthiasm@8:     % complexsd part
matthiasm@8:         dev=princarg(theta-2*theta1+theta2);
matthiasm@8: %  		keyboard
matthiasm@8:         meas=oldmag - (mag.*exp(i.*dev));
matthiasm@8:         df{1}(k) = sum(sqrt((real(meas)).^2+(imag(meas)).^2));
matthiasm@8:         
matthiasm@8: %      % spectral difference part    
matthiasm@8: %          dev2=sqrt(mag.^2-oldmag.^2);
matthiasm@8: %          df{2}(k) = sum(abs(dev2(mag>0.1)));
matthiasm@8: %      
matthiasm@8: %      
matthiasm@8: %      % phase deviation part
matthiasm@8: %          df{3}(k) = sum(abs(dev((mag>0.1))));
matthiasm@8: %  %        df{3}(k) = sum(abs(dev(end-5:end)));
matthiasm@8: %          
matthiasm@8: %      % energy part
matthiasm@8: %          df{4}(k) = sum(abs(segment));
matthiasm@8: %          
matthiasm@8: %      % hfc part    
matthiasm@8: %          df{5}(k) = sum(mag.*(length(mag):-1:1)');
matthiasm@8: %      %    df{5}(k) = sum(mag.*rand(hlfwin,1));
matthiasm@8: %          
matthiasm@8:     % update vectors
matthiasm@8:     theta2 = theta1;
matthiasm@8:     theta1 = theta;
matthiasm@8:     oldmag = mag;
matthiasm@8:     pin = pin+o_step;
matthiasm@8: end
matthiasm@8: 
matthiasm@8: % now interpolate each detection function by a factor of 2,
matthiasm@8: % to get resolution of 11.6ms
matthiasm@8: 
matthiasm@8: for j=1:1,%5
matthiasm@8:     df{j} = interp(df{j},2);
matthiasm@8:     % also check there are no negative elements
matthiasm@8:     df{j}(df{j}<0) = 0;
matthiasm@8:     % and scale to sum to unity
matthiasm@8: %    df{j} = (df{j}+eps)/sum(df{j}+eps);
matthiasm@8: end
matthiasm@8: