hw_1.m

文件同时包含了国外人脑分析的详细数据和处理程序.其中包括寻找匹配矩阵,数据遮盖,密度分析等等.

% created by Xiaoxu Kang on Feb. 6th, 2008
% for HW_1 in System Bioengineering 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc;
record = zeros(33,100);       
pos_mark = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]; 
load onset

%first generate X without considerations of K
stimuli_count = 1;
for row = 1:33
    begin = onset(stimuli_count);
    if ((begin+19)<100)
        record(row, begin:(begin+19)) = pos_mark;
    else 
        record(row, begin:100) = pos_mark(1:(length(pos_mark)-(begin+19-100)));
    end
    stimuli_count = stimuli_count+1;
end
record = record';
X_ini = zeros(100,20); 
for row = 1:100
    valid_part = find(record(row,:)>0);
    for i = 1:length(valid_part)
        X_ini(row, record(row,valid_part(i))) = 1;
    end
end

%generate K 
K = zeros(100,3);
K(:,1)= ones(100,1); 
for i = 0:99
    col_2(i+1) = cos((pi/99)*i);
end
K(:,2) = col_2';   
for i = 0:50
    col_3_fir(i+1) = cos((pi/50)*i);
end
i = 1;
for j = 48:(-1):0
    col_3_sec(i) = cos((pi/49)*j);
    i = i+1;
end
col_3 = [col_3_fir, col_3_sec];
K(:,3) = col_3';   

%generate X, which is Q1_X here.
Q1_X = [X_ini,K];
Q1_K = K;
figure(1);imagesc(Q1_X);colormap(jet);colorbar;title('Matrix X','fontsize',14);
figure(2);imagesc(Q1_K);colormap(jet);colorbar;title('Matrix K','fontsize',14);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load data;

%generate all b_k into one matrix B_K
B_K = zeros(3,100);
for col = 1:5017
    line = X(:,col);
    B_K(:,col) = inv(Q1_K'*Q1_K)*Q1_K'*line;
end

%generate all b into one matrix B
B = zeros(23,100);
for col = 1:5017
    line = X(:,col);
    B(:,col) = inv(Q1_X'*Q1_X)*Q1_X'*line;
end

%estimate Y0
Y0 = Q1_K*B_K;

%estimate Y1
Y1 = Q1_X*B;
figure(3);imagesc(Y1);title('The reconstructed data Y1','fontsize',14);colorbar;colormap(jet);

%calculate F and F Distribution
df1 = rank(full(Q1_X));
df0 = rank(full(Q1_K));
F = zeros(1,5017);
for col = 1:5017
    y = X(:,col);
    y0 = Y0(:,col);
    y1 = Y1(:,col);
    FS = (y1-y0)'*(y1-y0);
    RS = (y-y1)'*(y-y1);
    F(:,col) = FS*(100-df1)/(df1-df0)/RS;
end
figure(4);hist(F,50);colormap(jet);title('Distribution of F','fontsize',14);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load mask; 
count = 1;
[i,j] = find(mask==1);
for k = 1:length(j)
    mask(i(k),j(k)) = F(count);
    count = count+1;
end
figure(5);imagesc(mask);colormap(jet);colorbar;title('F match','fontsize',14);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load data;

B_Temp = B(1:20,:);
Target = B_Temp(:,find(F>4));
aver_HRF = mean(Target,2);
figure(6);plot(aver_HRF,'r','LineWidth',2,'Marker','*','MarkerSize',8);title('Mean HRF function');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load anatomical;
for i = 1:128
    for j = 1:128
        if(mask(i,j)>4)
            anatomical(i,j)=700;
        else 
            ;
        end
    end
end
figure(7);imagesc(anatomical);title('Highlightened higher activity','fontsize',14);colormap(jet);colorbar;