--- a
+++ b/preprocessOfApneaECG/BioSigKit/Algorithms/projective.m
@@ -0,0 +1,382 @@
+function [px] = projective (x, t, d, m, r) 
+% PROJECTIVE  
+%
+% CALL: [px] = projective (x, t, d, m, r);
+%
+% INPUTS:
+%   x: the vector containing the scalar observations
+%   t: delay parameter used in embedding
+%   d: embedding dimension to consider
+%   m: dimensiona of null space
+%   r: number of nearest neighbors to consider
+%
+% OUTPUT:
+%   px: cleaned scalar observations
+%
+% copiright (c) and written by David Chelidze 4/20/2011
+% [px] = projective (foetal_ecg(:,2), 1, 50, 45, 1500);
+ 
+if nargin < 5
+    error('need five input variables')
+end
+if nargin > 5
+    error('too many input variables, needs only five')
+end
+ 
+% check data and make a raw vector of scalar data
+if min(size(x)) ~= 1
+   error('first input should be a vector of scalar data points.')
+else
+    x=x(:)'; % just the way we need it
+end
+px = x;
+ 
+n = length(x) - d*t; % maximum number of reconstructed points possible
+
+indx = 1:n; % index of points in the embedding
+y = zeros(d+1,n);
+dy = y; % initialize corrective array
+
+% reconstruct the phase space
+for i = 1:d+1
+    y(i,:) = x(indx+(i-1)*t);
+end
+
+% define the weighting matrix
+W = eye(d+1,d+1); W(1,1) = 10^3; W(d+1,d+1) = 10^3;
+
+fprintf('Partitioning data\n')
+[kd_tree, q] = kd_partition(y, 256); % partition into kd-tree
+
+fprintf('Filtering data\n')
+for k = 1:n % find the nearest neighbors and their properties
+    [nns, ~, ~] = kd_search(y(:,k), r, kd_tree, y(:,q));
+    nns = nns{1};
+    for i = 1:d
+        nns(i,:) = nns(i,:) - mean(nns(i,:)); % remove the mean
+    end
+    D = nns*nns.'/r; % find correlation matrix
+    G = W*D*W; % transform it
+    %[E,~] = eig(G); % assuming the MATLAB gives them in increasing order
+    %E = E(:,end-1:end); 
+    %E = E(1:m);
+    % if eig does not sort correctly use the following instead:
+    [E,g] = eig(G);
+    [~,srt] = sort(diag(g));
+    E = E(:,srt(1:m));
+    nns(:,1) = (nns(:,1) - mean(nns,2));
+    dy(:,k) = W\(E*E.'*W*nns(:,1)); % corrective vector
+    dy(:,k) = mean(nns,2) + dy(:,k);
+    px(i) = mean(dy(:,k));
+end
+
+% clean the time series
+indx = [0:t:d*t, n+(0:t:d*t)];
+for q = [1:d+1, d:-1:1]
+    for l = (1+indx(q)):indx(q+1)
+        cx = 0;
+        for k = 1:q % go through all of possible corrections
+        cx = cx + dy(k,l-(k-1)*t); % all possible corrections to x(l)
+        end
+        px(l) = x(l) -  cx/q; % cleaned values
+    end
+end
+
+%--------------------------------------------------------------------------            
+function [kd_tree, r] = kd_partition(y, b, c)
+% KD_PARTITION  Create a kd-tree and partitioned database for efficiently 
+%               finding the nearest neighbors to a point in a 
+%               d-dimensional space.
+%
+% USE: [kd_tree, r] = kd_partition(y, b, c);
+%
+% INPUT:
+%   y: original multivariate data (points arranged columnwise, size(y,1)=d)
+%   b: maximum number of distinct points for each bin (default is 100)
+%   c: minimum range of point cluster within a final leaf (default is 0)
+%
+% OUTPUT:
+%   kd_tree structure with the following variables:
+%       splitdim: dimension used in splitting the node
+%       splitval: corresponding cutting point
+%       first & last: indices of points in the node
+%       left & right: node #s of consequent branches to the current node
+%   r: sorted index of points in the original y corresponding to the leafs
+%
+% to find k-nearest neighbors use kd_search.m
+%
+% copyrighted (c) and written by David Chelidze, January 28, 2009.
+ 
+% check the inputs
+if nargin == 0
+    error('Need to input at least the data to partition')
+elseif nargin > 3
+    error('Too many inputs')
+end
+ 
+ 
+% initializes default variables if needed
+if nargin < 2
+    b = 100;
+end
+if nargin < 3
+    c = 0;
+end
+ 
+[d, n] = size(y); % get the dimension and the number of points in y
+ 
+r = 1:n; % initialize original index of points in y
+ 
+% initializes variables for the number of nodes and the last node
+node = 1;
+last = 1;
+ 
+% initializes the first node's cut dimension and value in the kd_tree
+kd_tree.splitdim = 0;
+kd_tree.splitval = 0;
+ 
+% initializes the bounds on the index of all points
+kd_tree.first = 1;
+kd_tree.last = n;
+ 
+% initializes location of consequent branches in the kd_tree
+kd_tree.left = 0;
+kd_tree.right = 0;
+
+while node <= last % do until the tree is complete
+    
+    % specify the index of all the points that are partitioned in this node
+    segment = kd_tree.first(node):kd_tree.last(node);
+    
+    % determines range of data in each dimension and sorts it
+    [rng, index] = sort(range(y(:,segment),2));
+    
+    % now determine if this segment needs splitting (cutting)
+    if rng(d) > c && length(segment)>= b % then split
+        yt = y(:,segment); 
+        rt = r(segment);
+        [sorted, sorted_index] = sort(yt(index(d),:));
+        % estimate where the cut should go
+        lng = size(yt,2);
+        cut = (sorted(ceil(lng/2)) + sorted(floor(lng/2+1)))/2;
+        L = (sorted <= cut); % points to the left of cut
+        if sum(L) == lng % right node is empty
+            L = (sorted < cut); % decrease points on the left
+            cut = (cut + max(sorted(L)))/2; % adjust the cut
+        end
+        
+        % adjust the order of the data in this node
+        y(:,segment) = yt(:,sorted_index); 
+        r(segment) = rt(sorted_index);
+ 
+        % assign appropriate split dimension and split value
+        kd_tree.splitdim(node) = index(d);
+        kd_tree.splitval(node) = cut;
+        
+        % assign the location of consequent bins and 
+        kd_tree.left(node) = last + 1;
+        kd_tree.right(node) = last + 2;
+        
+        % specify which is the last node at this moment
+        last = last + 2;
+        
+        % initialize next nodes cut dimension and value in the kd_tree
+        % assuming they are terminal at this point
+        kd_tree.splitdim = [kd_tree.splitdim 0 0];
+        kd_tree.splitval = [kd_tree.splitval 0 0];
+ 
+        % initialize the bounds on the index of the next nodes
+        kd_tree.first = [kd_tree.first segment(1) segment(1)+sum(L)];
+        kd_tree.last = [kd_tree.last segment(1)+sum(L)-1 segment(lng)];
+ 
+        % initialize location of consequent branches in the kd_tree
+        % assuming they are terminal at this point
+        kd_tree.left = [kd_tree.left 0 0];
+        kd_tree.right = [kd_tree.right 0 0];
+        
+    end % the splitting process
+ 
+    % increment the node
+    node = node + 1;
+ 
+end % the partitioning
+
+%--------------------------------------------------------------------------            
+function [pqr, pqd, pqi] = kd_search(y,r,tree,yp)
+% KD_SEARCH     search kd_tree for r nearest neighbors of point yq.
+%               Need to partition original data using kd_partition.m
+%
+% USE: [pqr, pqd, pqi] = kd_search(y,r,tree,yp);
+%
+% INPUT:
+%       y: array of query points in columnwise form
+%       r: requested number of nearest neighbors to the query point yq
+%       tree: kd_tree constructed using kd_partition.m
+%       yp: partitioned (ordered) set of data that needs to be searched
+%          (using my kd_partirion you want to input ym(:,indx), where ym 
+%           is the data used for partitioning and indx sorted index of ym)
+%
+% OUTPUT:
+%       pqr: cell array of the r nearest neighbors of y in yp 
+%       pqd: cell array of the corresponding distances
+%       pqi: cell array of the indices of r nearest neighbors of y in yp 
+%
+% copyright (c) and written by David Chelidze, February 02, 2009.
+ 
+% check inputs
+if nargin < 4
+    error('Need four input variables to work')
+end
+ 
+% declare global variables for subfunctions
+global yq qri qr qrd finish b_lower b_upper
+ 
+[~, n] = size(y);
+pqr = cell(n,1);
+pqd = cell(n,1);
+pqi = cell(n,1);
+ 
+for k = 1:n
+    yq = y(:,k);
+    qrd = Inf*ones(1,r); % initialize array for r distances
+    qr = zeros(size(yq,1),r); % initialize r nearest neighbor points
+    qri = zeros(1,r); % initialize index of r nearest neighbors
+    finish = 0; % becomes 1 after search is complete
+ 
+    % set up the box bounds, which start at +/- infinity (whole kd space)
+    b_upper = Inf*ones(size(yq));
+    b_lower = -b_upper;
+ 
+    kdsrch(1,r,tree,yp); % start the search from the first node
+    pqr{k} = qr;
+    pqd{k} = sqrt(qrd);
+    pqi{k} = qri;
+end
+ 
+%--------------------------------------------------------------------------            
+function kdsrch(node,r,tree,yp)
+% KDSRCH    actual kd search
+% this drills down the tree to the end node and updates the 
+% nearest neighbors list with new points
+%
+% INPUT: starting node number, and kd_partition data
+%
+% copyright (c) and written by David Chelidze, February 02, 2009.
+ 
+global yq qri qr qrd finish b_lower b_upper
+ 
+% first find the terminal node containing yq
+if tree.left(node) ~= 0 % not a terminal node, search deeper
+ 
+    % first determin which child node to search
+    if yq(tree.splitdim(node)) <= tree.splitval(node)
+        % need to search left child node
+        tmp = b_upper(tree.splitdim(node));
+        b_upper(tree.splitdim(node)) = tree.splitval(node);
+        kdsrch(tree.left(node),r,tree,yp);
+        b_upper(tree.splitdim(node)) = tmp;
+    else % need to search the right child node
+        tmp = b_lower(tree.splitdim(node));
+        b_lower(tree.splitdim(node)) = tree.splitval(node);
+        kdsrch(tree.right(node),r,tree,yp);
+        b_lower(tree.splitdim(node)) = tmp;
+    end % when the terminal node (leaf) containing yq is reached
+    if finish % done searching
+        return
+    end
+ 
+    % check if other nodes need to be searched
+    if yq(tree.splitdim(node)) <= tree.splitval(node)
+        tmp = b_lower(tree.splitdim(node));
+        b_lower(tree.splitdim(node)) = tree.splitval(node);
+        if overlap(yq, b_lower, b_upper, qrd(r)) 
+            % need to search the right child node
+            kdsrch(tree.right(node),r,tree,yp);
+        end
+        b_lower(tree.splitdim(node)) = tmp;
+    else
+        tmp = b_upper(tree.splitdim(node));
+        b_upper(tree.splitdim(node)) = tree.splitval(node);
+        if overlap(yq, b_lower, b_upper, qrd(r)) 
+            % need to search the left child node
+            kdsrch(tree.left(node),r,tree,yp);
+        end
+        b_upper(tree.splitdim(node)) = tmp;
+    end % when all the other nodes are searched
+    if finish % done searching
+        return
+    end
+ 
+else % this is a terminal node: update the nearest neighbors
+ 
+    yt = yp(:,tree.first(node):tree.last(node)); % all points in node
+    dstnc = zeros(1,size(yt,2));
+    for k = 1:size(yt,1)
+        dstnc = dstnc + (yt(k,:)-yq(k)).^2; 
+    end
+    %dstnc = sum((yt-yq*ones(1,size(yt,2))).^2,1); % distances squared
+    qrd = [qrd dstnc]; % current list of distances squared
+    qr = [qr yt]; % current list of nearest neighbors
+    qri = [qri tree.first(node):tree.last(node)];
+    [qrd, indx] = sort(qrd); % sorted distances squared and their index
+    qr = qr(:,indx); % sorted list of nearest neighbors
+    qri = qri(indx); % sorted list of indexes
+    if size(qr,2) > r % truncate to the first r points
+        qrd = qrd(1:r);
+        qr = qr(:,1:r);
+        qri = qri(1:r);
+    end
+    % be done if all points are with this box on the first run
+    if within(yq, b_lower, b_upper, qrd(r));
+        finish = 1;
+    end % otherwise (during backtracking) WITHIN will always return 0.
+    
+end % if
+ 
+%--------------------------------------------------------------------------            
+function flag = within(yq, b_lower, b_upper, ball)
+% WITHIN    check if additional nodes need to be searched (i.e. if the ball
+%  centered at yq and containing all current nearest neighbors overlaps the
+%  boundary of the leaf box containing yq)
+%
+% INPUT:
+%   yq: query point
+%   b_lower, b_upper: lower and upper bounds on the leaf box
+%   ball: square of the radius of the ball centered at yq and containing
+%         all current r nearest neighbors
+% OUTPUT:
+%   flag: 1 if ball does not intersect the box, 0 if it does
+%
+% Modified by David Chelidze on 02/03/2009.
+ 
+if ball <= min([abs(yq-b_lower)', abs(yq-b_upper)'])^2
+    % ball containing all the current nn is inside the leaf box (finish)
+    flag = 1;
+else % ball overlaps other leaf boxes (continue recursive search)
+    flag = 0; 
+end
+ 
+%--------------------------------------------------------------------------            
+function flag = overlap(yq, b_lower, b_upper, ball)
+% OVERLAP   check if the current box overlaps with the ball centered at yq
+%   and containing all current r nearest neighbors’.
+%
+% INPUT:
+%   yq: query point
+%   b_lower, b_upper: lower and upper bounds on the current box
+%   ball: square of the radius of the ball centered at yq and containing
+%         all current r nearest neighbors
+% OUTPUT:
+%   flag: 0 if ball does not overlap the box, 1 if it does
+%
+% Modified by David Chelidze on 02/03/2009.
+ 
+il = find(yq < b_lower); % index of yq coordinates that are lower the box 
+iu = find(yq > b_upper); % index of yq coordinates that are upper the box
+% distance squared from yq to the edge of the box
+dst = sum((yq(il)-b_lower(il)).^2,1)+sum((yq(iu)-b_upper(iu)).^2,1);
+if dst >= ball % there is no overlap (finish this search)    
+    flag = 0;
+else % there is overlap and the box needs to be searched for nn
+    flag = 1;
+end
\ No newline at end of file