%% 原始数据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;