mne_source_spectral_analysis

 [res] = mne_source_spectral_analysis(fname_rawdata, fname_output, cfg)

 Estimate frequency spectra in the source space and optionally write out
 stc files which have frequencies along the "time" axis.

 fname_data     - Name of the data file

 MNE inversion
 cfg.inv        - Inverse operator structure or file name
 cfg.lambda2    - The regularization factor
 cfg.dSPM       - enable dSPM: 0 or 1

 Spectral estimation
 cfg.mode       - output quantity; 'amplitude', 'power', 'phase'
 cfg.window     - window type: 'hanning', 'hamming' or any other window
                  function available on Matlab
 cfg.fft_length - FFT length in samples (half-overlapping windows used)
 cfg.foi        - Frequencies of interest as [f_min f_max]

 Output
 cfg.outfile    - The stem of the output STC file name holding the spectra

 (C)opyright Lauri Parkkonen, 2012

 $Log$

0001 function [res] = mne_source_spectral_analysis(fname_rawdata, cfg)
0002 %
0003 % [res] = mne_source_spectral_analysis(fname_rawdata, fname_output, cfg)
0004 %
0005 % Estimate frequency spectra in the source space and optionally write out
0006 % stc files which have frequencies along the "time" axis.
0007 %
0008 % fname_data     - Name of the data file
0009 %
0010 % MNE inversion
0011 % cfg.inv        - Inverse operator structure or file name
0012 % cfg.lambda2    - The regularization factor
0013 % cfg.dSPM       - enable dSPM: 0 or 1
0014 %
0015 % Spectral estimation
0016 % cfg.mode       - output quantity; 'amplitude', 'power', 'phase'
0017 % cfg.window     - window type: 'hanning', 'hamming' or any other window
0018 %                  function available on Matlab
0019 % cfg.fft_length - FFT length in samples (half-overlapping windows used)
0020 % cfg.foi        - Frequencies of interest as [f_min f_max]
0021 %
0022 % Output
0023 % cfg.outfile    - The stem of the output STC file name holding the spectra
0024 %
0025 % (C)opyright Lauri Parkkonen, 2012
0026 %
0027 % $Log$
0028 %
0029 
0030 me='MNE:mne_spectral_analysis';
0031 
0032 global FIFF;
0033 if isempty(FIFF)
0034    FIFF = fiff_define_constants();
0035 end
0036 
0037 if nargin ~= 2
0038    error(me,'Incorrect number of arguments'); 
0039 end
0040 
0041 % Set defaults
0042 if ~isfield(cfg, 'lambda2')
0043     cfg.lambda2 = 1.0;
0044 end
0045 if ~isfield(cfg, 'dSPM')
0046     cfg.dSPM = 0;
0047 end
0048 if ~isfield(cfg, 'mode')
0049    cfg.mode = 'amplitude';
0050 end
0051 if ~isfield(cfg, 'window')
0052     cfg.window = 'hanning';
0053 end
0054 if ~isfield(cfg, 'fft_length')
0055     cfg.fft_length = 1024;
0056 end
0057 
0058 % Set up raw data input
0059 rawfile = fiff_setup_read_raw(fname_rawdata);
0060 if ~isfield(cfg, 'first_samp') || cfg.first_samp <= 0
0061     cfg.first_samp = rawfile.first_samp;
0062 end
0063 if ~isfield(cfg, 'last_samp') || cfg.last_samp == Inf
0064     cfg.last_samp = rawfile.last_samp;
0065 end
0066 
0067 % Make a selector for frequencies
0068 faxis = ((1:cfg.fft_length/2) - 1) / cfg.fft_length * rawfile.info.sfreq;
0069 if ~isfield(cfg, 'foi')
0070     cfg.foi = [faxis(1) faxis(end)];
0071 end
0072 foi_i = find(faxis >= cfg.foi(1) & faxis <= cfg.foi(2));
0073 fprintf(1, 'Considering frequencies %g ... %g Hz\n', faxis(foi_i(1)), faxis(foi_i(end)));
0074 
0075 % Compose the inverse operator
0076 if isstruct(cfg.inv)
0077     inv = cfg.inv;
0078 else
0079     inv = mne_read_inverse_operator(cfg.inv);
0080 end
0081 inv = mne_prepare_inverse_operator(inv, 1, cfg.lambda2, cfg.dSPM);
0082 chsel = fiff_pick_channels(rawfile.info.ch_names, inv.noise_cov.names);
0083 fprintf(1,'Using %d channels for the inverse solution\n', length(chsel));
0084 fprintf(1,'Computing inverse...');
0085 inverse_noleads = diag(sparse(inv.reginv))*inv.eigen_fields.data*inv.whitener*inv.proj;
0086 if inv.eigen_leads_weighted
0087    %
0088    %     R^0.5 has been already factored in
0089    %
0090    fprintf(1,'(eigenleads already weighted)...');
0091    inverse = inv.eigen_leads.data*inverse_noleads;
0092 else
0093    %
0094    %     R^0.5 has to factored in
0095    %
0096    fprintf(1,'(eigenleads need to be weighted)...');
0097    inverse = diag(sparse(sqrt(inv.source_cov.data)))*inv.eigen_leads.data*inverse_noleads;
0098 end
0099 fprintf(1, '[done]\n');
0100 
0101 % Compute the average spectra
0102 ave = [];
0103 nspectra = 0;
0104 windowfun =feval(cfg.window, cfg.fft_length);
0105 for w = cfg.first_samp : cfg.fft_length/2 : cfg.last_samp
0106     data_sensor = fiff_read_raw_segment(rawfile, w, w + cfg.fft_length - 1, chsel);
0107     % Check for an incomplete last segment
0108     if size(data_sensor, 2) ~= cfg.fft_length
0109         fprintf(1, 'Skipping the remaining incomplete window\n');
0110         break;
0111     end
0112     nspectra = nspectra + 1;
0113     for ch = 1:size(data_sensor, 1)
0114         data_sensor(ch, :) = data_sensor(ch, :) .* windowfun';
0115     end
0116     data_fft = fft(data_sensor, cfg.fft_length, 2);
0117     data_fft = data_fft(:, foi_i);
0118     data_fft_source = inverse * double(data_fft);
0119     if isempty(ave)
0120         ave = zeros(size(data_fft_source));
0121     end
0122     ave = ave + abs(data_fft_source);
0123     fprintf('Window %d, %3.3g%% completed\n', nspectra, double(w - cfg.first_samp + cfg.fft_length/2) / double(cfg.last_samp - cfg.first_samp) * 100);
0124 end
0125 ave = ave / nspectra;
0126 
0127 % Combine the three source orientations if needed
0128 if inv.source_ori == FIFF.FIFFV_MNE_FREE_ORI
0129     % Do not do this for phase maps; just pick the normal component
0130     fprintf(1,'Combining the spectra of the current components\n');
0131     ave1 = zeros(size(ave,1)/3, size(ave,2));
0132     for k = 1:size(ave,2)
0133         ave1(:,k) = sqrt(mne_combine_xyz(ave(:,k)));
0134     end
0135     ave = ave1;
0136 end
0137 
0138 % Do dSPM weighting if requested
0139 if cfg.dSPM
0140     fprintf(1,'Applying dSPM weighting\n');
0141     ave = inv.noisenorm * ave;
0142 end
0143 
0144 % Compose the result structure
0145 res.inv     = inv;
0146 res.rawfile = rawfile;
0147 res.spectra = ave;
0148 res.fmin    = faxis(foi_i(1));
0149 res.fmax    = faxis(foi_i(end));
0150 res.fstep   = faxis(2) - faxis(1);
0151 res.faxis   = faxis;
0152 
0153 % Write out the STC file
0154 if isfield(cfg, 'outfile')
0155     fprintf('Writing out STC files...\n');
0156     mne_write_inverse_sol_stc(cfg.outfile, res.inv, res.spectra, 1e-3 * res.fmin, 1e-3 * res.fstep);
0157 end
0158 
0159 fprintf(1,'[done]\n');
0160 
0161 return;
0162 end