频谱分析代码片段

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,:);