频谱分析代码片段
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,:); |
作者:木木
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步