diff --git a/fastlib/u/niche/functional/CAR.m b/fastlib/u/niche/functional/CAR.m
new file mode 100644
index 0000000000..502d0e8889
--- /dev/null
+++ b/fastlib/u/niche/functional/CAR.m
@@ -0,0 +1,13 @@
+function [filtered_data] = CAR(data)
+%CAR Apply CAR filter to EEG data
+%
+%            [filtered_data] = CAR(DATA)
+%            (time,numChannels,numEpochs)
+
+% Nishant Mehta 2007
+
+disp('applying CAR filter');
+
+for z = 1:28
+  filtered_data(:,z,:) = data(:,z,:) - (sum(data,2) / 28);
+end
diff --git a/fastlib/u/niche/functional/find_opt_subrotation.m b/fastlib/u/niche/functional/find_opt_subrotation.m
new file mode 100644
index 0000000000..e625c92ae1
--- /dev/null
+++ b/fastlib/u/niche/functional/find_opt_subrotation.m
@@ -0,0 +1,25 @@
+% find_opt_subrotation(X)
+function [theta_star, rotator_star] = find_opt_subrotation(X);
+
+thetas = linspace(-pi/4, pi/4, 200);
+
+max_h_sum = -Inf;
+
+for theta = thetas
+  rotator =[cos(theta) -sin(theta); sin(theta) cos(theta)];
+  rotated_X = rotator * X;
+  
+  h_sum = get_vasicek_entropy_estimate(rotated_X(1,:)) + ...
+	  get_vasicek_entropy_estimate(rotated_X(2,:));
+  
+  if h_sum > max_h_sum
+    max_h_sum = h_sum;
+    theta_star = theta;
+  end
+end
+
+rotator_star = ...
+    [cos(theta_star) -sin(theta_star); ...
+     sin(theta_star) cos(theta_star)];
+
+theta_star
diff --git a/fastlib/u/niche/functional/find_opt_unmixing_matrix.m b/fastlib/u/niche/functional/find_opt_unmixing_matrix.m
new file mode 100644
index 0000000000..4dd36488bc
--- /dev/null
+++ b/fastlib/u/niche/functional/find_opt_unmixing_matrix.m
@@ -0,0 +1,47 @@
+% find_opt_unmixing_matrix(X) - find the unmixing matrix of W that
+% best separates the components of X
+% X is a p x n matrix
+function [Y_pos,Y_neg,W_pos,W_neg] = find_opt_unmixing_matrix(X);
+
+p = size(X, 1);
+
+current_X = X;
+
+total_rotator = eye(p);
+
+for epoch = 1:1
+
+  for i = 1:p
+    for j = i+1:p
+      
+      disp(sprintf('[i,j] = [%d,%d]', i, j));
+      
+      
+      subspace = [current_X(i,:) ; current_X(j,:)];
+      [theta_star, rotator_star] = find_opt_subrotation(subspace);
+      
+      new_rotator=eye(p);
+      new_rotator(i,i)=cos(theta_star);
+      new_rotator(i,j)=-sin(theta_star);
+      new_rotator(j,i)=sin(theta_star);
+      new_rotator(j,j)=cos(theta_star);
+      
+      total_rotator = new_rotator * total_rotator;
+      current_X = total_rotator * X;
+      
+    end  
+  end
+end
+
+Y_pos = current_X;
+
+neg_total_rotator=eye(p);
+neg_total_rotator(1,1)=cos(theta_star + pi);
+neg_total_rotator(1,2)=-sin(theta_star + pi);
+neg_total_rotator(2,1)=sin(theta_star + pi);
+neg_total_rotator(2,2)=cos(theta_star + pi);
+
+Y_neg = neg_total_rotator * X;
+W_pos = total_rotator;
+W_neg = neg_total_rotator;
+
diff --git a/fastlib/u/niche/functional/fun.m b/fastlib/u/niche/functional/fun.m
index 797dbc864b..b4f58c9610 100644
--- a/fastlib/u/niche/functional/fun.m
+++ b/fastlib/u/niche/functional/fun.m
@@ -15,8 +15,11 @@ rand('state', sum(100*clock))
 
 % sample from laplacian
 
+clear;
+
 D = 2;
 N = 10000;
+p = 30;
 
 mu = 0;
 sigma = 1;
@@ -58,22 +61,22 @@ l_x = l_x - repmat(mean(l_x')', 1, N);
 
 % generate b-spline basis curves
 t = linspace(0,1,1000);
-mybasis = create_bspline_basis([0 1], 30, 4);
+mybasis = create_bspline_basis([0 1], p, 4);
 basis_curves = eval_basis(t, mybasis);
 
 
-load e1e2;
-e = [e1(t); e2(t)]';
+load s1s2;
+s = [s1(t); s2(t)]';
 
 
-data = e * l_x;
+data = s * l_x;
 
 
 
 myfd_data = data2fd(data, t, mybasis);
 coef = getcoef(myfd_data);
 %data1 = basis_curves * coef(:,1);
-pca_results = pca_fd(myfd_data, 30);
+pca_results = pca_fd(myfd_data, p);
 pc_coef = getcoef(pca_results.harmfd);
 pc_curves = basis_curves * pc_coef;
 pc_scores = pca_results.harmscr;
@@ -83,31 +86,75 @@ pc_scores = pca_results.harmscr;
 
 
 % encode our source functions e1 and e2 using the pc basis
-for i=1:30
-  e1_weights(i) = ...
-      diff(ppval(fnint(spline(t, e1(t) .* pc_curves(:,i)')), ...
+for i=1:p
+  s1_weights(i) = ...
+      diff(ppval(fnint(spline(t, s1(t) .* pc_curves(:,i)')), ...
 		 [0 1]));
-  e2_weights(i) = ...
-      diff(ppval(fnint(spline(t, e2(t) .* pc_curves(:,i)')), ...
+  s2_weights(i) = ...
+      diff(ppval(fnint(spline(t, s2(t) .* pc_curves(:,i)')), ...
 		 [0 1]));
 end
 
 for i=1:N
-  e1_scores(i) = dot(e1_weights, pc_scores(i,:));
-  e2_scores(i) = dot(e2_weights, pc_scores(i,:));
+  s1_scores(i) = dot(s1_weights, pc_scores(i,:));
+  s2_scores(i) = dot(s2_weights, pc_scores(i,:));
 end
 
 
+p_small = 2;
+
+sub_pc_coef = pc_coef(:,1:p_small);
+E = pc_scores(:,1:p_small)';
+
+[Y_pos,Y_neg,W_pos,W_neg] = find_opt_unmixing_matrix(E);
+
+
+for i=1:p_small
+  h_E(i) = get_vasicek_entropy_estimate(E(i,:));
+  h_Y_pos(i) = get_vasicek_entropy_estimate(Y_pos(i,:));
+  h_Y_neg(i) = get_vasicek_entropy_estimate(Y_neg(i,:));
+end
+
+
+
+ic_coef_pos = (W_pos * sub_pc_coef')';
+ic_coef_neg = (W_neg * sub_pc_coef')';
+
+
+sub_pc_curves = basis_curves * sub_pc_coef;
+ic_curves_pos = basis_curves * ic_coef_pos;
+ic_curves_neg = basis_curves * ic_coef_neg;
+
+
+figure(1);
+clf;
+hold on;
+plot(s, 'b');
+plot(sub_pc_curves, 'r');
+plot(ic_curves_pos, 'g');
+plot(ic_curves_neg, 'c');
+
+
+
+% using sub_pc_coef', recover the data
+
+
+  
+
+
+
+% now we want to find a matrix W that unmixes well
+
+
 
 
 % let f be some candidate solution
 
-let the data be encoded as your mother
-f1_weights = rand(30,1);
-f1_weights = f1_weights / norm(f1_weights);
+%f1_weights = rand(p,1);
+%f1_weights = f1_weights / norm(f1_weights);
 
 
-f1 = pc_curves * f1_weights;
+%f1 = pc_curves * f1_weights;
 
 
 
@@ -121,11 +168,15 @@ f1 = pc_curves * f1_weights;
 %    curves with splines, then use quadrature to evaluate the
 %    integral in [0,1]
 
-f1_scores = zeros(N,1);
+%f1_scores = zeros(N,1);
+
+%for i=1:N
+%  f1_scores(i) = dot(f1_weights, pc_scores(i,:));
+%end
+
+
+
 
-for i=1:N
-  f1_scores(i) = dot(f1_weights, pc_scores(i,:));
-end
 
 
 
diff --git a/fastlib/u/niche/functional/funica.m b/fastlib/u/niche/functional/funica.m
new file mode 100644
index 0000000000..b4f58c9610
--- /dev/null
+++ b/fastlib/u/niche/functional/funica.m
@@ -0,0 +1,205 @@
+% Given
+%   distribution F
+%   observations X (dims N x T)
+
+% construct random vector w by sampling from F
+
+% sort w descending
+
+% attribute highest values of w to observations with maximal variance
+
+
+% initialize random number generator
+rand('state', sum(100*clock))
+
+
+% sample from laplacian
+
+clear;
+
+D = 2;
+N = 10000;
+p = 30;
+
+mu = 0;
+sigma = 1;
+b = sigma/2;
+
+
+
+% laplacian
+clear l_x l_px;
+for i=1:D
+  l_x(i,:) = laplacinv(rand(N, 1), mu, b);
+  l_px(i,:) = laplacpdf(l_x(i,:), mu, b);
+end
+
+% center the sampling distribution
+l_x = l_x - repmat(mean(l_x')', 1, N);
+
+
+% gaussian
+%clear g_x g_px;
+%for i=1:D
+%  g_x(i,:) = norminv(rand(N, 1), mu, sigma);
+%  g_px(i,:) = normpdf(g_x(i,:), mu, sigma);
+%end
+
+% center the sampling distribution
+%g_x = g_x - repmat(mean(g_x')', 1, N);
+
+
+
+% mix the source signals
+
+%a = rand(D,D);
+%X = a * l_x;
+
+%shuffling = shuffle(1:N);
+%X_shuffled = X(:,shuffling);
+
+
+% generate b-spline basis curves
+t = linspace(0,1,1000);
+mybasis = create_bspline_basis([0 1], p, 4);
+basis_curves = eval_basis(t, mybasis);
+
+
+load s1s2;
+s = [s1(t); s2(t)]';
+
+
+data = s * l_x;
+
+
+
+myfd_data = data2fd(data, t, mybasis);
+coef = getcoef(myfd_data);
+%data1 = basis_curves * coef(:,1);
+pca_results = pca_fd(myfd_data, p);
+pc_coef = getcoef(pca_results.harmfd);
+pc_curves = basis_curves * pc_coef;
+pc_scores = pca_results.harmscr;
+
+
+
+
+
+% encode our source functions e1 and e2 using the pc basis
+for i=1:p
+  s1_weights(i) = ...
+      diff(ppval(fnint(spline(t, s1(t) .* pc_curves(:,i)')), ...
+		 [0 1]));
+  s2_weights(i) = ...
+      diff(ppval(fnint(spline(t, s2(t) .* pc_curves(:,i)')), ...
+		 [0 1]));
+end
+
+for i=1:N
+  s1_scores(i) = dot(s1_weights, pc_scores(i,:));
+  s2_scores(i) = dot(s2_weights, pc_scores(i,:));
+end
+
+
+p_small = 2;
+
+sub_pc_coef = pc_coef(:,1:p_small);
+E = pc_scores(:,1:p_small)';
+
+[Y_pos,Y_neg,W_pos,W_neg] = find_opt_unmixing_matrix(E);
+
+
+for i=1:p_small
+  h_E(i) = get_vasicek_entropy_estimate(E(i,:));
+  h_Y_pos(i) = get_vasicek_entropy_estimate(Y_pos(i,:));
+  h_Y_neg(i) = get_vasicek_entropy_estimate(Y_neg(i,:));
+end
+
+
+
+ic_coef_pos = (W_pos * sub_pc_coef')';
+ic_coef_neg = (W_neg * sub_pc_coef')';
+
+
+sub_pc_curves = basis_curves * sub_pc_coef;
+ic_curves_pos = basis_curves * ic_coef_pos;
+ic_curves_neg = basis_curves * ic_coef_neg;
+
+
+figure(1);
+clf;
+hold on;
+plot(s, 'b');
+plot(sub_pc_curves, 'r');
+plot(ic_curves_pos, 'g');
+plot(ic_curves_neg, 'c');
+
+
+
+% using sub_pc_coef', recover the data
+
+
+  
+
+
+
+% now we want to find a matrix W that unmixes well
+
+
+
+
+% let f be some candidate solution
+
+%f1_weights = rand(p,1);
+%f1_weights = f1_weights / norm(f1_weights);
+
+
+%f1 = pc_curves * f1_weights;
+
+
+
+
+
+% we evaluate some f by considering projections P of the data onto f
+
+% define the l2 norm for functional space:
+%    given some vector a and another vector b, we dot multiply the
+%    two vectors at the specified values, then approximate the
+%    curves with splines, then use quadrature to evaluate the
+%    integral in [0,1]
+
+%f1_scores = zeros(N,1);
+
+%for i=1:N
+%  f1_scores(i) = dot(f1_weights, pc_scores(i,:));
+%end
+
+
+
+
+
+
+
+%given the f1_scores, what to do now?
+
+
+
+
+% objective function
+% min sigma H(X_i)
+%  X    i
+%for a given input variable X, we seek to minimize the sum of the ...
+%      entropies of the marginal distributions we consider the sum
+%      of the entropies of the marginal distributions
+
+% in the case of one dimension, we are given a set of scalar values
+% - we can study the distribution of these values
+
+% in the case of two dimensions, we are given a set of 2-vector
+% values
+% we want to know the entropy of this distribution
+
+
+% the m spacing estimator studies the spacing between the sample
+% points
+
diff --git a/fastlib/u/niche/functional/get_vasicek_entropy_estimate.m b/fastlib/u/niche/functional/get_vasicek_entropy_estimate.m
new file mode 100644
index 0000000000..dc2bc25c8e
--- /dev/null
+++ b/fastlib/u/niche/functional/get_vasicek_entropy_estimate.m
@@ -0,0 +1,23 @@
+% get_vasicek_entropy_estimate() - estimate entropy using vasicek's
+% entropy estimator
+function h = get_vasicek_entropy_estimate(X);
+
+n = length(X);
+m = round(sqrt(n)/2);
+Z = sort(X, 'ascend');
+
+sum_logs = 0;
+
+for i = 1:n
+  if (i + m) > n
+    sum_logs = sum_logs + log(Z(n) - Z(i-m));
+  elseif (i-m) < 1
+    sum_logs = sum_logs + log(Z(i+m) - Z(1));
+  else    
+    sum_logs = sum_logs + log(Z(i+m) - Z(i-m));
+  end
+end
+
+  
+
+h = (sum_logs/n) + log(n/(2*m));
diff --git a/fastlib/u/niche/functional/l2_fnorm.m b/fastlib/u/niche/functional/l2_fnorm.m
index 2e733a2a16..cda16e054b 100644
--- a/fastlib/u/niche/functional/l2_fnorm.m
+++ b/fastlib/u/niche/functional/l2_fnorm.m
@@ -5,4 +5,4 @@ if size(f) == size(g')
   g = g';
 end
 
-norm = quad(@ppval, domain(1), domain(end), [], [], spline(domain, f .* g));
+norm = diff(ppval(fnint(spline(domain, f .* g)), [domain(1) domain(end)]));
\ No newline at end of file
diff --git a/fastlib/u/niche/functional/test.m b/fastlib/u/niche/functional/test.m
index 76637bd752..7cfb4d9500 100644
--- a/fastlib/u/niche/functional/test.m
+++ b/fastlib/u/niche/functional/test.m
@@ -5,7 +5,9 @@ c4 = 18;
 
 mybasis = create_bspline_basis([0 .5], 30, 4);
 
-load ../bci_comp_2003_dataSet_IV/sp1s_aa_1000Hz;
+load ~/brains/bci_comp_2003_dataSet_IV/sp1s_aa_1000Hz;
+
+
 
 x_train = CAR(x_train);
 ground_truth = y_train;
@@ -71,4 +73,4 @@ right_correct = ...
     abs(score_test(rights_test, pc) - lmean);
 
 [sum(left_correct) / length(left_correct)
- sum(right_correct) / length(right_correct)]
\ No newline at end of file
+ sum(right_correct) / length(right_correct)]