function [rfto GX] = FourierFiltTO(rf, param)
%FourierFilTO   FourierFilTO produce RF image with Transverse Oscillations
%   [rfto] = FourierFiltTO(rf, param) introduce Transverse
%   Oscillations into the image stored in the array rf. rfto contains RF
%   image with transverse ocillation. sigmax and fx0 are the paramteres of
%   the double gaussian function allowing the introduction of Transverse
%   Oscillations.GX contains the calculated "mask" by the function
%   FourierFiltTO. rfto is the filtered RF image containing Transverse
%   oscillation 
% 
%
%   %   PARAM is a structure that contains the following fields:
%   -------------------------------------------------------
% 
%   1) PARAM.sigmax: The wanted lateral  full width at half maximum (FWHM) (in mm, by default = 1)                                                
%   2) PARAM.fx0: The wanted lateral wavelengths (in mm, by default = 1.5)
%
%
%   
%   Examples:
%   --------
%   % 1-D example
%   load phantom.mat
%   param.sigmax=1;
%   param.lambdax=1.5;
%   [rfto] = FourierFiltTO(migRF1(450:600,:,1), param);
%   figure(1);
%   subplot(121);imagesc(migRF1(450:600,:,1));
%   subplot(122);imagesc(rfto);
%
%   Reference
%   ---------
%
%   Sebastien salles, Simon Lai, Damien Garcia, Alfred Yu, Didier Vray,
%   Hervé Liebgott1 "Estimation of arterial wall motion using ultrafast
%   imaging with transverse oscillations" Ultrasonics Symposium (IUS),
%   2014 IEEE International 
%
%   Salles, Sebastien; Garcia, Damien; Bou-Said, Benyebka; Savary, Florent;
%   Serusclat, Andre; Vray, Didier; Liebgott, Herve, "Plane wave transverse
%   oscillation (PWTO): An ultra-fast transverse oscillation imaging mode
%   performed in the Fourier domain for 2D motion estimation of the carotid
%   artery," Biomedical Imaging (ISBI), 2014 IEEE 11th International
%   Symposium on , vol., no., pp.1409,1412, April 29 2014-May 2 2014
% 

%--- Check input parameters ---
if ~isfield(param,'sigmax')
    param.sigmax = 1.5; 
end
if ~isfield(param,'lambdax')
    param.lambdax = 1; 
end


sigmax=param.sigmax;
fx0=1/param.lambdax;
pitch=param.pitch;
[no_lignes, no_col, no_times] = size(rf);

fx = linspace(-1/(2*pitch),+1/(2*pitch),no_col);

gx =(exp(-2*(sigmax*pi*(fx-fx0)).^2) + ...
     exp(-2*(sigmax*pi*(fx+fx0)).^2));

rfto = real(ifft2(bsxfun(@times,fft2(rf),ifftshift(gx))));

 for i=1:no_lignes
     GX(i,:)=gx;
 end
 param.GX=GX;
end
