频谱分析代码片段

ldcca_tms = img_To_4D_array('C:\Users\Administrator\Desktop\contrast\2014-05-20-12-16.img');
spm_tms = img_To_4D_array('C:\Users\Administrator\Desktop\contrast\no_phycaa.img');

% x = size(spm_tms,1);
% y = size(spm_tms,2);
% z = size(spm_tms,3);
%
% %spm_tms_1d = reshape(spm_tms ,[1,x*y*z] );
%
%
% x = size(ldcca_tms,1);
% y = size(ldcca_tms,2);
% z = size(ldcca_tms,3);
%
% % ldcca_tms_1d = reshape( ldcca_tms , [1,x*y*z] );

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;

%% 先提取出被去噪去除的点的时间过程  实验结果，去除了435个点，

% left_spm_tms = spm_tms(inv_mask_ldcca_tms);

tmp_spm = spm_tms .* inv_mask_ldcca_tms;

%>>这是最后的模板
mask_big_left_spm_tms = tmp_spm>0;

% big_left_spm_tms = left_spm_tms(left_spm_tms>0);

% aa = size(big_left_spm_tms)

%% 利用我们的去噪算法,得到的新的点，一个点都没有增加

% left_ldcca_tms = ldcca_tms(ldcca_tms);

tmp_ldcca = ldcca_tms .* inv_mask_spm_tms;

%>>这是最后的模板
mask_big_left_ldcca_tms = tmp_ldcca > 0;

% big_left_ldcca_tms = left_ldcca_tms( left_ldcca_tms>0 );

% bb = size(big_left_ldcca_tms)

%% 提取被去除的点的时间过程
%--brain functional 4d data

datacell_4d   = load_untouch_nii('C:\Users\Administrator\Desktop\workspace\phycaa_plus_2104_03_27\func_3d.nii');

%datacell_4d   = load_untouch_nii('F:\4Dfile\func_3d_PHYCAA_step1+2.nii');

ldim = size(datacell_4d.img);

fullMsk = repmat( mask_big_left_spm_tms, [1,1,1,ldim(4)] );
left_spm_tms_2d = reshape( datacell_4d.img(fullMsk>0), [], ldim(4) );

inv_fullMsk =  repmat( mask_spm_tms.*(~mask_big_left_spm_tms), [1,1,1,ldim(4)] );
inv_left_spm_tms_2d = reshape( datacell_4d.img(inv_fullMsk>0), [], ldim(4) );

%% fft decomposition

TR =2;

% parameters for spectral estimation
Fny    = 0.5 * (1/ TR);                % nyquist frequency
NFFT   = 2^nextpow2(70);           % Next power of 2 from length of time-axis
f      = Fny*linspace(0,1,NFFT/2+1);  % fourier data corresponds to these frequency points
numlow = sum( f <= 0.1 );           % count the number of frequency points below threshold

% =========================================================================
%   get ordered spectral estimates

powHgh =[];
powLow =[];

inv_powHgh=[];
inv_powLow=[];

i=1;

%for i=1:size(left_spm_tms_2d,1)
% get sum of spectral power values, for f > dataInfo.FreqCut
powMat  = abs( fft( double(left_spm_tms_2d) , NFFT ,2) ) / 70;
powMat  = powMat(:,1:NFFT/2+1);
powHgh  = sum( sqrt(powMat(:,numlow+1:end)), 2 );
powLow  = sum( sqrt(powMat(:,1:numlow)), 2 );
%end

mean_powHgh = mean(powHgh);
median_powHgh = median(powHgh);
min_powHgh = min(powHgh);

mean_powLow = mean(powLow);
median_powLow = median(powLow);
min_powLow = min(powLow);

%for i=1:size(inv_left_spm_tms_2d,1)
% get sum of spectral power values, for f > dataInfo.FreqCut
inv_powMat  = abs( fft( double(inv_left_spm_tms_2d) , NFFT ,2) ) / 70;
inv_powMat  = inv_powMat(:,1:NFFT/2+1);
inv_powHgh  = sum( sqrt(inv_powMat(:,numlow+1:end)), 2 );
inv_powLow  = sum( sqrt(inv_powMat(:,1:numlow)), 2 );
%end

inv_mean_powHgh = mean(inv_powHgh);
inv_median_powHgh = median(inv_powHgh);
inv_min_powHgh = min(inv_powHgh);

inv_mean_powLow = mean(inv_powLow);
inv_median_powLow = median(inv_powLow);
inv_min_powLow = min(inv_powLow);

%vol_4d(:,:,slice,:);