lrkrol / SEREEGA

SEREEGA: Simulating Event-Related EEG Activity
64 stars 16 forks source link

generate_scalpdata from signal component error #22

Closed olivierroymd closed 3 months ago

olivierroymd commented 1 year ago

Hello, I was at the PracticalMEEG2022. I am trying the provided example scripts and systematically get an error while running the generate_scalpdata line concerning the signal component. Up to now, I tried scripts names code01.. to code05.. . Code01 was OK but the same errors come up at the same place in subsequent scripts. Here is the error:

Error using fir1 SWITCH expression must be a scalar or a character vector. Error in ersp_generate_signal (line 142) b = fir1(c{:}); Error in ersp_generate_signal_fromclass (line 124) signal = ersp_generate_signal(frequency, amplitude, phase, epochs.srate, epochs.length, ... Error in generate_signal_fromclass (line 73) signal = class_generate_signal_fromclass(class, epochs, 'epochNumber', epochNumber, 'baseonly', baseonly); Error in generate_signal_fromcomponent (line 72) signal = generate_signal_fromclass(component.signal{s}, epochs, 'epochNumber', epochNumber, 'baseonly', baseonly); Error in generate_scalpdata (line 131) parfor e = 1:epochs.n

On MacBook Pro M1 Max

Edit: format

Thanks

lrkrol commented 1 year ago

Dear Olivier,

This seems to be an issue with the fir1 function. Can you tell me the output of the following two commands?

which fir1
ver signal

They show where the fir1 function is located on your machine, and what version of the Signal Processing Toolbox you have installed.

olivierroymd commented 1 year ago

There it is:

/Applications/MATLAB_R2022b.app/toolbox/signal/signal/fir1.m
-----------------------------------------------------------------------------------------------------
MATLAB Version: 9.13.0.2105380 (R2022b) Update 2
MATLAB License Number: STUDENT
Operating System: macOS  Version: 13.1 Build: 22C65 
Java Version: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode
-----------------------------------------------------------------------------------------------------
Signal Processing Toolbox                             Version 9.1         (R2022b)

Thanks!

lrkrol commented 1 year ago

That all seems to be in order. However, when I install the same versions of MATALB and the Signal Processing Toolbox, I'm afraid I cannot reproduce the error you're getting: code02 runs through without errors. Are you sure you are using this exact file without any changes?

% the purpose of this script is to construct a simulation that has more
% than one source. furthermore, instead of an ERP, we will use an ERSP.

% configuring the simulation
epochs = struct( ...
        'n', 10, ...
        'srate', 1000, ...
        'length', 1000);

% loading a head model: the new york head
lf = lf_generate_fromnyhead('montage', 'BioSemi64');

% we take our ERSP to reflect a hand movement. therefore, we put its source
% in the motor cortex. (it may be easure to use the GUI to pick sources.)
% we furthermore add additional sources to achieve a total of 64 sources.
% we select these additional sources randomly, but such that each source is
% at least 10 mm from the others, including our hand-picked ersp source.
src_ersp = lf_get_source_nearest(lf, [20 0 65]);
src_all = lf_get_source_spaced(lf, 64, 10, 'sourceIdx', src_ersp);

% we now define an ERSP signal class representing event-related
% desynchronisation in the mu band, with some custom variability.
sig_ersp = struct( ...
        'type', 'ersp', ...
        'frequency', [8 9 12 13], ...
        'amplitude', 1, ...
        'modulation', 'invburst', ...
        'modLatency', 500, ...
        'modLatencyDv', 100, ...
        'modWidth', 500, ...
        'modWidthDv', 200, ...
        'modTaper', .5, ...
        'modMinRelAmplitude', .1);

% we verify and complete the signal class, and plot it to verify it
% visually.
sig_ersp = utl_check_class(sig_ersp);        
plot_signal_fromclass(sig_ersp, epochs);

% we now define a noise signal class. we use the uniform variant (-unif)
% for this tutorial simply because it is faster than the gaussian variant
% and does not rely on the DSP System Toolbox.
sig_noise = struct(...
        'type', 'noise', ...
        'color', 'brown-unif', ...
        'amplitude', 1);

% and again we assign the signals to the sources. we first assign the same
% noise signal to all 64 sources. because the first source in the src_all
% variable is our motor cortex source, the first component in the cmp
% variable will be our motor cortex component. therefore, we add the motor
% ERSP to this first component. we can again plot it to confirm this is
% correct.
cmp = utl_create_component(src_all, sig_noise, lf);
cmp(1) = utl_add_signal_tocomponent(sig_ersp, cmp(1));

plot_component(cmp(1:10), epochs, lf);

% as before, we simulate the data and use EEGLAB to analyse it.
data = generate_scalpdata(cmp, lf, epochs);
EEG = utl_create_eeglabdataset(data, epochs, lf);

pop_eegplot(EEG, 1, 0);
figure; pop_spectopo(EEG, 1, [0  999], 'EEG', 'freqrange', [2 30], 'freq', [11]);

% as you will see, we now have a dataset that already looks much more
% realistic. however, we have mixed background noise with our ERSP signal
% in a way that is difficult to predict or control: the relative amplitudes
% at the scalp of the noise and our ERSP depend simultaneously on the
% number of sources, their individual signal amplitudes, and the lead
% field. we will fix this in the next script.
olivierroymd commented 1 year ago

With a copy-paste of the above code, I get the same error. Strange.

lrkrol commented 1 year ago

Well, let's do some tests... Here is some sample code directly from MATLAB's own kaiserord documentation.

fsamp = 8000;
fcuts = [1000 1500];
mags = [1 0];
devs = [0.05 0.01];

[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');

freqz(hh,1,1024,fsamp)

Does this run for you?

And how about this?

c = kaiserord(fcuts,mags,devs,fsamp, 'cell');
hh = fir1(c{:});
olivierroymd commented 1 year ago

Works no problem. Same output for both methods.

hh =
  Columns 1 through 4
      -0.00241742493228756      -0.00307857977111449      3.06913893225434e-18       0.00546835860893794
  Columns 5 through 8
       0.00779124343076122        0.0020733853498149      -0.00934798216760109       -0.0159759611807864
  Columns 9 through 12
      -0.00763221294785892        0.0135238769719747        0.0297203708815585        0.0202508673081746
  Columns 13 through 16
       -0.0172807685397804       -0.0558722345972558       -0.0524277942694054        0.0198958188129006
  Columns 17 through 20
         0.144481692153349          0.26350827896674                    0.3125          0.26350827896674
  Columns 21 through 24
         0.144481692153349        0.0198958188129006       -0.0524277942694054       -0.0558722345972558
  Columns 25 through 28
       -0.0172807685397804        0.0202508673081746        0.0297203708815585        0.0135238769719747
  Columns 29 through 32
      -0.00763221294785892       -0.0159759611807864      -0.00934798216760109        0.0020733853498149
  Columns 33 through 36
       0.00779124343076122       0.00546835860893794      3.06913893225434e-18      -0.00307857977111449
  Column 37
      -0.00241742493228756
TylorJHarlow commented 3 months ago

Has this issue been resolved? I am getting the same error when running the copied & pasted script. However, additionally when I am trying to generate scalp data I see:

In parallel_function>make_general_channel/channel_general (line 837) In parallel.internal.parfor.cppRemoteParallelFunction (line 53) Warning: all-zero orientation for component 1

Additionally, when running:

which fir1 ver signal

I get similar results to the above:

/Applications/Toolbox/fieldtrip/external/signal/fir1.m

MATLAB Version: 23.2.0.2515942 (R2023b) Update 7 MATLAB License Number: 1144620 Operating System: macOS Version: 14.4 Build: 23E214 Java Version: Java 1.8.0_382-b05 with Amazon.com Inc. OpenJDK 64-Bit Server VM mixed mode

Signal Processing Toolbox Version 23.2 (R2023b)

I am running macOS Sonoma on an M2 Max with MATLAB 2023b.

TylorJHarlow commented 3 months ago

I was able to resolve this issue with the addition of the following line in fir1 - alternatively, replacing fir1 with the updated version from the signal processing toolbox works as well.

% sort arglist, normalize any string length(varargin) for i=1:length(varargin) arg = varargin{i}; if ischar(arg), arg=lower(arg);end if isempty(arg), continue; end
if isa(arg,'double'), window = arg; continue; end switch arg case {'low','stop','dc-1'}, ftype = 1; case {'high','pass','bandpass','dc-0'}, ftype = 0; case {'scale'}, scale = 1; case {'noscale'}, scale = 0; end end

The former issue of all-zero orientation was due to my not setting orientations for my sources. This was resolved by doing such.

lrkrol commented 3 months ago

Thanks, Tylor, indeed the issue had remained stuck at being not reproducible from my side. I see that your original fir1 stems from Fieldtrip, not the Signal Processing Toolbox, so I feel this confirms my earlier suspicion. Thanks for posting the solution you found. Assuming this would solve it for Olivier as well, I shall close this issue.