%% 原始数据
datacell_4d = load_untouch_nii('C:\Users\Administrator\Desktop\workspace\phycaa_plus_2104_03_27\func_4d.nii');
ldim = size(datacell_4d.img);
%% 制造模板
ldcca_tms = img_To_4D_array('C:\Users\Administrator\Desktop\contrast\2014-05-26-20-00.img');
spm_tms = img_To_4D_array('C:\Users\Administrator\Desktop\contrast\no_phycaa.img');
mask_ldcca_tms = ldcca_tms > 0;
inv_mask_ldcca_tms = ~mask_ldcca_tms;
mask_spm_tms = spm_tms > 0;
inv_mask_spm_tms = ~mask_spm_tms;
tmp_spm = spm_tms .* inv_mask_ldcca_tms;
mask_big_left_spm_tms = tmp_spm>0;
tmp_ldcca = ldcca_tms .* inv_mask_spm_tms;
mask_big_left_ldcca_tms = tmp_ldcca > 0;
%% 原始数据激活区能量谱均值
original_fullMsk = repmat( mask_spm_tms, [1,1,1,ldim(4)] );
original2d_img = reshape( datacell_4d.img(original_fullMsk>0), [], ldim(4) );
TR =2;
Fny = 0.5 * (1/ TR);
NFFT = 2^nextpow2(70);
f = Fny*linspace(0,1,NFFT/2+1);
powMat = abs( fft( double(original2d_img) , NFFT ,2) ) / 70;
powMat = powMat(:,1:NFFT/2+1);
figure(1);
mean_powMat = mean(powMat);
plot(f,mean_powMat);
hold on;
%% 新数据激活区能量谱均值
new_fullMsk = repmat( mask_ldcca_tms, [1,1,1,ldim(4)] );
new2d_img = reshape( datacell_4d.img(new_fullMsk>0), [], ldim(4) );
TR =2;
Fny = 0.5 * (1/ TR);
NFFT = 2^nextpow2(70);
f = Fny*linspace(0,1,NFFT/2+1);
powMat = abs( fft( double(new2d_img) , NFFT ,2) ) / 70;
powMat = powMat(:,1:NFFT/2+1);
mean_powMat = mean(powMat);
plot(f,mean_powMat,'Color','red');
%% 去除激活区的能量谱均值
substruct_fullMsk = repmat( mask_big_left_spm_tms, [1,1,1,ldim(4)] );
substruct_img_2d = reshape( datacell_4d.img(substruct_fullMsk>0), [], ldim(4) );
TR =2;
Fny = 0.5 * (1/ TR);
NFFT = 2^nextpow2(70);
f = Fny*linspace(0,1,NFFT/2+1);
powMat = abs( fft( double(substruct_img_2d) , NFFT ,2) ) / 70;
powMat = powMat(:,1:NFFT/2+1);
mean_powMat = mean(powMat);
plot(f,mean_powMat,'Color','blue');
hold on;
%% 新增加的激活区的能量谱均值
add_fullMsk = repmat( mask_big_left_ldcca_tms, [1,1,1,ldim(4)] );
add_spm_tms_2d = reshape( datacell_4d.img(add_fullMsk>0), [], ldim(4) );
TR =2;
Fny = 0.5 * (1/ TR);
NFFT = 2^nextpow2(70);
f = Fny*linspace(0,1,NFFT/2+1);
powMat = abs( fft( double(add_spm_tms_2d) , NFFT ,2) ) / 70;
powMat = powMat(:,1:NFFT/2+1);
mean_powMat = mean(powMat);
plot(f,mean_powMat,'Color','yellow');
hold on;