function [an, bn, period, sighat, E, amp] = utl_sincos_2d(signal,dt, varargin)


E=[];
if nargin > 2
   E = varargin{1};
end


for i=1:length(signal(1,:))
signal(:,i) = signal(:,i) - mean(signal(:,i));
end

T = length(signal(:,1));
sig(:,:)=signal;


w=2*pi/(T-1);
t=[0:T-1]';
ishift=T/2;

% Set up
% E [   cos( wt1)   cos( 2*wt1)    
%   [   cos( wt2)   cos( 2*wt1)    
%   [    ....
%   [   sin( wt1)   cos( 2*wt1)    
%   [   sin( wt2)   cos( 2*wt1)    
%   [    ....
% Dimensions of E number of time X number of frequencies
%                 number of observations X number of paramters

% x' = [ b1 b2 ...... a1 a2 .... ];

if isempty(E) == 1
  E=zeros(T,T/2*2);
  disp('Computing E ..');
icos=0;
for n=1:T/2
    icos=icos+1;
    E(:,icos) = cos(w*n*t);
    isin = ishift + icos;
    E(:,isin) = sin(w*n*t);
end
end

% x = amp in this case 

W = diag( ones(T,1)*1.0e-10 );
CME = inv(E'*E + W);
amp = CME*E'*sig;
sighat = E*amp;  % this is yhat = E*xhat

an = amp(1:ishift,:);
bn = amp(ishift+1:end,:);


% comppute the period
n=[1:T/2];
period=2*pi./(w*n)*dt;



in=1:length(period);
% I filter here 
%in = find( period/360 > 4.5 & period/360 < 5.5 ); size(in)
in = find( period/360 > 1.1);
disp('filtering');
%in = find( period/360 > 3);
%in=[];
anhat = an*0;
bnhat = bn*0;

anhat(in,:) = an(in,:);
bnhat(in,:) = bn(in,:);
amp_filter = [anhat; bnhat];

sig_filter = E*amp_filter;
sighat=sig_filter;

return