filtfilt error

[ShahinSafa]

When I use .resp log file in PhysIO, I get an error (copied below). Would you please let me know how to fix the issue? (Btw, the signal processing toolbox has been already installed on my PC and I can use the "filtfilt" function.)

Thanks in advance for your help.

Bests, Shahin

Error using tapas_physio_filter_respiratory (line 137) Usage: y=filtfilt(b,a,x) In file "/Applications/MATLAB_R2021a.app/toolbox/tapas-master/PhysIO/code/preproc/tapas_physio_filter_respiratory.m" (???), function "tapas_physio_filter_respiratory" at line 137. In file "/Applications/MATLAB_R2021a.app/toolbox/tapas-master/PhysIO/code/tapas_physio_main_create_regressors.m" (???), function "tapas_physio_main_create_regressors" at line 225. In file "/Applications/MATLAB_R2021a.app/toolbox/tapas-master/PhysIO/code/tapas_physio_cfg_matlabbatch.m" (???), function "run_physio" at line 1661.

[mrikasper]

Dear Shahin,

Thank you for using PhysIO! It might be that filtfilt is still the culprit, if you are using an older version of SPM (i.e., not the current version from the SPM GitHub) and have the SPM subfolders of FieldTrip added to your path.

Please see this issue where I describe how to use the diagnostics of PhysIO and add/remove the right paths in the right order.

A quick check would be to type which filtfilt...if it shows an SPM-related path, this would be an indicator you have the same issue.

All the best, Lars

[ShahinSafa]

Hi Lars,

Thanks for the quick respond.

The 'filtfilt" path is not within the SPM folder. It is here:

/Applications/MATLAB_R2021a.app/toolbox/fieldtrip/fieldtrip-20220113/external/signal/filtfilt.m

As I followed the error path, it comes from line 56 in the "filtfilt" script:

if (min(size(b)) > 1) warning('SOS structures not supported for FieldTrip drop-in filtfilt'); error('Usage: y=filtfilt(b,a,x)'); end

In the workspace, I see that b is a 10x6 array, while (if I understood the code properly), it is supposed to be a single value to be used as the filter order, right?

Could the .resp file be the issue? (The log files are attached).

Thanks again!

Archive.zip

[mrikasper]

Dear Shahin,

If I understand the Fieldtrip documentation correctly, the issue should have been fixed on Jan 12, 2022, so your version should be fine (see linked thread in Guillaume's comment). However, just as an easy test, could you

1. Remove the fieldtrip path you mentioned from your path and run PhysIO again (check with which that Matlab's own filtfilt is now used)
2. Update Fieldtrip to the most recent version (and make sure it's used in your path)

If neither of this works, could you please send me your batch or Matlab .m script as well to test it on your data. From looking at the raw .resp data, it seems to be fine.

All the best, Lars

[ShahinSafa]

Dear Lars,

The second method (updating the Fieldtrip and using it "filtfilt") did not work. However, the first method (removing the mentioned path and using the Matlab "filtfilt" code) works. But, can I trust the generated result not using the "filtfilt" from Fieldtrip?

Btw, I attached the batch file.

Thanks very much! Shahin

matlabbatch.m.zip

[ShahinSafa]

Hi Aagin,

The generated regressors for the respiratory are zeros! It seems that using the Matlab "filtfilt", the PhysIO tollbox does not generate accurate regressors! The generated regressors are attached!

I would be thankful if you could generate the regressors in your computer and let me know if the problem is due to the log files or the installed toolbox on my computer.

Thanks very much!

Regards, Shahin

multiple_regressors_1.txt

[mrikasper]

Dear Shahin,

just to let you know: I am working on this. I have attached a batch that works for the respiratory data (just paste it in an empty .m file). I think the reason it was all zeros was that you specified the .ecg file for the respiratory log file as well.

However, I will have to update the read-in for the ECG data, it's a new version of the Siemens logfile with 4 channels.

With respect to whether you can trust Matlab's filtfilt, yes, of course, this is the assumed function for PhysIO. I would rather be careful using the fieldtrip version (because I don't check compatibility with it upon releases), though I don't think it's the issue you are encountering.

BTW: is your TR really 10 seconds? That seems long to me.

All the best, Lars

%-----------------------------------------------------------------------
% Job saved on 27-Jul-2022 12:56:46 by cfg_util (rev $Rev: 7345 $)
% spm SPM - SPM12 (7771)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.spm.tools.physio.save_dir = {'/Physio'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Siemens';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'C:\Users\kasperla\polybox\Shared\PhysioUserData\gh204-filtfilt\run1.ecg'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'C:\Users\kasperla\polybox\Shared\PhysioUserData\gh204-filtfilt\run1.resp'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {''};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = [];
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 29;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 10;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 127;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 24;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.nominal = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'ECG';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.max_heart_rate_bpm = 90;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.respiratory.filter.passband = [0.01 2];
matlabbatch{1}.spm.tools.physio.preproc.respiratory.despike = false;
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = '[]';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

[ShahinSafa]

Dear Lars,

The respiratory file in the batch that I sent you was accidently chosen a .ecg file. However, my script selects the .resp batch correctly.

Additionally, the TR is 10 seconds because of the spare sampling mode. This is an fMRI study on the auditory stystem of the brain. So, I have to put scanning gaps when I am presenting the sound to the subject. Otherwise, the scanner noise and the presenting sound will be overlapped and the BOLD response becomes messy!

However, I think I found the problem. In your batch, the "sampling interval" is empty whereas in my batch it was 0.0025. Now, the regressors for the respiratory are not zeros anymore. I hope the the generated values could correct the physiological motion. I will test it today.

And finally, can I still use this batch for the new version of the Siemens log files? Or I have to wait for an update on PhsyIO?

Best regards, Shahin

[mrikasper]

Dear Shahin,

Regarding the sampling rate, yes, it makes sense that this created your error, because you will have to specify an array of two values for cardiac and respiratory logfiles, if the sampling rates differ (which they did at least for older logfile versions on Siemens, where ECG was sampled every 2.5 ms, but pulse and resp every 20 ms, and external TTL triggers every 5 ms). In most cases for standard scanner log files, it's best to leave it up to PhysIO to determine these rates, which is the default if your provide empty values.

The output respiratory regressors should be fine with the current version of PhysIO and the correct batch, but the cardiac regressors are wrong. However, it should only be a couple of days until I post the updated version of the reader.

All the best, Lars

[ShahinSafa]

Hi Lars,

Then, I will look forward to downloading the new version of the reader.

Thanks again for your kind support.

Best regards, Shahin

[mrikasper]

Dear Shahin,

I have now updated the reader and retrieve decent results for your example data with batch .m file pasted below:

1. I will make the code available in the physio branch of TAPAS as a beta version, probably tomorrow, and it will be included and fully tested in the next TAPAS release in September.

2. As an aside, while the ECG peak detection looks OK, there is some noise variation (probably whenever the gradients are on vs silent periods in your paradigm)

   • It seems the default choice of PhysIO to take the mean of the ECG channels is not optimal (because the first channel is very noisy). Ideally, you would take the 3rd channel or so (see screenshot)
   • There is no way to change this as an input parameter right now, but you could hard-code the variable cardiac_modality in tapas_physio_read_physlogfiles_siemens to a value v3 for the ECG case in line 71.

3. Would you be willing to share your example data as one of the standard examples shipped with PhysIO? This way, I could integrate it better in my read-in tests and make sure all logfile versions of Siemens are supported in future versions of PhysIO.

I hope that helps.

All the best, Lars

matlabbatch{1}.spm.tools.physio.save_dir = {'/Physio'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Siemens';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'run1.ecg'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'run1.resp'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {''};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = [];
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 29;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 10;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 127;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 24;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.nominal = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'ECG';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.max_heart_rate_bpm = 90;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.respiratory.filter.passband = [0.01 2];
matlabbatch{1}.spm.tools.physio.preproc.respiratory.despike = false;
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = '[]';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

[ShahinSafa]

Hi Lars,

Thank you so much for working on this problem and updating the scripts.

I edited the scripts on my computer and everything goes smoothly.

Just regarding the ECG channels, I replaced "mean" with "v3". But, it seems "v3" is not an option for this variable!

_switch cardiacmodality case 'ECG' defaults.ecgChannel = 'v3'; %'mean'; 'v1'; 'v2' otherwise defaults.ecgChannel = 'v1'; end

And sure. Feel free to use that log file as an example in the new version.

I really appreciate your help to fix this issue.

Best regards, Shahin

[mrikasper]

Dear Shahin,

Sorry, I must have overlooked this comment. Did you try setting defaults.ecgChannel='v3'; and it didn't work? Or do you mean that the option is not mentioned in the comment of that line or function header? I believe it should work if you hard-code the channel in that file, but I agree, the high-level interface of PhysIO does not support this selection yet.

Let me know if this fixes the issue for now, please!

All the best, Lars

[ShahinSafa]

Dear Lars,

Thanks a lot for looking back and replying to this email. I really appreciate your kind support.

I did hardcode the "v3" part and it runs without error. Although, I am not sure about the result. The result of the fMRI data does not seem to be changed with and without applying PhysIO. Maybe, that is what it should be.

Thanks a lot.

Bests, Shahin

On Tue, Oct 25, 2022 at 3:33 PM Lars Kasper @.***> wrote:

Dear Shahin,

Sorry, I must have overlooked this comment. Did you try setting defaults.ecgChannel='v3'; and it didn't work? Or do you mean that the option is not mentioned in the comment of that line or function header? I believe it should work if you hard-code the channel in that file, but I agree, the high-level interface of PhysIO does not support this selection yet.

Let me know if this fixes the issue for now, please!

All the best, Lars

— Reply to this email directly, view it on GitHub https://github.com/translationalneuromodeling/tapas/issues/204#issuecomment-1290565772, or unsubscribe https://github.com/notifications/unsubscribe-auth/AQGQYBCQKWVNXPQI23SULLDWE7OUJANCNFSM54RUSMLA . You are receiving this because you authored the thread.Message ID: @.***>

-- S. Safazadeh @.***