function [allfiltered] = applyfilterfftms(I, allfilter, sxm, sym);
% filter image by a bank of filters
% I: input images
% allfilter: filter bank
% threshold: filter response is set to zero if below threshold

nimage = size(I, 2);    % number of training images
norient = size(allfilter, 2);     % number of orientations
h = (size(allfilter{1}, 1)-1)/2;  % half size of filters
for (i = 1 : nimage)
    [sx, sy] = size(I{i});
    Im = zeros(sxm, sym); 
    Im(1:sx, 1:sy) = I{i}; 
    Jm{i} = Im; 
end

allfiltered = cell(nimage, norient);  % filtered images
for (i = 1 : nimage)
   fftI{i} = fft2(Jm{i}, sxm+h+h, sym+h+h);
end
for  (o = 1 : norient)
   fftf{o} = fft2(allfilter{1, o}, sxm+h+h, sym+h+h); 
end

for (i = 1 : nimage)
   tot = 0.;  [sx, sy] = size(I{i});
   fftIi = fftI{i}; 
   for (o = 1 : norient)
      fftfo = fftf{o}; 
      out = ifft2(fftIi.*fftfo); 
      filtered = out(h+1:h+sx, h+1:h+sy); 
       re = real(filtered); im = imag(filtered); 
       energy = re.*re + im.*im;  % compute the local energy 
       allfiltered{i, o} = energy; 
       tot = tot + sum(sum(energy((h+1):(sx-h-1), (h+1):(sy-h-1))))/(sx-2*h-1)/(sy-2*h-1); 
   end
   ave = tot/norient; 
   for o = 1 : norient
       allfiltered{i, o} = allfiltered{i, o}/ave;  % normalizing by whitening
   end
end
    
 
