[FieldTrip] Frequency shift after source reconstruction

Wed May 27 09:47:30 CEST 2020

Caro Luca,

I haven’t read your code in detail, but I suspect that you perform a spectral analysis on the source.avg.pow, which is the magnitude(-squared) version of the signal’s time course.
As a consequence of this, your frequency ‘doubles’. Consider a sine wave with a given frequency (i.e. a period of 1/f), say 1 Hz. If you take that absolute(-squared) of that sine wave, you get a signal with ’two bumps’ per second, which reflects a periodicity of 2 Hz.

You should work with the ‘mom’ field at the source level.

Best wishes,

Jan-Mathijs

On 27 May 2020, at 09:33, Luca La Fisca <534843 at umons.ac.be<mailto:534843 at umons.ac.be>> wrote:

Dear community,

I'm Luca La Fisca, PhD student in Belgium and I'm validating the source reconstruction algorithm using "ft_dipolesimulation".

After having performed a frequency analysis, I noticed the reconstructed sources have twice the initial frequencies as shown in the following figures (initial dipole frequencies = 10 and 85 Hz) :

Timelock:
<image.png>

Reconstructed sources:
<image.png>

Here is the code used to get the timelock signal:
fs = 2048;
fs_down = 512;

% generate EEG pseudo-signal
cfg      = [];
cfg.headmodel = vol;
cfg.elec = elec_aligned;
cfg.channel = elec_aligned.label(1:64);

cfg.dip.pos = [
   -1 65 36           % dipole left
   1 -62 43          % dipole right
   ];

unif_dip = 3^(1/2)/3;
cfg.dip.mom = ...
  [ unif_dip unif_dip unif_dip 0 0 0 ]' + ...
  [ 0 0 0 unif_dip unif_dip unif_dip ]';

event = ft_read_event(EEG_FILE);
event = event(3:end);
signal1 = zeros(length(eeg_data.trial));
signal2 = zeros(length(eeg_data.trial));
val = [event.value];
for s = drange([event(val==10).sample, event(val==11).sample, event(val==12).sample])
    time = (0:ceil(fs*0.5))/fs;
    signal1(s:s+ceil(fs*0.5)) = 10*sin(10*time*2*pi);
    signal2(s:s+ceil(fs*0.5)) = 30*cos(85*time*2*pi);
end

cfg.dip.signal = {[signal1; signal2]};
cfg.fsample = fs;
raw1 = ft_dipolesimulation(cfg);

% Timelock analysis
% define trials
cfg = [];
cfg.dataset = EEG_FILE;
cfg.trialdef.eventtype  = 'STATUS';
cfg.trialdef.eventvalue = [10 11 12];
cfg.trialdef.prestim    = 0.5;
cfg.trialdef.poststim   = 1;
cfg = ft_definetrial(cfg);
data_trial = ft_redefinetrial(cfg, raw1);

% downsample
cfg = [];
cfg.resamplefs = fs_down;
data_trial = ft_resampledata(cfg,data_trial);

cfg = [];
timelock = ft_timelockanalysis(cfg, data_trial);

Here the code related to the source reconstruction:
% Estimating the leadfield
cfg = [];
cfg.sourcemodel = sourcemodel;    %% where are the sources?
cfg.headmodel   = vol;      %% how do currents spread?
cfg.elec        = elec_aligned; %% where are the sensors?
cfg.channel = elec_aligned.label(1:64);
sourcemodel_and_leadfield = ft_prepare_leadfield(cfg, timelock);

% Source reconstruction
cfg                     = [];
cfg.method              = 'mne';
cfg.sourcemodel         = sourcemodel_and_leadfield;        %the precomputed leadfield
cfg.headmodel           = vol;                      %the head model
cfg.elec                = elec_aligned;             %the electrodes
cfg.channel             = elec_aligned.label(1:64); %the useful channels
cfg.mne.prewhiten       = 'yes';                    %prewhiten data
cfg.mne.lambda          = 0.01;                     %regularisation parameter
cfg.mne.scalesourcecov  = 'yes';                %scaling the source covariance matrix
minimum_norm_eeg        = ft_sourceanalysis(cfg,timelock);

Here the code to perform the frequency analysis:
tmp_pow = minimum_norm_eeg.avg.pow;
pow_fft = fft(tmp_pow'); %freq analysis of reconstructed sources
% pow_fft = fft(timelock.avg'); %freq analysis of timelock signal
% (un)comment depending on the desired one

L = length(timelock.avg);
f = fs_down*(0:(L/2))/L;

P2 = abs(pow_fft/L);
P1 = P2(1:L/2+1,:);
P1 = P1.*f';

figure()
plot(f,P1)
title('Single-Sided Amplitude Spectrum of Reconstructed Sources')
xlabel('f (Hz)')
ylabel('|P1(f)|')

Do you have any idea of what is responsible to this double frequency?
Do you think it is a general issue?
Is it a major issue leading to wrong analysis of the reconstructed signal (e.g. connectivity analysis)?

I thank you in advance.

Best regards,

[https://docs.google.com/uc?export=download&id=1kIixsLbG6VMY0eKRmD_V2qGfWdocai3Y&revid=0B-dw3epKQGxwWWNSYThqWVJCZjZYOGJtNnljTHg1OGR5WFNJPQ]
Luca La Fisca
PhD Student
Service d'Information, Signal et Intelligence Artificielle
Boulevard Dolez, 31
7000 Mons
+32 (0)65/37.40.83
Luca.LAFISCA at umons.ac.be<mailto:Luca.LAFISCA at umons.ac.be>

[https://docs.google.com/uc?export=download&id=1Os8PgOWiHRZwQs_oiXqKDg2gkH2mcJDb&revid=0B-dw3epKQGxweUxueTZMZzRNd2pYR1docldnQzJzWExCWGpFPQ]

[https://ipmcdn.avast.com/images/icons/icon-envelope-tick-round-orange-animated-no-repeat-v1.gif]<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail>  Garanti sans virus. www.avast.com<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail>
<x-msg://5/#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://doi.org/10.1371/journal.pcbi.1002202

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20200527/e82236f7/attachment.htm>