nt_dft_filter

y=nt_dft_filter(x,transfer,N) - apply filter using DFT

  y: result
  
  x: data (time*channels)
  transfer: transfer function (size=N/2)
  N: number of samples to use in DFT (default: next power of 2)

 Data are zero-padded to N.

0001 function y=nt_dft_filter(x,transfer,N)
0002 %y=nt_dft_filter(x,transfer,N) - apply filter using DFT
0003 %
0004 %  y: result
0005 %
0006 %  x: data (time*channels)
0007 %  transfer: transfer function (size=N/2)
0008 %  N: number of samples to use in DFT (default: next power of 2)
0009 %
0010 % Data are zero-padded to N.
0011 
0012 if nargin<3||isempty(N); N=2.^nextpow2(size(x,1)); end
0013 if nargin<2; error('!'); end
0014 
0015 if mod(N,2)==1; error('N must be multiple of 2'); end % for convenience
0016 
0017 % transfer function - if pair of numbers interpret as bandpass
0018 transfer=transfer(:);
0019 if numel(transfer)==2;
0020     lo=1+round(N/2*transfer(1)); 
0021     hi=1+round(N/2*transfer(2));
0022     if min(lo,hi)<1; error('band limits too low'); end
0023     if max(lo,hi)>N/2; error('band limits too high'); end
0024     transfer=zeros(N/2,1);
0025     transfer(lo:hi)=1;
0026 end
0027 
0028 [m,n,o]=size(x);
0029 x=reshape(x,m,n*o); % all trials and channels in parallel
0030 
0031 % pad with zeros
0032 x=[x;zeros(N-m,n*o)];
0033 
0034 % DFT, apply transfer function, inverse DFT
0035 y=fft(x);
0036 transfer=repmat(transfer,1,n*o);
0037 size(transfer);
0038 y=y.*[transfer;flipud(transfer)];
0039 y=real(ifft(y));
0040 
0041 % trim, reshape to original geometry
0042 y=y(1:m,:);
0043 y=reshape(y,[m,n,o]);
0044 
0045