Back to the main page.
Bug 1408 - make the demeaning of data prior to filtering consistent
| Status | CLOSED FIXED |
| Reported | 2012-04-04 15:49:00 +0200 |
| Modified | 2012-10-10 13:27:32 +0200 |
| Product: | FieldTrip |
| Component: | preproc |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P4 minor |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Robert Oostenveld - 2012-04-04 15:49:13 +0200
as discussed in the FT meeting: preproc calls demean/detrend/polyremoval both before and after the filtering. There are good reasons for it, but the code can (and should) be cleaned up 1) it is desirable to have all ft_preproc_xxxfilter functions demean and (if applicable) re-mean the data 1b) should figure out what the fastest way is of adding a constant to all channels 2) the preproc code should go like this if polyremove polyremove with the desired order elseif detrend polyremove with order 1 elseif demean polyremove with order 0 end and this should happen for the whole data (not just the baseline, and not just the data minus padding).
Robert Oostenveld - 2012-04-04 15:49:30 +0200
also revisit bug #780
Roemer van der Meij - 2012-04-04 15:55:00 +0200
I emailed this for fun to JM, is useful here: The non-for-loops use maximally around 9,2GB, the for-loops around 4,8GB (difference between the first two depends on the order and other things, over several runs they were the same) clear; dat = rand(200,3e6); datmean = rand(200,1); tic; dat = dat - datmean(:,ones(1,3e6));toc Elapsed time is 8.134499 seconds. clear; dat = rand(200,3e6); datmean = rand(200,1); tic; dat = dat - repmat(datmean,[1 3e6]);toc Elapsed time is 7.630849 seconds. clear; dat = rand(200,3e6); datmean = rand(200,1); tic for ichan = 1:200 dat(ichan,:) = dat(ichan,:) - datmean(1,ones(1,3e6)); end toc Elapsed time is 33.508252 seconds. clear; dat = rand(200,3e6); datmean = rand(200,1); tic for ichan = 1:200 dat(ichan,:) = dat(ichan,:) - repmat(datmean(1),[1 3e6]); end toc Elapsed time is 24.960733 seconds.
Jan-Mathijs Schoffelen - 2012-04-04 15:56:20 +0200
I replied this to Roemer: Een betere vergelijking zou zijn: clear dat = rand(200,3e6); datmean = rand(200,1); tic for ichan = 1:200 dat(ichan,:) = dat(ichan,:) - datmean(ichan); end toc clear dat = rand(200,3e6); datmean = rand(200,1); tic for ichan = 1:3e6 dat(:,ichan) = dat(:,ichan) - datmean; end toc
Roemer van der Meij - 2012-04-04 16:01:58 +0200
My third and fourth were of course rubbish, I was sleeping at the moment. The ones JM mentioned: clear; dat = rand(200,3e6); datmean = rand(200,1); tic for ichan = 1:200 dat(ichan,:) = dat(ichan,:) - datmean(ichan); end toc Elapsed time is 17.273900 seconds. clear; dat = rand(200,3e6); datmean = rand(200,1); tic dat = dat'; for ichan = 1:200 dat(:,ichan) = dat(:,ichan) - datmean(ichan); end dat = dat'; toc Elapsed time is 15.349796 seconds. The last one, with the transpose, uses 9,2GB of course, the first 4,8GB of course. The GB is total memory use, the data matrix is 4,8GB.
Jan-Mathijs Schoffelen - 2012-04-04 16:05:08 +0200
still sleeping, indexing over the columns should run between 1 and 3e6 (rather than 200). Also, in that case you don't need to index datmean.
Robert Oostenveld - 2012-04-04 16:10:25 +0200
With a slightly smaller problem size I get >> nchans = 200; >> nsamples = 1e6; >> dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic dat = bsxfun(@minus, dat, datmean); toc Elapsed time is 2.418713 seconds. >> dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic dat = dat - datmean(:,ones(1,nsamples)); toc Elapsed time is 3.118488 seconds. so bsxfun is the winner sofar.
Robert Oostenveld - 2012-04-04 16:16:10 +0200
(In reply to comment #6) > so bsxfun is the winner sofar dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic dat = bsxfun(@minus, dat, datmean); toc Elapsed time is 2.365440 seconds. dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic for isample = 1:nsamples dat(:,isample) = dat(:, isample) - datmean; end toc Elapsed time is 1.118270 seconds. Ha, found an even faster one. The penalty of the many-more-iterations loop over 1e6 samples is still outweighed by the most efficient memory access! Note that my problem size is 3x smaller, so to make it a better comparison here is the one Roemer started with... dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic dat = dat - datmean(:,ones(1,nsamples)); toc Elapsed time is 3.155108 seconds
Robert Oostenveld - 2012-04-04 16:22:31 +0200
To wrap up: I made this (which I'll add as test script) ------------------------------------------------------- nchans = 200; nsamples = 1e6; dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic dat = dat - datmean(:,ones(1,nsamples)); t(1) = toc dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic dat = dat - repmat(datmean,[1 nsamples]); t(2) = toc dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic for ichan = 1:nchans dat(ichan,:) = dat(ichan,:) - datmean(ichan); end t(3) = toc dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic dat = dat'; for ichan = 1:nchans dat(:,ichan) = dat(:,ichan) - datmean(ichan); end dat = dat'; t(4) = toc dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic dat = bsxfun(@minus, dat, datmean); t(5) = toc dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic for isample = 1:nsamples dat(:,isample) = dat(:, isample) - datmean; end t(6) = toc [minval, minindx] = min(t); if minindx~=6 error('unexpected winner of the speed test'); end ------------------------------------------------------- which results in t = 3.0803 3.1377 7.9075 8.8355 2.3918 1.1138
Robert Oostenveld - 2012-04-04 16:43:02 +0200
I found this old code (commented out) in ft_preproc_baselinecorrect % save this old code because it looks non-trivially optimized % % persistent hasbsxfun % if isempty(hasbsxfun) % hasbsxfun = exist('bsxfun', 'builtin')==5; % end % % % determine the size of the data % [Nchans, Nsamples] = size(dat); % % % determine the interval to use for baseline correction % if nargin<2 || isempty(begsample) % begsample = 1; % end % if nargin<3 || isempty(endsample) % endsample = Nsamples; % end % % % estimate the baseline and subtract it % baseline = mean(dat(:,begsample:endsample), 2); % % % ensure that the data is not represented as integer, otherwise "minus" fails % dat = double(dat); % % if hasbsxfun % % it is even faster to do this % dat = bsxfun(@minus,dat,baseline); % else % % it is faster to loop over samples than over channels due to the internal memory representation of Matlab % % for chan=1:Nchans % % dat(chan,:) = dat(chan,:) - baseline(chan); % % end % % for sample=1:Nsamples % dat(:,sample) = dat(:,sample) - baseline; % end % end It suggests that bsxfun is the fastest, which is inconsistent with my findings today. Perhaps it is Matlab+OS version dependent...
Robert Oostenveld - 2012-04-04 16:53:10 +0200
I finished part 2, i.e. the consistent handling in preproc. Part 1 still needs to be done, the answer to 1b is known. manzana> svn commit private/preproc.m preproc/ test/test_bug1408.m Sending preproc/ft_preproc_baselinecorrect.m Sending preproc/ft_preproc_detrend.m Sending preproc/ft_preproc_polyremoval.m Sending private/preproc.m Adding test/test_bug1408.m Transmitting file data ..... Committed revision 5591.
Robert Oostenveld - 2012-04-04 17:18:40 +0200
there are >> ls *filter* ft_preproc_bandpassfilter.m ft_preproc_highpassfilter.m ft_preproc_online_filter_apply.m ft_preproc_bandstopfilter.m ft_preproc_lowpassfilter.m ft_preproc_online_filter_init.m ft_preproc_dftfilter.m ft_preproc_medianfilter.m The online filters don't have to be done, nor does the medianfilter. I have added a consistent demean/remean section to each of the other functions and added a few lines to the test script. manzana> svn commit preproc/ test/test_bug1408.m Sending preproc/ft_preproc_bandpassfilter.m Sending preproc/ft_preproc_bandstopfilter.m Sending preproc/ft_preproc_dftfilter.m Sending preproc/ft_preproc_highpassfilter.m Sending preproc/ft_preproc_lowpassfilter.m Sending preproc/ft_preproc_medianfilter.m Sending test/test_bug1408.m Transmitting file data ....... Committed revision 5593.
Roemer van der Meij - 2012-04-04 17:37:56 +0200
I am completely surprised and stunned by this: ******** dat = rand(nchans,nsamples); datmean = rand(nchans,1); tic for isample = 1:nsamples dat(:,isample) = dat(:, isample) - datmean; end toc Elapsed time is 1.118270 seconds. Ha, found an even faster one. The penalty of the many-more-iterations loop over 1e6 samples is still outweighed by the most efficient memory access! ******** That is 1 million loops! Something I would never even think of doing. And 1 million samples at 2000Hz is still 500 seconds of 200-channel data! A simulated recording of more than 8(!) hours. And decreased computation time is not insignificant either.
Robert Oostenveld - 2012-04-04 18:07:35 +0200
(In reply to comment #12) I think that the performance on this one depends quite a bit on the MATLAB version. Recent versions (I was using 2010b here) are considerably more efficient in for loops than the older versions. That is also why the bsxfun might have given better performance in the past. So an optimal solution might be if matlabversion<7.x do the one else do the other optimize solution end I leave that as exercise for the reader...
Robert Oostenveld - 2012-04-11 16:48:37 +0200
I cleaned up my bugzilla list by changing the status from resolved (either fixed or wontfix) into closed. If you don't agree, please reopen the bug. Robert