rfchop_start

 rfchop_start.m

0001 % rfchop_start.m
0002 
0003 global dataroot proc_subs
0004 
0005 SCAN_OFFSET = 0;
0006 
0007 CONVERT_FORMAT = []; %
0008 GLOBAL_SPECT = 1;
0009 DO_CALC = 1;
0010 PLOT_SPECT = 1;
0011 LOWPASS = 1;
0012 
0013 mask_threshold = 500;
0014 Fs = 1/2;
0015 
0016 proc_subs = 1:2;
0017 
0018 s = get_rfchop_sinfo;
0019 
0020 %
0021 % Convert data from GE format into AVW
0022 %
0023 
0024 if CONVERT_FORMAT
0025   convert_format(s(CONVERT_FORMAT),dataroot,dataroot, SCAN_OFFSET);
0026 end
0027 
0028 if GLOBAL_SPECT
0029   if DO_CALC
0030     clear avgmag stdmag ss
0031   end
0032   for isub = 1:length(proc_subs)
0033     subidx = proc_subs(isub);
0034 
0035     fprintf('Subject: %s\n', s(subidx).id);
0036     
0037     subdir = fullfile(dataroot, s(subidx).id);
0038     nruns = length(s(subidx).use_runs);
0039     
0040     if DO_CALC
0041       for irun = 1:nruns
0042     runidx = s(subidx).use_runs(irun);
0043     rundir = fullfile(subdir,'epi',sprintf('run%d',runidx));
0044     
0045     fprintf('Run (%d/%d): %d\n', irun, nruns, runidx)
0046     
0047     % Load all of the images
0048     flist = dir(fullfile(rundir,sprintf('%s*i0*.img',s(subidx).id)));
0049     V = spm_vol(fullfile(rundir,flist(1).name));
0050     if length(flist) ~= s(subidx).nvol(runidx)
0051       error('Wrong number of files in directory')
0052     end
0053     
0054     fprintf('Loading data\n')
0055     data = zeros(s(subidx).nvol(runidx),prod(V.dim(1:3)));
0056     nvol = s(subidx).nvol(runidx);
0057     for ivol = 1:nvol
0058       fprintf('\t%s\n', flist(ivol).name);
0059       fid = fopen(fullfile(rundir,flist(ivol).name),'rb');
0060       data(ivol,:) = fread(fid,spm_type(V.dim(4)))';
0061       fclose(fid);
0062     end % for ivol=
0063     
0064     % Get voxels that exceed threshold
0065     vox = find(data(1,:) > mask_threshold);
0066 
0067     % Remove mean from each voxel
0068     data = detrend(data,'constant');
0069     
0070     % Filter if desired
0071     if LOWPASS
0072       [b,a] = butter(5,0.8);
0073       for ivox = 1:length(vox)
0074         data(:,vox(ivox)) = filtfilt(b,a,data(:,vox(ivox)));
0075       end
0076     end
0077     
0078     % Compute FFTs
0079     fprintf('Computing spectra ...\n');
0080     magdata = abs(fft(data(:,vox)))/nvol;
0081     avgmag{isub}(:,irun) = mean(magdata,2);
0082     stdmag{isub}(:,irun) = std(magdata,[],2);
0083 
0084     
0085     % Get SS across all voxels
0086     ss{isub}(irun,:) = sum(data.^2);
0087     avg_ss{isub}(irun) = mean(ss{isub}(irun,:),2);
0088     std_ss{isub}(irun) = std(ss{isub}(irun,:),[],2);
0089       end % for irun=
0090     end % if DO_CALC
0091     
0092     if PLOT_SPECT
0093       figure(isub),clf
0094 
0095       nvol = s(subidx).nvol(runidx);
0096 
0097       fscale = (0:fix(nvol/2)-1)/fix(nvol/2)*(Fs/2);
0098       plot(fscale,(avgmag{isub}(1:fix(nvol/2),:).^2)./repmat(avg_ss{isub},fix(nvol/2),1))
0099       xlabel('Frequency (Hz)')
0100       ylabel('Power')
0101       legend(strrep(cond_names(s(subidx).cond_order(s(subidx).use_runs)),'_','\_'))
0102     end
0103   end % for isub=
0104 end % if GLOBAL_SPECT
0105