频谱分析代码片段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;
作者:木木