Reconstruct ICA components in Matlab command window

alirezatalesh · November 7, 2017, 8:58am

Hi,
I computed my data ICA components using brainstorm. Now, I want to have time series of components in Matlab.

How can I do that? I think I should save active projectors by using option in “Select active projectors”, then apply pinv(Projector.Components)*data. Is it a correct manner?

Best,
Talesh

Francois · November 7, 2017, 9:53am

You build the projectors with function process_ssp2(‘BuildProjector’, ListProj, 1):

github.com

brainstorm-tools/brainstorm3/blob/master/toolbox/process/functions/process_ssp2.m#L867


nProj    = size(U,2);
proj.Components = zeros(nChannel, nProj);
proj.Components(chanmask == 1,:) = U;
% Other fields
proj.SingVal  = diag(S)';
proj.CompMask = zeros(1,nProj);
proj.Status   = 1;
end




%% ===== BUILD PROJECTOR ====
% Combine all the projectors in decomposed form to create a [nChan x nChan] matrix
%
% USAGE: Projector = process_ssp2('BuildProjector', ListProj, ProjStatus)
%
% INPUT:
%    - ListProj   : Array of db_template('projector')
%    - ProjStatus : List of the projector status to include 
%                   0 = not applied, not used
%                   1 = have to be used, but still have to be applied on the fly
%                   2 = already applied, saved in the file, not revertible

ListProj: Get it from the channel file (field .Projector).
This creates a Projector matrix, that you can then simply multiply with the data. Example:

github.com

brainstorm-tools/brainstorm3/blob/master/toolbox/io/in_fread.m#L234


%% ===== GRADIENT CORRECTION =====

% 3rd-order gradient correction

if ~isempty(ImportOptions) && ImportOptions.UseCtfComp && ~strcmpi(sFile.format, 'BST-DATA') && ~isempty(ChannelMat) && ~isempty(ChannelMat.MegRefCoef) && (sFile.prop.currCtfComp ~= sFile.prop.destCtfComp)

iMeg = good_channel(ChannelMat.Channel,[],'MEG');

iRef = good_channel(ChannelMat.Channel,[],'MEG REF');

if ~isempty(iChannels) && (length(iChannels) ~= length(ChannelMat.Channel))

    error('CTF compensators require that you read all the channels at the same time.');

else

    F(iMeg,:) = F(iMeg,:) - ChannelMat.MegRefCoef * F(iRef,:);

end

end



%% ===== SSP PROJECTORS =====

if ~isempty(ImportOptions) && ImportOptions.UseSsp && ~strcmpi(sFile.format, 'BST-DATA') && ~isempty(ChannelMat) && ~isempty(ChannelMat.Projector)

% Build projector matrix

Projector = process_ssp2('BuildProjector', ChannelMat.Projector, 1);

% Get bad channels

iBadChan = find(sFile.channelflag == -1);

% Apply projector

if ~isempty(Projector)

    % Remove bad channels from the projector (similar as in process_megreg)

Alternatively, you can display the IC time series from the window “Select active projectors”, then export the displayed time series (unselect the [DC] button in the Record tab): right-click on the figure > Snapshot > Export to Matlab.

alirezatalesh · November 7, 2017, 10:30am

Dear Francois,

I can plot the ICA component in brainstorm from the window “Select active projectors”, but I also want to see every component PSD. How can I do that?

Francois · November 7, 2017, 11:06am

You can try the following:

Open the recordings, from the “Select active projectors”, display the IC time series
In the Record tab, click on the drop-down menu with the montage list > Edit montage
Duplicate the montage “ICA components[tmp]”, rename it into something like “IC” (the important is to remove the [tmp] tag from the name)
Import the recordings in the database (right-click on the Link to raw file > Import in database), possibly as a long continuous segment
Select the imported file in Process1, run the process Standardize > Apply montage > “IC”
You should be able to compute the PSD from this imported file
If you need to review this as a continuous file, right-click on it > Review as raw

alirezatalesh · November 7, 2017, 1:42pm

Thanks so much dear Francois, your suggestions have always been useful for me. But, when I process the step 6, the error "No sensors are selected " occurred. My data is recorded with 256 channel and my components are 251.

The warm regards,
Talesh

Francois · November 8, 2017, 10:10am

I agree, there should be a button for this. It “just has to be done”. I added this to our infinite to-do list: "Save IC time series in database "
http://neuroimage.usc.edu/brainstorm/Next#Processes

About your error: you must configure your subjects with the option “Default channel : No, Use one channel file per acquisition run”. Then the process “Apply montage” should create a new channel file, with a new type of sensors “MONTAGE” (right-click on the channel file > Edit, to verify the type of the channels). You can change it if you want.
When computing the PSD, you must indicate the type of the channels you want to process (“MONTAGE”), or leave the option “Sensor names or types” empty