0001 function [varargout]=rest_powerspectrum(AOperation, varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 if nargin<1, help rest_powerspectrum; return; end
0016
0017 persistent ReHo_WaveGraph_Cfg;
0018
0019 if ~mislocked(mfilename),mlock; end
0020
0021 switch upper(AOperation),
0022 case 'SHOWFLUCTUATION',
0023 if nargin~=5, error('Usage: rest_powerspectrum(''ShowFluctuation'', ADataDir, AVoxelPosition, ASamplePeriod, ABandRange);'); end
0024 ADataDir =varargin{1};
0025 AVoxelPosition =varargin{2};
0026 ASamplePeriod =varargin{3};
0027 ABandRange =varargin{4};
0028 theFig =ExistDisplayFigure(ReHo_WaveGraph_Cfg, ADataDir);
0029 isExistFig =rest_misc( 'ForceCheckExistFigure' , theFig);
0030 if ~isExistFig
0031
0032 ReHo_WaveGraph_Cfg.Config(1+GetDisplayCount(ReHo_WaveGraph_Cfg)) =InitControls(ADataDir, AVoxelPosition, ASamplePeriod, ABandRange);
0033
0034 theFig =ReHo_WaveGraph_Cfg.Config(GetDisplayCount(ReHo_WaveGraph_Cfg)).hFig;
0035 else
0036
0037 theCardinal =ExistDisplay(ReHo_WaveGraph_Cfg, ADataDir);
0038 if theCardinal>0,
0039 ReHo_WaveGraph_Cfg.Config(theCardinal).VoxelPosition =AVoxelPosition;
0040
0041 ReHo_WaveGraph_Cfg.Config(theCardinal).VoxelArray =[AVoxelPosition; ReHo_WaveGraph_Cfg.Config(theCardinal).VoxelArray];
0042
0043 ReHo_WaveGraph_Cfg.Config(theCardinal).SamplePeriod =ASamplePeriod;
0044 ReHo_WaveGraph_Cfg.Config(theCardinal).BandRange =ABandRange;
0045
0046
0047
0048
0049
0050 PlotFluctuation(ReHo_WaveGraph_Cfg.Config(theCardinal));
0051 end
0052 end
0053 figure(theFig);
0054 varargout{1} =theFig;
0055
0056 case 'ONFIGUREDELETE',
0057 if nargin~=2, error('Usage: rest_powerspectrum(''OnFigureDelete'', ADataDir);'); end
0058
0059 ADataDir =varargin{1};
0060 ReHo_WaveGraph_Cfg =DeleteFigure(ReHo_WaveGraph_Cfg, ADataDir);
0061
0062 case 'ONFIGURERESIZE',
0063 if nargin~=2, error('Usage: rest_powerspectrum(''ResizeFigure'', ADataDir);'); end
0064 ADataDir =varargin{1};
0065 theFig =ExistDisplayFigure(ReHo_WaveGraph_Cfg, ADataDir);
0066 isExistFig =rest_misc( 'ForceCheckExistFigure' , theFig);
0067 if isExistFig
0068 theCardinal =ExistDisplay(ReHo_WaveGraph_Cfg, ADataDir);
0069 if theCardinal>0
0070 OnFigureResize(ReHo_WaveGraph_Cfg.Config(theCardinal));
0071 end
0072 end
0073
0074 case 'ONOPTIONDETRENDBEFOREFFT',
0075 if nargin~=2, error('Usage: rest_powerspectrum(''OnOptionDetrendBeforeFFT'', ADataDir);'); end
0076 ADataDir =varargin{1};
0077 theFig =ExistDisplayFigure(ReHo_WaveGraph_Cfg, ADataDir);
0078 isExistFig =rest_misc( 'ForceCheckExistFigure' , theFig);
0079 if isExistFig
0080 theCardinal =ExistDisplay(ReHo_WaveGraph_Cfg, ADataDir);
0081 if theCardinal>0,
0082 if get(ReHo_WaveGraph_Cfg.Config(theCardinal).hDetrendBeforeFFT, 'Value'),
0083 ReHo_WaveGraph_Cfg.Config(theCardinal).DetrendBeforeFFT ='Yes';
0084 else
0085 ReHo_WaveGraph_Cfg.Config(theCardinal).DetrendBeforeFFT ='No';
0086 end
0087 PlotFluctuation(ReHo_WaveGraph_Cfg.Config(theCardinal));
0088 end
0089 end
0090 varargout{1} =theFig;
0091
0092 case 'QUITALLPOWERSPECTRUM',
0093 if nargin~=1, error('Usage: rest_powerspectrum(''QuitAllPowerSpectrum'');'); end
0094 for x=1:GetDisplayCount(ReHo_WaveGraph_Cfg),
0095 rest_powerspectrum('OnFigureDelete', ReHo_WaveGraph_Cfg.Config(x).DataDir);
0096 end
0097 clear ReHo_WaveGraph_Cfg;
0098
0099
0100 otherwise
0101 end
0102
0103 function Result =InitControls(ADataDir, AVoxelPosition, ASamplePeriod, ABandRange)
0104
0105 [theDataset, theVoxelSize, theImgFileList, Header]=rest_to4d(ADataDir);
0106 [nDim1, nDim2, nDim3, nDim4] =size(theDataset);
0107 rest_waitbar;
0108
0109
0110 theFig =figure('Units', 'pixel', 'MenuBar', 'none', 'Toolbar', 'figure', ...
0111 'NumberTitle', 'off', 'Name', ADataDir, ...
0112 'DeleteFcn', sprintf('rest_powerspectrum(''OnFigureDelete'', ''%s'');', rest_misc('ReplaceSingleQuota', ADataDir)) , ...
0113 'ResizeFcn', sprintf('rest_powerspectrum(''OnFigureResize'', ''%s'');', rest_misc('ReplaceSingleQuota', ADataDir)) );
0114 MarginX =10; MarginY =10;
0115 OffsetX =3*MarginX; OffsetY =MarginY;
0116
0117
0118 hVoxelPosition =uicontrol(theFig, 'Style','radiobutton', 'Units','pixels', ...
0119 'String', sprintf('(%d,%d,%d)',AVoxelPosition), ...
0120 'BackgroundColor', get(theFig,'Color'), ...
0121 'HorizontalAlignment', 'left', ...
0122 'Position',[OffsetX, OffsetY, 150,15]);
0123
0124
0125 OffsetX =3*MarginX; OffsetY =MarginY +25;
0126 hDetrendBeforeFFT =uicontrol(theFig, 'Style','checkbox', 'Units','pixels', ...
0127 'String', 'Remove Linear Trend Before Power Spectrum', ...
0128 'BackgroundColor', get(theFig,'Color'), ...
0129 'HorizontalAlignment', 'left', ...
0130 'Value', 1, ...
0131 'Callback', ...
0132 sprintf('rest_powerspectrum(''OnOptionDetrendBeforeFFT'', ''%s'');', rest_misc('ReplaceSingleQuota', ADataDir)), ...
0133 'Position',[OffsetX, OffsetY, 250,15]);
0134
0135
0136 OffsetX =6*MarginX; OffsetY =3*MarginY +15 +3*MarginY;
0137 hAxesTimeCourse =axes('Parent', theFig, ...
0138 'Units', 'pixel', 'DrawMode','fast', ...
0139 'Position', [OffsetX OffsetY 3*nDim4 150]);
0140
0141 OffsetX =6*MarginX; OffsetY =3*MarginY +15 +6*MarginY +150 +6*MarginY;
0142 hAxesAmplitude =axes('Parent', theFig, ...
0143 'Units', 'pixel', 'DrawMode','fast', ...
0144 'Position', [OffsetX OffsetY 3*nDim4 150], ...
0145 'NextPlot', 'replacechildren');
0146
0147 AConfig.hFig =theFig;
0148
0149 AConfig.hAxesTimeCourse =hAxesTimeCourse;
0150 AConfig.hAxesAmplitude =hAxesAmplitude;
0151
0152
0153 AConfig.hVoxelPosition =hVoxelPosition;
0154
0155 AConfig.hDetrendBeforeFFT =hDetrendBeforeFFT;
0156
0157
0158 AConfig.DataDir =ADataDir;
0159 AConfig.VoxelPosition =AVoxelPosition;
0160 AConfig.VoxelArray =AVoxelPosition;
0161
0162 AConfig.Dataset =theDataset;
0163 AConfig.SamplePeriod =ASamplePeriod;
0164 AConfig.BandRange =ABandRange;
0165 AConfig.DetrendBeforeFFT='Yes';
0166 Result =AConfig;
0167
0168
0169
0170 PlotFluctuation(AConfig);
0171
0172
0173 FigWidth =6*MarginX + 3*nDim4 + 6*MarginX;
0174 FigHeight =3*MarginY +15 +6*MarginY +150 +6*MarginY + 150+ 3*MarginY ;
0175 thePos =get(theFig, 'Position');
0176 theScreenSize =get(0,'ScreenSize');
0177 if thePos(1)>= theScreenSize(3)
0178 thePos(1) =theScreenSize(1);
0179 end
0180 if (thePos(2) +FigHeight) +100>= theScreenSize(4)
0181 thePos(2) =theScreenSize(4) -FigHeight -100;
0182 end
0183 thePos =[thePos(1), thePos(2), FigWidth,FigHeight];
0184 set(theFig, 'Position', thePos);
0185
0186 OnFigureResize(AConfig);
0187
0188 return;
0189
0190
0191
0192 function Result =DeleteFigure(AGlobalConfig, ADataDir)
0193 x =ExistDisplay(AGlobalConfig, ADataDir);
0194 if x>0
0195 theDisplayCount =GetDisplayCount(AGlobalConfig);
0196 isExistFig =rest_misc( 'ForceCheckExistFigure' , AGlobalConfig.Config(x).hFig);
0197 if isExistFig
0198 delete(AGlobalConfig.Config(x).hFig);
0199 if theDisplayCount>x
0200 for y=x:theDisplayCount-1
0201 AGlobalConfig.Config(x) =AGlobalConfig.Config(x+1);
0202 end
0203 end
0204 AGlobalConfig.Config(theDisplayCount)=[];
0205 end
0206 end
0207 Result =AGlobalConfig;
0208 function Result =GetDisplayCount(AGlobalConfig)
0209
0210 if isempty(AGlobalConfig) || isempty(AGlobalConfig.Config),
0211 Result =0;
0212 else
0213 Result =length(AGlobalConfig.Config);
0214 end
0215 return;
0216
0217 function Result =ExistDisplayFigure(AGlobalConfig, ADataDir)
0218
0219 Result =-1;
0220 theCardinal =ExistDisplay(AGlobalConfig, ADataDir);
0221 if theCardinal>0,
0222 Result =AGlobalConfig.Config(theCardinal).hFig;
0223 end
0224
0225 function Result =ExistDisplay(AGlobalConfig, ADataDir)
0226
0227 Result =0;
0228 if (isstruct(AGlobalConfig) && isstruct(AGlobalConfig.Config))
0229 for x=1:length(AGlobalConfig.Config)
0230 if strcmpi( AGlobalConfig.Config(x).DataDir, ADataDir)
0231 Result =x;
0232 return;
0233 end
0234 end
0235 else
0236 return;
0237 end
0238
0239
0240
0241 function PlotFluctuation(AConfig);
0242 theX =AConfig.VoxelPosition(1);
0243 theY =AConfig.VoxelPosition(2);
0244 theZ =AConfig.VoxelPosition(3);
0245
0246
0247 if all([AConfig.VoxelPosition,0]<=size(AConfig.Dataset)) && all(AConfig.VoxelPosition>=[1 1 1]),
0248 else
0249 error(sprintf('Illegal voxel position: (%s)', num2str(AConfig.VoxelPosition)));
0250 end
0251
0252
0253 theTimeCourse =squeeze(AConfig.Dataset(theX, theY, theZ, :));
0254 axes(AConfig.hAxesTimeCourse); cla;
0255 if rest_misc('GetMatlabVersion')>=7.3,
0256 plot([1:length(theTimeCourse)] *AConfig.SamplePeriod, ...
0257 theTimeCourse,'blue', 'DisplayName', 'Time Course');
0258 else
0259 plot([1:length(theTimeCourse)] *AConfig.SamplePeriod, ...
0260 theTimeCourse,'blue');
0261 end
0262 xlim([1, length(theTimeCourse)] *AConfig.SamplePeriod);
0263 theYLim =[min(theTimeCourse), max(theTimeCourse)];
0264 if ~isfinite(theYLim(1)), theYLim(1)=0; end
0265 if ~isfinite(theYLim(2)), theYLim(2)=0; end
0266 if theYLim(2)>theYLim(1), ylim(theYLim); end
0267 set(gca, 'Title',text('String','Time Course', 'Color', 'magenta'));
0268 xlabel('Time( seconds)');
0269 ylabel('Intensity');
0270
0271
0272 thePaddedLen =rest_nextpow2_one35(length(theTimeCourse));
0273 if strcmpi(AConfig.DetrendBeforeFFT, 'Yes'),
0274
0275 theTimeCourseNoTrend =detrend(double(theTimeCourse)) ;
0276
0277
0278 axes(AConfig.hAxesTimeCourse); hold on;
0279 thePlotTimeCourse =theTimeCourseNoTrend + repmat(mean(double(theTimeCourse)), [length(theTimeCourse), 1]);
0280 if rest_misc('GetMatlabVersion')>=7.3,
0281 plot([1:length(theTimeCourse)] *AConfig.SamplePeriod, ...
0282 thePlotTimeCourse, 'g:', 'DisplayName', 'Detrended Time Course');
0283 else
0284 plot([1:length(theTimeCourse)] *AConfig.SamplePeriod, ...
0285 thePlotTimeCourse, 'g:');
0286 end
0287 theYLim =[min(theYLim(1),min(thePlotTimeCourse)), max(theYLim(2),max(thePlotTimeCourse))];
0288 if ~isfinite(theYLim(1)), theYLim(1)=0; end
0289 if ~isfinite(theYLim(2)), theYLim(2)=0; end
0290 if theYLim(2)>theYLim(1), ylim(theYLim); end
0291 set(gca, 'Title',text('String','Time Course(Green dot line is after removing linear trend)', 'Color', 'magenta'));
0292
0293
0294 thePowerTitle ='Power Spectrum after removing linear trend';
0295
0296 theFreqSeries =fft(theTimeCourseNoTrend, thePaddedLen);
0297 else
0298
0299 thePowerTitle ='Power Spectrum';
0300
0301 theFreqSeries =fft(double(theTimeCourse), thePaddedLen);
0302 end
0303
0304 theSampleFreq =1/AConfig.SamplePeriod ;
0305 theFreqPrecision =theSampleFreq/thePaddedLen;
0306 theFreqLim =[theFreqPrecision: theFreqPrecision :theSampleFreq/2];
0307 theXLim =[2,(thePaddedLen/2 +1)];
0308
0309
0310 theFreqSeries =abs(theFreqSeries([theXLim(1):theXLim(2)]));
0311 theFreqSeries(1:end) =theFreqSeries(1:end).^2 /length(theTimeCourse);
0312
0313
0314
0315
0316
0317 theFreqSeries(1:end-1) =theFreqSeries(1:end-1) *2;
0318
0319
0320 axes(AConfig.hAxesAmplitude); cla;
0321 if rest_misc('GetMatlabVersion')>=7.3,
0322 plot(theFreqLim, theFreqSeries, ...
0323 'Color', 'blue', 'DisplayName', 'Power Spectrum');
0324 else
0325 plot(theFreqLim, theFreqSeries, 'Color', 'blue');
0326 end
0327
0328 xlim([theFreqLim(1) , theFreqLim(end)]);
0329 theYLim =[min(theFreqSeries(1:end)), max(theFreqSeries(1:end))];
0330 if ~isfinite(theYLim(1)), theYLim(1)=0; end
0331 if ~isfinite(theYLim(2)), theYLim(2)=0; end
0332 if theYLim(2)>theYLim(1), ylim(theYLim); end
0333 set(gca, 'Title',text('String',thePowerTitle, 'Color', 'magenta'));
0334 xlabel(sprintf('Frequency( Hz), Sample Period( TR)=%g seconds', AConfig.SamplePeriod));
0335 ylabel('Amplitude');
0336
0337
0338 if rest_misc('GetMatlabVersion')>=7.0,
0339 hDataCursor = datacursormode(AConfig.hFig);
0340 set(hDataCursor,'DisplayStyle','datatip', 'SnapToDataVertex','on','Enable','on')
0341 set(hDataCursor,'UpdateFcn', @SetDataTip);
0342 end
0343
0344
0345 hold on;
0346 if rest_misc('GetMatlabVersion')>=7.3,
0347 plot([1, 1]*AConfig.BandRange(1), get(gca,'Ylim'), 'r:', 'DisplayName', 'Band Limit');
0348 else
0349 plot([1, 1]*AConfig.BandRange(1), get(gca,'Ylim'), 'r:');
0350 end
0351 text(AConfig.BandRange(1),theYLim(2)-theYLim(2)/10, ...
0352 sprintf('\\leftarrow %g Hz', AConfig.BandRange(1)),...
0353 'HorizontalAlignment','left', 'Color', 'red');
0354 if rest_misc('GetMatlabVersion')>=7.3,
0355 plot([1, 1]*AConfig.BandRange(2), get(gca,'Ylim'), 'r:', 'DisplayName', 'Band Limit');
0356 else
0357 plot([1, 1]*AConfig.BandRange(2), get(gca,'Ylim'), 'r:');
0358 end
0359 text(AConfig.BandRange(2),theYLim(2)-theYLim(2)/10, ...
0360 sprintf('\\leftarrow %g Hz', AConfig.BandRange(2)),...
0361 'HorizontalAlignment','left', 'Color', 'red');
0362
0363
0364 theName =sprintf('(%d,%d,%d)',AConfig.VoxelPosition);
0365
0366
0367 set(AConfig.hVoxelPosition, ...
0368 'String', theName, ...
0369 'TooltipString', sprintf('Time Course: mean=%g, std=%g\n', ...
0370 mean(double(theTimeCourse)), std(double(theTimeCourse))) );
0371
0372 OnFigureResize(AConfig);
0373
0374 function OnFigureResize(AConfig)
0375 MarginX =10; MarginY =10;
0376
0377 theFigPos =get(AConfig.hFig, 'Position');
0378 FigWidth =theFigPos(3);
0379 FigHeight=theFigPos(4);
0380
0381
0382
0383 theAxesTimeCoursePos =get(AConfig.hAxesTimeCourse, 'Position');
0384 theAxesTimeCoursePos(3) =FigWidth -6*MarginX -3*MarginX;
0385
0386 theAxesTimeCoursePos(2) =6*MarginY +15;
0387 theAxesTimeCoursePos(4) =(FigHeight -3*MarginY)/2 -9*MarginY;
0388 set(AConfig.hAxesTimeCourse, 'Position', theAxesTimeCoursePos);
0389
0390
0391 theAxesAmplitudePos =get(AConfig.hAxesAmplitude, 'Position');
0392 theAxesAmplitudePos(3) =FigWidth -6*MarginX -3*MarginX;
0393
0394 theAxesAmplitudePos(2) =theAxesTimeCoursePos(2) +theAxesTimeCoursePos(4) +9*MarginY;
0395 theAxesAmplitudePos(4) =(FigHeight -3*MarginY)/2-9*MarginY;
0396 set(AConfig.hAxesAmplitude, 'Position', theAxesAmplitudePos);
0397
0398
0399 thePos =get(AConfig.hVoxelPosition,'Position');
0400 thePos(3) =FigWidth -thePos(1) -MarginX;
0401 set(AConfig.hVoxelPosition,'Position', thePos);
0402
0403 drawnow;
0404
0405 function Result =SetDataTip(empt, event_obj)
0406 pos = get(event_obj,'Position');
0407 hTarget = get(event_obj,'Target');
0408 theName = get(hTarget,'DisplayName');
0409 theIdx =get(event_obj,'DataIndex');
0410 if strcmpi(theName(1:5), 'Power') ,
0411 Result = { [theName], ...
0412 sprintf('\nFrequency: \t%g Hz', pos(1)), ...
0413 sprintf('Power: \t%g',pos(2)), ...
0414 sprintf('\nIndex: \t%d', theIdx) };
0415 elseif any(findstr(theName, 'Time')),
0416 Result = { sprintf('%s\n',theName), ...
0417 sprintf('Time: %g seconds',pos(1)), ...
0418 ['Intensity: ',num2str(pos(2))] , ...
0419 sprintf('\nIndex: \t%d', theIdx) };
0420 else
0421 Result = { sprintf('%s\n',theName), ...
0422 ['X: ',num2str(pos(1))], ...
0423 ['Y: ',num2str(pos(2))] };
0424 end