gerard@1: % Copyright 2013 MUSIC TECHNOLOGY GROUP, UNIVERSITAT POMPEU FABRA
gerard@1: %
gerard@1: % Written by Gerard Roma <gerard.roma@upf.edu>
gerard@1: %
gerard@1: % This program is free software: you can redistribute it and/or modify
gerard@1: % it under the terms of the GNU Affero General Public License as published by
gerard@1: % the Free Software Foundation, either version 3 of the License, or
gerard@1: % (at your option) any later version.
gerard@1: %
gerard@1: % This program is distributed in the hope that it will be useful,
gerard@1: % but WITHOUT ANY WARRANTY; without even the implied warranty of
gerard@1: % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
gerard@1: % GNU Affero General Public License for more details.
gerard@1: %
gerard@1: % You should have received a copy of the GNU Affero General Public License
gerard@1: % along with this program.  If not, see <http://www.gnu.org/licenses/>.
gerard@1: 
gerard@1: function [ftr] = RQA(D, lmin,vmin)
gerard@1:   N = max(size(D));
gerard@1:   pts =  sum(D(:));
gerard@1: 
gerard@1:   det = 0;
gerard@1:   diags=[];
gerard@1:   for i = 1:N-1     
gerard@1:       line = [0;diag(D,i);0];
gerard@1:       starts=find(diff(line)==1); % line starts
gerard@1:       ends=find(diff(line)==-1); % line ends
gerard@1:       diags = [diags; ends-starts];
gerard@1:   end
gerard@1: 
gerard@1:   verts = [];
gerard@1:   intervals = [];
gerard@1:   for i = 1:N   
gerard@1:       line = [0;D(:,i);0];
gerard@1:       starts=find(diff(line)==1); % line starts
gerard@1:       ends=find(diff(line)==-1); % line ends
gerard@1:       verts = [verts;ends-starts];          
gerard@1:   end
gerard@1: 
gerard@1:   rr = pts/N^2; % recurrence rate
gerard@1:   det = sum(diags(diags>=lmin))/pts; % determinism
gerard@1:   adll = mean(diags(diags>=lmin)); % average diagonal line length
gerard@1:   sddll = std(diags(diags>=lmin)); % standard deviation
gerard@1:   
gerard@1:   dratio = N^2 * sum(diags(diags>=lmin))/(pts^2);
gerard@1:  
gerard@1:   lam = sum(verts(verts>=vmin))/pts; % laminarity
gerard@1:   tt = mean(verts(verts>=vmin)); % trappnig time
gerard@1:   sdvll = std(verts(verts>=lmin)); % standard deviation
gerard@1:   if(isnan(tt))
gerard@1:       tt = 0;
gerard@1:   end
gerard@1:   
gerard@1:   vratio = N^2 * sum(verts(verts>=vmin))/(pts^2); % ratio
gerard@1:  
gerard@1:   lmax = max(1,max(diags)); % maximum diagonal line length
gerard@1:   vmax = max(1,max(verts)); % maximum vertical line length
gerard@1:   
gerard@1:   if(isempty(lmax))
gerard@1:       lmax = 1;
gerard@1:   end
gerard@1:   
gerard@1:   
gerard@1:   if(isempty(vmax))
gerard@1:       vmax = 1;
gerard@1:   end
gerard@1:   
gerard@1:   ddiv = 1/lmax; % divergence
gerard@1:   vdiv = 1/vmax;
gerard@1:  
gerard@1:   diagsH = hist(diags,max(diags));
gerard@1:   diagsH = diagsH/sum(diagsH);
gerard@1:   z = find(diagsH==0);
gerard@1:   temp = diagsH.*log(diagsH);
gerard@1:   temp(z)=0; % remove NaNs from zero prob
gerard@1:   dentropy = sum(temp); % shannon entropy (diagonal)
gerard@1:  
gerard@1:  
gerard@1:   vertsH = hist(verts,max(verts));
gerard@1:   vertsH = vertsH/sum(vertsH);
gerard@1:   z = find(vertsH==0);
gerard@1:   temp = vertsH.*log(vertsH);
gerard@1:   temp(z)=0; % remove NaNs from zero prob
gerard@1:   ventropy = sum(temp); % shannon entropy (vertical)
gerard@1:  
gerard@1:   ftr = [rr,det,adll,lam,tt,dratio,ddiv,dentropy,vratio,ventropy,vdiv];
gerard@1:   end