forked from ossadtchi/PSIICOS
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathProjectAwayFromPowerComplete.m
More file actions
60 lines (53 loc) · 1.46 KB
/
ProjectAwayFromPowerComplete.m
File metadata and controls
60 lines (53 loc) · 1.46 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
function [CTp,Upwr,ds,nrmre,nrmim,nrmvc] = ProjectAwayFromPowerComplete( CT, G2dU, PwrRnk)
if(nargin<3)
PwrRnk = 350;
end;
if(nargin<2)
disp('use error\n');
CTp = [];
Upwr = [];
return;
end;
Ns = size(G2dU,2)/2; % two topography columns per each source of the grid
% perform projection of the coherence matrix away from the power only
Nch = size(G2dU,1);
% component
% project for each potential source
fprintf('Collecting power subspace of coherence span...\n');
A = zeros(Nch^2,Ns*3);
for i=1:Ns
gi = G2dU(:,2*i-1);
v = gi*gi';
A(:,3*i-2) = v(:)/norm(v(:));
gj = G2dU(:,2*i);
v = gj*gj';
A(:,3*i-1) = v(:)/norm(v(:));
v = gi*gj'+gj*gi';
A(:,3*i) = v(:)/norm(v(:));
end;
fprintf('Finding eigen space...\n');
AA = A*A';
[u s] = eigs(AA,PwrRnk);
ds = diag(s);
Upwr = u(:,1:PwrRnk);
if(size(CT,1)~=size(Upwr,1))
CTpvec = CT(:)-Upwr*(Upwr'*CT(:));
CTp = reshape(CTpvec,size(CT));
else
CTp = CT-Upwr*(Upwr'*CT);
end;
nrmin = [];
nrmre = [];
nrmvc = [];
if(nargout==6)
k =1;
for rnk = 1:10:size(u,2)
CTImp = imag(CT)-u(:,1:rnk)*(u(:,1:rnk)'*imag(CT));
CTRep = real(CT)-u(:,1:rnk)*(u(:,1:rnk)'*real(CT));
CTVCp = A-u(:,1:rnk)*(u(:,1:rnk)'*A);
nrmim(k) = norm(CTImp(:))/norm(imag(CT(:)));
nrmre(k) = norm(CTRep(:))/norm(real(CT(:)));
nrmvc(k) = norm(CTVCp(:))/norm(A(:));
k = k+1;
end;
end;