Back to the main page.
Bug 2931 - Problem with antialiasing filter in resample
| Status | CLOSED FIXED |
| Reported | 2015-07-15 16:24:00 +0200 |
| Modified | 2019-08-10 12:30:49 +0200 |
| Product: | FieldTrip |
| Component: | preproc |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Windows |
| Importance: | P5 normal |
| Assigned to: | |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Vladimir Litvak - 2015-07-15 16:24:11 +0200
Created attachment 721 SPM resampling code I got the report below from one of our SPM users. Until now we haven't used ft_preproc_resample but there resample sigproc function is used as well and I think it's the default. One idea we discussed with Guillaume is that if our in-house downsampling function which uses fft can be added to FT as an option we can then call ft_preproc_resample in SPM code and give all the choices to the user so that people who are worried about this behaviour of resample can choose another option. Adding our code will be good for FT as well because all the present options require sigproc toolbox. I attach the code to this bug. Vladimir -------- I am writing to you because Alain (cc’ed) brought an issue to my attention, concerning the way Matlab handles downsampling. From the following bit of code, am I correct in assuming that SPM uses Matlab’s resample() to downsample? if flag_tbx % Signal Proc. Toolbox Dnew(blkchan,:) = resample(Dtemp', P, Q)'; else Dnew(blkchan,:) = spm_timeseries_resample(Dtemp,P/Q); end The problem is that resample() does actually not have an appropriate antialiasing filter. Here is a snippet of code that illustrates this. I create 100 trials of a random signal at 600Hz to which I add a few sinusoids. I then downsample the signal using both resample() and decimate(). You can see that the peaks just above Nyquist frequency get folded back into resample()’s downsampled signal. clear close all fs = 600; % sampling frequency of 1kHz Ntrials=50; % Create random data x = rand(Ntrials,3*fs) ; % Add a few peaks below and above new Nyquist frequency x = x+repmat(.2*sin(2*pi*50*[1:3*fs]/fs),Ntrials,1); x = x+repmat(.5*sin(2*pi*60*[1:3*fs]/fs),Ntrials,1); x = x+repmat(.6*sin(2*pi*160*[1:3*fs]/fs),Ntrials,1); x = x+repmat(.8*sin(2*pi*220*[1:3*fs]/fs),Ntrials,1); x = x+repmat(1*sin(2*pi*155*[1:3*fs]/fs),Ntrials,1); x = [zeros(Ntrials,3*fs) x zeros(Ntrials,3*fs)]; % pad with zeros x=x'; % time* channels y = resample(x,fs/2,fs); z=[]; % decimate does not work for multidimensional arrays for i=1:size(x,2) z(:,i) = decimate(x(:,i),2); end figure; % psd plots subplot(3,1,1) pwelch(mean(x,2),[],[],[],fs); legend('original') subplot(3,1,2); pwelch(mean(y,2),[],[],[],fs/2); legend('resample') subplot(3,1,3); pwelch(mean(z,2),[],[],[],fs/2); legend('decimate') This can have serious implications for EEG/MEG data, especially for those looking at high gamma effects.
Eelke Spaak - 2015-07-15 16:33:53 +0200
The results for the test script indeed look pretty bad. I'd never tested Matlab's resample() before, but just trusted the documentation which explicitly states "resample applies an antialiasing FIR lowpass filter".
Vladimir Litvak - 2015-07-15 16:41:24 +0200
The guys who found this out suggested adding this code to reduce the cutoff frequency for resample [p2,q2] = rat( P/Q, 1e-12 ); pqmax = max(p2,q2); bta=5; fc = 1/2/pqmax; FCFACTOR=0.8; fc=fc*FCFACTOR; L = 2*10*pqmax + 1; h = p2*firls( L-1, [0 2*fc 2*fc 1], [1 1 0 0]).*kaiser(L,bta)' ; y = resample(x,P,Q,h) It might be a good idea to add it to FT code for the resample option. Vladimir
Vladimir Litvak - 2015-07-15 16:42:13 +0200
See also this related discussion http://www.auditory.org/mhonarc/2014/msg00420.html
Jan-Mathijs Schoffelen - 2015-08-07 15:20:15 +0200
what about adding an optional lowpassfiltering step (using Alex Widmann's firws filter), with an appropriate cutoff frequency? I would be reluctant on having the fourth input argument to the resample-function, which looks a bit magical to me (and moreover it uses firls, which is strongly discouraged by Alex).
Vladimir Litvak - 2015-08-08 06:55:38 +0200
Yes, this could be done as well. But adding our code which doesn't use the sigproc toolbox will also be useful.
Jan-Mathijs Schoffelen - 2015-08-09 20:10:35 +0200
Fair point. I was not aware that the resampling functions are signal processing toolbox functions. It would anyhow be good to have a single function that is doing the resampling. As far as I am aware ft_preproc_resample is not used in FT, because the downsampling runs through ft_resampledata, rather than through preproc.m which is the only function that calls it. I am perfectly fine with adding an optional method that is achieving the resampling through an fft, and back again. My feeling would be to get rid of ft_preproc_resample, and implement all functionality in a clean way in ft_resampledata. Unless you guys would prefer to keep ft_preproc_resample.
Vladimir Litvak - 2015-08-09 20:44:10 +0200
(In reply to Jan-Mathijs Schoffelen from comment #6) As I understand the code philosophy, SPM should call the low-level module functions rather than high-level functions. So if ft_resampledata uses FT data structures rather than generic matrices, it would not conform to the present standards. What would make sense is that ft_resampledata calls ft_preproc_resample and spm_eeg_downsample calls it as well. Vladimir
Jan-Mathijs Schoffelen - 2015-08-09 21:04:05 +0200
I was afraid of that. This is fine with me.
Vladimir Litvak - 2015-11-20 20:14:29 +0100
(In reply to Jan-Mathijs Schoffelen from comment #8) It looks like the last comment meant that you expected me to implement it, but I thought that someone at the Donders would. Anyway, I did it now and committed. Should not affect FT too much as you don't even use this function.