频谱分析代码片段2

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