频谱分析代码片段

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
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,:);

  

posted @   二郎那个三郎  阅读(502)  评论(0编辑  收藏  举报
(评论功能已被禁用)
点击右上角即可分享
微信分享提示