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


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

signal = signal - mean(signal);
T = length(signal);
sig=zeros(T,1);
sig(:,1)=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*wt2)    
%   [    ....
%   [   sin( wt1)   sin( 2*wt1)    
%   [   sin( wt2)   sin( 2*wt2)    
%   [    ....
% 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)
%sin = find( period/360 < 2);
%in = find( period/360 > 3);
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