BML_VARAVGREF_COEF calculates the crosstalk coefficients between electrodes Use as [coef, fitted_COV] = bml_varavgref_coef(COV) Calculates crosstal coefficients as defined by electrical crosstalk model Arguments COV - variance-covariance matrix between channels Returns coef - vector of crosstalk coefficients 0 <= coef <= 1 fitted_COV - fitted variance-covariance matrix
0001 function [coef, fitted_COV] = bml_varavgref_coef(COV) 0002 0003 % BML_VARAVGREF_COEF calculates the crosstalk coefficients between electrodes 0004 % 0005 % Use as 0006 % [coef, fitted_COV] = bml_varavgref_coef(COV) 0007 % 0008 % Calculates crosstal coefficients as defined by electrical crosstalk model 0009 % 0010 % Arguments 0011 % COV - variance-covariance matrix between channels 0012 % 0013 % Returns 0014 % coef - vector of crosstalk coefficients 0 <= coef <= 1 0015 % fitted_COV - fitted variance-covariance matrix 0016 0017 0018 assert(isnumeric(COV)); 0019 assert(size(COV,1)==size(COV,2)); 0020 0021 n = size(COV,1); 0022 0023 %calculate model values from crosstalk parameters 0024 function [V0V0,sigma2_S,VkVj,c] = varavgref_model(p) 0025 0026 %p=zeros(1,n) 0027 VkVj=diag(diag(COV)); % matrix of predicted variances-covariances 0028 sigma2_S = zeros(1,n); %initializing singal variance 0029 c = exp(p)./(exp(p)+1); % vector of crosstalk coefficients 0030 sum_c = sum(c); 0031 0032 %calculating V0 variance 0033 V0V0=0; 0034 for l=1:n 0035 for j=1:n 0036 V0V0 = V0V0 + exp(p(l)+p(j))*COV(l,j); 0037 end 0038 end 0039 V0V0 = V0V0*(sum(exp(p))^(-2)); 0040 0041 %calculating signal variances sigma2_S 0042 for k=1:n 0043 sigma2_S(k) = (COV(k,k)-(c(k)^2)*V0V0)/(((1-c(k))^2)+2*(c(k)^2)*(1-c(k))/sum_c); 0044 end 0045 0046 %completing predicted variance-covariance matrix 0047 for k=1:(n-1) 0048 for j=(k+1):n 0049 VkVj(k,j)= c(k)*c(j)*( (1-c(k))*sigma2_S(k)/sum_c + (1-c(j))*sigma2_S(j)/sum_c + V0V0); 0050 VkVj(j,k)=VkVj(k,j); 0051 end 0052 end 0053 end 0054 0055 % optimization routine 0056 p = zeros(1,n); 0057 [f_V0V0,f_sigma2_S,f_VkVj,f_c] = varavgref_model(p); 0058 0059 %first guess based on solving approximated cuadratic expression 0060 p1 = zeros(1,n); 0061 for i=1:n 0062 a_= -f_sigma2_S(i); 0063 b_= f_sigma2_S(i); 0064 b_= b_ + (sum(f_c .* (1-f_c) .* f_sigma2_S)-f_c(i)*(1-f_c(i))*f_sigma2_S(i))/sum(f_c); 0065 b_= b_ + f_V0V0*sum(f_c); 0066 c_= -sum(COV(i,:))+COV(i,i); 0067 0068 d_ = b_^2-4*a_*c_; 0069 if d_ >= 0 0070 c_k1=(-b_ + sqrt(d_))/(2*a_); 0071 c_k2=(-b_ - sqrt(d_))/(2*a_); 0072 if 0<=c_k1 && c_k1<=1 0073 c_k=c_k1; 0074 elseif 0<=c_k2 && c_k2<=1 0075 c_k=c_k2; 0076 else 0077 c_k=(-b_/(2*a_)); 0078 end 0079 else 0080 c_k=(-b_/(2*a_)); 0081 end 0082 0083 if c_k <=0 0084 p1(i) = -10; 0085 elseif c_k >= 1 0086 p1(i) = 10; 0087 else 0088 p1(i) = log(abs(c_k/(1-c_k))); 0089 end 0090 end 0091 0092 %cost function 0093 function [cost, grad] = costfun(p) 0094 [~,~,VkVj,c] = varavgref_model(p); 0095 0096 %calculating cost in variance metric 0097 cost=0; 0098 for k=1:n 0099 for j=1:n 0100 cost = cost + (VkVj(k,j)-COV(k,j))^2; 0101 end 0102 end 0103 cost = cost / sum(sum(COV.^2)); 0104 0105 %ridge regression on crosstalk coefficients 0106 cost = cost + (1e-3)*sum(c.^2); 0107 end 0108 0109 %optimizing from first guess 0110 p=fminunc(@costfun,p1,optimoptions('fminunc','MaxFunctionEvaluations',1e6,'OptimalityTolerance',1e-6)); 0111 0112 % p=p1; 0113 % 0114 % %cost function 0115 % function [cost, grad] = costfun2(p) 0116 % [~,~,VkVj,~] = varavgref_model(p); 0117 % 0118 % %calculating cost in correlation metric 0119 % 0120 % %Transformin VkVj to positive semi definte 0121 % VkVj(ismissing(VkVj))=0; 0122 % [V,D] = eig(VkVj); 0123 % D(D<0)=0; 0124 % psd_VkVj = real(V*D/(V)); 0125 % cor_VkVj = corrcov(psd_VkVj); 0126 % cor_COV = corrcov(COV); 0127 % 0128 % cost=0; 0129 % for k=1:n 0130 % for j=1:n 0131 % cost = cost + (cor_VkVj(k,j)-cor_COV(k,j))^2; 0132 % %cost = cost + (clip_inv_sigmoid(cor_VkVj(k,j))-clip_inv_sigmoid(cor_COV(k,j)))^2; 0133 % end 0134 % end 0135 % %cost = cost / sum(sum(COV.^2)); 0136 % end 0137 % 0138 % p=fminunc(@costfun2,p,optimoptions('fminunc','MaxFunctionEvaluations',1e9,'OptimalityTolerance',1e-9)); 0139 % 0140 % function x = clip_inv_sigmoid(y) 0141 % if y > 0.999954602131298 % exp(10)/(exp(10)+1) 0142 % x = 10; 0143 % elseif y < 4.539786870243439e-05 % exp(-10)/(exp(-10)+1) 0144 % x = -10; 0145 % else 0146 % x = log(y/(1-y)); 0147 % end 0148 % end 0149 % 0150 % %cost function 0151 % function [cost, grad] = costfun3(p) 0152 % [~,~,VkVj,~] = varavgref_model(p); 0153 % 0154 % %calculating cost in correlation metric 0155 % 0156 % %Transformin VkVj to positive semi definte 0157 % VkVj(ismissing(VkVj))=0; 0158 % [V,D] = eig(VkVj); 0159 % D(D<0)=0; 0160 % psd_VkVj = real(V*D/(V)); 0161 % cor_VkVj = corrcov(psd_VkVj); 0162 % cor_COV = corrcov(COV); 0163 % 0164 % cost=0; 0165 % for k=1:n 0166 % for j=1:n 0167 % %cost = cost + (cor_VkVj(k,j)-cor_COV(k,j))^2; 0168 % cost = cost + (clip_inv_sigmoid(cor_VkVj(k,j))-clip_inv_sigmoid(cor_COV(k,j)))^2; 0169 % end 0170 % end 0171 % %cost = cost / sum(sum(COV.^2)); 0172 % end 0173 % 0174 % p=fminunc(@costfun3,p,optimoptions('fminunc','MaxFunctionEvaluations',1e9,'OptimalityTolerance',1e-9)); 0175 0176 0177 [~,~,fitted_COV,coef] = varavgref_model(p); 0178 %bounding coefficients 0179 for i=1:n 0180 if coef(i) > 0.999 0181 coef(i) = 0.999; 0182 end 0183 end 0184 0185 % figure 0186 % pcolor(fitted_COV) 0187 % figure 0188 % pcolor(COV) 0189 0190 end