function filter1d(action,varargin)
%   1D window-based finite impluse response filter
%
%   This is a digital symetric non-recursive filter
%   which can be Lowpass, Highpass, Bandpass or Bandstop
%   by using a windowed sinc function.
%   Non-Recursive filters are based on convolving 
%   the input series with a weight function.
%   Their advantage compared to recursive filters
%   is that they are always stable and matematically simpler.
%   Their disadvantage is that they may require a large
%   number of weights to achieve a desired response.
%   Some of these windows are available with fir1 like
%   hamming but others are not. 
%   
%   Additional tools include:
%       1. Outlier detection
%       2. Manual deletion of points
%       3. Time-series calculator 
%           Allows the user to modify the independant (t) 
%           and dependant (x) variables
%       4. Interpolation and resampling
% Version: 2.0
% Date: November 1, 2012
% Author: Alejandro Sanchez, USACE-CHL
if nargin==0
    Initialize
    return
end
feval(action,varargin{:}); 
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Initialize
%default parameters
ud.pathfilter1d = cd;
ud.headerlines = 6;
ud.fc = 6;
ud.col = 3;
ud.dt = 1/100;
ud.len = 3*fix(1/ud.fc/ud.dt); %minimum
ud.fo = 1;
ud.t = [];
ud.x = [];
ud.dtf = [];
ud.tf = [];
ud.xf = [];
ud.maxgap = [];
ud.win = 'Hanning';
ud.type = 'Low';
ud.mode = 'ASCII';
ud.name = '';
ud.curvetype = '';
ud.fmt = 't z';
ud.maxdt = 3;
ud.fin = '';
ud.pathin = '';
ud.fout = '';
ud.pathout = '';
bckgndcol = 1*[1,1,1];
sizefont = 10;
sizefont2 = 8;
sizefont3 = 12;
x1 = 0.03;

% screensize = get(0,'ScreenSize');
% figpos = screensize*0.8 + [0.1*screensize(3),0.1*screensize(4),0,0];
% figpos = [172   147   692   479];
%figpos = [0.2    0.24    0.58    0.5];
figpos = [0.27    0.28    0.57    0.56];
ud.hFig = figure('Name','filter1d','Units','Normalized',...
    'NumberTitle','Off','Visible','Off',...
    'Menubar','figure','Toolbar','figure',... 
    'HandleVisibility','on','Position',figpos, ...
    'DoubleBuffer','On','color',bckgndcol,'Interruptible','off');

%Input File
uipanel('Title','Input File','FontSize',12,...
    'BackgroundColor',bckgndcol,...
    'Position',[x1 0.88 0.44 0.11]);
% uicontrol('Style','Text','String','Input File',...
%     'FontSize',sizefont,'BackgroundColor',bckgndcol,...
%     'Units','Normalized','Position',[0.02 0.91 0.15 0.05]);
ud.hfin = uicontrol('Style','Text','String','-',...
    'FontSize',sizefont3,'BackgroundColor','w',...
    'Units','Normalized','Position',[x1+0.01, 0.9 0.28 0.04]);
uicontrol('Style','pushbutton','String','Browse',...
    'FontSize',sizefont2,...
    'Units','normalized','Position',[x1+0.30, 0.895 0.12 0.05],...
    'Callback','filter1d(''getdata'')'); 

%============ Filter ===============================
uipanel('Title','Windowed Sinc Filter','FontSize',12,...
    'BackgroundColor',bckgndcol,...
    'Position',[x1 0.62 0.44 0.25]);
%Windows
uicontrol('Style','Text','String','Window',...
    'FontSize',sizefont,'BackgroundColor',bckgndcol,...
    'Units','Normalized','Position',[x1+0.01, 0.78 0.10 0.05]);
windows = 'Hanning|Hamming|Lanczos|Blackman|Gauss|Welch|Boxcar';
ud.hwin = uicontrol('Style','Popup','String',windows,...
    'Units','Normalized','Position',[x1+0.01, 0.73 0.10 0.06],...
    'BackgroundColor','w','FontSize',sizefont2,...
    'Enable','On','Callback','filter1d(''SetWindow'')');

%Filter Types
uicontrol('Style','Text','String','Type',...
    'FontSize',sizefont,'BackgroundColor',bckgndcol,...
    'Units','Normalized','Position',[x1+0.16, 0.78 0.10 0.05]);
TL = 'Low-Pass|High-Pass|Band-Pass|Band-Stop';
ud.htype = uicontrol('Style','Popup','String',TL,...
    'Units','Normalized','Position',[x1+0.16, 0.73 0.10 0.06],...
    'BackgroundColor','w','FontSize',sizefont2,...
    'Enable','On','Callback','filter1d(''SetType'')');

%Analize
uicontrol('Style','pushbutton','String','Apply',...
    'Fontweight','Normal','Fontsize',sizefont,...
    'Units','normalized','Position',[x1+0.30, 0.74 0.12 0.06],...
    'Callback','filter1d(''Analize'')'); 

%Cut-Off frequency fc
uicontrol('Style','Text','String','Cut-Off Freq.',...
    'FontSize',sizefont,'BackgroundColor',bckgndcol,...
    'Horizontalalignment','center',...
    'Units','Normalized','Position',[x1+0.01, 0.68 0.12 0.05]);
ud.hfc = uicontrol('Style','Edit','String',num2str(ud.fc),...
    'FontSize',sizefont2,'BackgroundColor','w','Enable','On',...
    'Units','Normalized','Position',[x1+0.01, 0.64 0.1 0.05],...
    'Callback','filter1d(''Setfc'')');

%Frequency width
uicontrol('Style','Text','String','Width',...
    'FontSize',sizefont,'BackgroundColor',bckgndcol,...
    'Horizontalalignment','center',...
    'Units','Normalized','Position',[x1+0.16, 0.68 0.10 0.05]);
ud.hfo = uicontrol('Style','Edit','String',num2str(ud.fc),...
    'FontSize',sizefont2,'BackgroundColor','w','Enable','On',...
    'Units','Normalized','Position',[x1+0.16, 0.64 0.1 0.05],...
    'Callback','filter1d(''Setfo'')');

set(ud.hfo,'Enable','Off','BackgroundColor',0.8*[1,1,1])

%Filter Length
uicontrol('Style','Text','String','Length',...
    'FontSize',sizefont,'BackgroundColor',bckgndcol,...
    'Horizontalalignment','center',...
    'Units','Normalized','Position',[x1+0.31, 0.68 0.10 0.05]);
ud.hLen = uicontrol('Style','Edit','String',num2str(ud.len),...
    'FontSize',sizefont2,'BackgroundColor','w','Enable','On',...
    'Units','Normalized','Position',[x1+0.31, 0.64 0.1 0.05],...
    'Callback','filter1d(''SetLen'')');

%=========== Output =========================
uipanel('Title','Output File','FontSize',12,...
    'BackgroundColor',bckgndcol,...
    'Position',[x1 0.50 0.44 0.11]);
%Output options
% uicontrol('Style','Text','String','Output File',...
%     'FontSize',sizefont,'BackgroundColor',bckgndcol,...
%     'Units','Normalized','Position',[0.05 0.48 0.15 0.05]);
% ud.hfout = uicontrol('Style','Edit','String','-',...
%     'FontSize',sizefont2,'BackgroundColor','w',...
%     'Units','Normalized','Position',[0.02 0.42 0.31 0.05]);
ud.hfout = uicontrol('Style','Text','String','-',...
    'FontSize',sizefont3,'BackgroundColor','w',...
    'Units','Normalized','Position',[x1+0.01, 0.52 0.31 0.04]);
uicontrol('Style','pushbutton','String','Save',...
    'FontSize',sizefont2,...
    'Units','normalized','Position',[x1+0.33, 0.515 0.10 0.05],...
    'Callback','filter1d(''savedata'')'); 

%Starting Time 
% uicontrol('Style','Text','String','t0',...
%     'FontSize',sizefont,'BackgroundColor',bckgndcol,...
%     'Horizontalalignment','center',...
%     'Units','Normalized','Position',[0.02 0.46 0.10 0.05]);
% ud.hfo = uicontrol('Style','Edit','String',num2str(ud.t(1)),...
%     'FontSize',sizefont2,'BackgroundColor','w','Enable','On',...
%     'Units','Normalized','Position',[0.02 0.41 0.1 0.05],...
%     'Callback','filter1d(''Setfo'')');

%End time
% uicontrol('Style','Text','String','tn',...
%     'FontSize',sizefont,'BackgroundColor',bckgndcol,...
%     'Horizontalalignment','center',...
%     'Units','Normalized','Position',[0.17 0.46 0.10 0.05]);
% ud.hfo = uicontrol('Style','Edit','String',num2str(ud.t(1)),...
%     'FontSize',sizefont2,'BackgroundColor','w','Enable','On',...
%     'Units','Normalized','Position',[0.17 0.41 0.1 0.05],...
%     'Callback','filter1d(''Setfo'')');

%Offset
% uicontrol('Style','Text','String','dz',...
%     'FontSize',sizefont,'BackgroundColor',bckgndcol,...
%     'Horizontalalignment','center',...
%     'Units','Normalized','Position',[0.32 0.46 0.1 0.05]);
% ud.hLen = uicontrol('Style','Edit','String',num2str(ud.len),...
%     'FontSize',sizefont2,'BackgroundColor','w','Enable','On',...
%     'Units','Normalized','Position',[0.32 0.41 0.1 0.05],...
%     'Callback','filter1d(''SetLen'')');

%=================================================
% uimenu('Label','File');

%=================================================
hEdit = uimenu('Label','F1D-Edit');
uimenu(hEdit,'Label','Delete Points',...
    'Callback','filter1d(''DeletePoints'')');
ud.hRestore = uimenu(hEdit,'Label','Restore Dataset',...
    'Enable','Off','Callback','filter1d(''Restore'')');

%=================================================
hView = uimenu('Label','F1D-View');
uimenu(hView,'Label','X-Label', ...
     'Checked','off','Callback','filter1d(''setxlabel'')');
uimenu(hView,'Label','Y-Label', ...
     'Checked','off','Callback','filter1d(''setylabel'')');
ud.hLineMark = uimenu(hView,'Label','Line Markers', ...
    'Checked','off','Callback','filter1d(''LineMarkers'')');

% hTS = uimenu(hView,'Label','Time-Series Plot');
% uimenu(hTS,'Label','X-Label', ...
%      'Checked','off','Callback','filter1d(''setxlabel'')');
% uimenu(hTS,'Label','Y-Label', ...
%      'Checked','off','Callback','filter1d(''setylabel'')');
% ud.hLineMark = uimenu(hTS,'Label','Line Markers', ...
%     'Checked','off','Callback','filter1d(''LineMarkers'')');
% hPR = uimenu(hView,'Label','Spectra-Response Plot');
% hPR_XAxes = uimenu(hPR,'Label','X-Axes');
% hPR_YAxes = uimenu(hPR,'Label','Y-Axes');
% hPR_XAxes = uimenu(hPR,'Label','X-Axes');
% hPR_YAxes = uimenu(hPR,'Label','Y-Axes');

%=================================================
hExtras = uimenu('Label','F1D-Tools');
uimenu(hExtras,'Label','Find Outliers',...
    'Separator','On','Callback','filter1d(''FindOutliers'')');
ud.hOutliers = uimenu(hExtras,'Label','Clear Outlier Markers',...
    'Enable','Off',...
    'Callback','filter1d(''ClearOutMark'')');
ud.hResInt = uimenu(hExtras,'Label','Resample/Interpolate',...
    'Separator','On','Enable','Off',...
    'Callback','filter1d(''ResampleInterpolate'')');
uimenu(hExtras,'Label','Calculator',...
    'Callback','filter1d(''Calculator'')');

%===================================================
hHelp = uimenu('Label','F1D-Help');
uimenu(hHelp,'Label','About filter1d',...
    'Callback','filter1d(''AboutFilter1d'')');
uimenu(hHelp,'Label','Documentation',...
    'Callback','filter1d(''Documentation'')');
uimenu(hHelp,'Label','Disclaimer',...
    'Callback','filter1d(''Disclaimer'')');

%[ud.hRestore,ud.hOutliers,ud.hResInt]

%==== Initialize variables ===========
ud.t = NaN;
ud.x = NaN;
out.f = NaN;
out.a  = NaN;

%===== Make Plots =================
pos1 = [0.08 0.1 0.85 0.35];
pos3 = [0.54 0.58  0.39 0.35];

%---- Time-series ------------
ud.ax(1) = axes('position',pos1,'XMinorTick','On');
ud.horig = plot(ud.ax(1),ud.t,ud.x,'b-');
hold on
ud.hfilt = plot(ud.ax(1),NaN,NaN,'r-');
ud.hHi = plot(ud.ax(1),NaN,NaN,'ko','Markerfacecolor','g');
ud.hout = plot(ud.ax(1),NaN,NaN,'cs');
hold off
xlabel('Elapsed Time')
set(ud.ax(1),'XMinorTick','On')

%----- Filter response curve -----------
ud.ax(2) = axes('position',pos3);
ud.hRes = semilogx(ud.ax(2),NaN,NaN,'g');
hold on
ud.hCut = plot(ud.ax(2),NaN,NaN,'m--');
hold off
xlabel(ud.ax(2),'Frequency')
ylabel(ud.ax(2),'Response')
set(ud.ax(2),'YAxisLocation','Right')

%------ Spectra --------------------
ud.ax(3) = axes('position',pos3);
ud.hSpec = semilogx(ud.ax(3),NaN,NaN,'b');
hold on
ud.hSpecf = semilogx(ud.ax(3),NaN,NaN,'r');
hold off
ylabel(ud.ax(3),'Amplitude')
set(ud.ax(3),'Color','none','box','off')

set(ud.hFig,'Visible','On')

%=== Clears file with name of output files ===========
%so that everytime filter1d starts it has zero output files
fout = 'filter1d_outputfiles.txt';
fid = fopen(fout,'w');
fclose(fid);

%===== Obtain initial data =========
fin = 'filter1d_inputfile.txt';
s = -1;
if exist(fin,'file')
    fid = fopen(fin,'r');
    s = fgetl(fid);
    fclose(fid);
end
if isequal(s,-1)
    ud = initial_timeseries(ud); %sinthetic time series
    set(ud.horig,'xdata',ud.t,'ydata',ud.x) 
    %Input Amplitude Spectra
    [f,a] = ampspectra(ud.x,ud.dt);
    set(ud.hSpec,'Xdata',f,'Ydata',a)
else
    ud = getdata(s,ud);
    [fpath,fname] = fileparts(s);
    ud.fin = fname;
    ud.pathin = fpath;
end
ud.to = ud.t;
ud.xo = ud.x;

if any(isfinite(ud.x))
    if isempty(ud.name)
        legend(ud.ax(1),'Original',0)
    else
        legend(ud.ax(1),ud.name,0)
    end
end

set(ud.hFig,'Userdata',ud)

%Analize(ud.hFig)
return
%---------------------------------------------------
function AxesOptions

hFig = gcbf;
Button = get(hFig,'SelectionType');
ud = get(hFig,'UserData');



return
%---------------------------------------------------
function LineMarkers
hFig = gcbf;
ud = get(hFig,'UserData');
stat = get(ud.hLineMark,'Checked');
if strcmpi(stat,'On')
    set([ud.horig,ud.hfilt],'Marker','none')
    set(ud.hLineMark,'Checked','Off');
else
    set([ud.horig,ud.hfilt],'Marker','.')
    set(ud.hLineMark,'Checked','On');
end
return
%---------------------------------------------------
function Calculator
hFig = gcbf;
ud = get(hFig,'UserData');
z = ud.x;
t = ud.t;
s{1} = 'Enter expression below:';
s{2} = 'Examples: ';
s{3} = '     z = 0.34*z + 0.01';
s{4} = '     z = z.*(1-sech(1.5*t).^4)';
s{5} = '     t = t(t<232)';
s{6} = '     t(t<datenum(199,5,1))=[]';
s{7} = 'where t - Independent Variable and z - Dependent Variable';
f = char(s);
options.Resize = 'On';
txt = inputdlg(f,'Calculator',1,{'z = z'},options);
if isempty(txt)
    return
end
%ramp = tanh(4.5*t/Tramp);
%ramp = 1 - 1./cosh(3*t/Tramp).^4;
try
    eval([txt{1},';'])
catch
   warndlg('Invalid expression','filter1d') 
   return
end
%check for clipped data
if(length(z)~=length(t))
   z = interp1(ud.t,ud.x,t);
end
%Sort data
[t,ind] = sort(t);
z = z(ind);
ud.x = z;
ud.t = t;
set(ud.horig,'ydata',ud.x,'xdata',ud.t)
set(ud.hRestore,'Enable','on')
set(hFig,'Userdata',ud)
Analize
return
%---------------------------------------------------
function ResampleInterpolate
hFig = gcbf;
ud = get(hFig,'UserData');
options.Resize = 'On';
if (strcmpi(ud.mode,'OTHER') || strcmpi(ud.mode,'XYS')) 
    if isempty(ud.dtf)
        return
    else
        dtf = ud.dtf;
    end
    dto = ud.dt;
    s{1} = ['Original interval is ',num2str(dto)];
    s{2} = ['Filtered interval is currently ',num2str(dtf)];
    s{3} = 'Enter new interval for filtered data';
    f = char(s);    
    txt = inputdlg(f,'Resample/Interpolate',1,{num2str(dto)},options);
    if isempty(txt)
        return
    end
    eval(['ud.dtf = ',txt{1},';'])
else
    if isempty(ud.dtf)
        return
    else
        dtf = ud.dtf*24*60;
    end
    dto = ud.dt*24*60;
    s{1} = ['Original time-series interval is ',num2str(dto),' minutes'];
    s{2} = ['Filtered time-series interval is currently ',num2str(dtf),' minutes'];
    s{3} = 'Enter new time-interval for filtered data';
    f = char(s);
    txt = inputdlg(f,'Time-interval',1,{num2str(dto)},options);
    if isempty(txt)
        return
    end
    eval(['ud.dtf = ',txt{1},'/24/60;'])
end
tf = ud.tf;
ud.tf = (ud.t(1):ud.dtf:ud.t(end))';
ud.xf = interp1(tf,ud.xf,ud.tf,'cubic');
set(ud.hfilt,'xdata',ud.tf,'ydata',ud.xf)
set(hFig,'Userdata',ud)
return
%-------------------------------------------------
function FindOutliers
%global ud
hFig = gcbf;
ud = get(hFig,'UserData');
dt2 = ud.dt^2; %days
d2 = 0*ud.x;
d2(2:end-1) = diff(ud.x,2,1)/dt2; %second derivative
%defmaxd2 = num2str(rounddec(4*std(d2),3));
defmaxd2 = num2str(4*std(d2));
options.Resize = 'On';
if isempty(ud.xf) % || defmaxd2<=1e-3
    txt = inputdlg('Maximum absolute value of second derivative',...
        'Threshold',1,{defmaxd2},options);
    if isempty(txt)
        return
    end
    d2max = str2double(txt{1});
    %Find outliers
    ind = find(abs(d2)>d2max);
else %Filter has been applied
    if ud.dt~=ud.dtf
        xf = interp1(ud.tf,ud.xf,ud.t,'linear','extrap');
    else
        xf = ud.xf;
    end
    r = ud.x-xf; %residuals
%    defmaxr = num2str(rounddec(3*std(r),3));
    defmaxr = num2str(3*std(r));
    s = cell(2,1);
    s{1} = 'Maximum absolute value for residual';
    s{2} = 'Maximum absolute value for second derivative';
    txt = inputdlg(s,'Thresholds',1,{defmaxr,defmaxd2},options);
    if isempty(txt)
        return
    end
    rmax = str2double(txt{1});
    d2max = str2double(txt{2});    
    %Find outliers
    ind = find(abs(r)>rmax | d2>d2max);
end

set(ud.hout,'Xdata',ud.t(ind),'Ydata',ud.x(ind))
set(ud.hOutliers,'Enable','On')
set(hFig,'Userdata',ud)
return
%---------------------------------------------------------
function ClearOutMark
hFig = gcbf;
ud = get(hFig,'UserData');
set(ud.hout,'Xdata',NaN,'Ydata',NaN)
set(ud.hOutliers,'Enable','off')
set(hFig,'Userdata',ud)
return
%---------------------------------------------------
function Restore
hFig = gcbf;
ud = get(hFig,'UserData');
ud.t = ud.to;
ud.x = ud.xo;
set(ud.horig,'xdata',ud.t,'ydata',ud.x)
set([ud.hfilt,ud.hRes,ud.hCut,ud.hSpecf],'xdata',NaN,'ydata',NaN)
set(ud.hRestore,'Enable','Off')
set(hFig,'Userdata',ud)
return
%-------------------------------------------------
function DeletePoints
hFig = gcbf;
ud = get(hFig,'UserData');
j = ud.t;
k = ud.x;
ind = true(size(j));
while 1
    [jp,kp,button] = ginput(1);
    if button==1
%        d = sqrt((jp-j).^2);
        d = sqrt((jp-j).^2 + (kp-k).^2);
        [m,c] = min(d);
        set(ud.hHi,'Xdata',j(c),'Ydata',k(c))
        [jp,kp,button] = ginput(1);
        if button==1
            ind(c) = false(1);
            %Interpolate point
            k(c) = interp1(j(ind),k(ind),j(c),'linear','extrap');
            set(ud.horig,'Xdata',j,'Ydata',k)
            ind = true(size(j));
        end %if
    else
        break
    end %if
    set(ud.hHi,'Xdata',NaN,'Ydata',NaN)
end
set(ud.hHi,'Xdata',NaN,'Ydata',NaN)
ud.t = j(ind);
ud.x = k(ind);
set(ud.hRestore,'Enable','On')
set(hFig,'Userdata',ud)
return
%------------------------------------------------
function ud = initial_timeseries(ud)

n = 1000;
if 0
    t = linspace(0,4,n)';
    dt = diff(t(1:2));
    ud.dt = dt;
    s1 = 10*sin(2*pi*t*1);
    s2 = 5*sin(2*pi*t*15 + 0.2);
    s3 = 3*sin(2*pi*t*30);
    ud.x = s1 + s2 + s3 + 1*randn(n,1);
else
    dt = 1/24/4; %15 min
    t = dt*(0:n-1)';
    ud.dt = dt;
    fSsa= 24/4383.0; %Solar Semi Annual, cycles per day
    fM2 = 24/12.42; %Main lunar semidiurnal component, cycles per day
    fO1 = 24/25.8203; %Main lunar diurnal component, cycles per day
    fM4 = 24/6.21;  %Compound, cycles per day
    s1 = 0.5*sin(2*pi*t*fSsa);
    s2 = 0.5*sin(2*pi*t*fM2 + 0.1);
    s3 = 1*sin(2*pi*t*fO1 + 0.5);
    s4 = 0.3*sin(2*pi*t*fM4 + 0.2);
    ns = 0.1*sin(2*pi*t*100 + 0.1) + 0.1*randn(n,1);
    ud.x = s1 + s2 + s3 + s4 + ns;
end
ud.t = t;
return

%-----------------------------------------------------
function Setfo
hFig = gcbf;
ud = get(hFig,'UserData');
s = get(ud.hfo,'string');
ud.fo = str2num(s);
set(hFig,'Userdata',ud)
Analize
return
%-----------------------------------------------------
function Setfc
hFig = gcbf;
ud = get(hFig,'UserData');
s = get(ud.hfc,'string');
ud.fc = str2num(s);
%Check cut-off frequency
fn = 1/(2*ud.dt); %Nyquist frequency
if (ud.fc>=fn) 
    h=warndlg('Cut-Off Frequency is larger than Nyquist frequency. No filtering applied.',...
        'filter1d');
    uiwait(h);
end
%Check Filter length
minlen = fix(1/min(ud.fc,fn)/ud.dt);
if ud.len<minlen && ud.len>0
    %ud.len = minlen;
    h=warndlg('Filter length may be too small.',...
        'filter1d');
    uiwait(h);
elseif ud.len==0
    h=warndlg('Filter length equal to 0. No filtering applied',...
        'filter1d');
    uiwait(h);
end
maxlen = 100*minlen;
m = length(ud.x) + ud.len;
if ud.len>maxlen && m>10000
    h=warndlg('Filter length is very large and may cause program to be slow.',...
        'filter1d');
    uiwait(h);
end
% maxlen = 20*minlen;
% if ud.len>maxlen
%     h=warndlg('Filter length is too large.','filter1d');
%     uiwait(h);    
%     ud.len = 3*minlen;
% end
set(ud.hLen,'string',num2str(ud.len));
set(hFig,'Userdata',ud)
drawnow
Analize
return
%-----------------------------------------------------
function SetCol
hFig = gcbf;
ud = get(hFig,'UserData');
s = get(ud.hCol,'string');
ud.col = str2num(s);
set(hFig,'Userdata',ud)
drawnow
Analize
return
%-----------------------------------------------------
function SetLen
hFig = gcbf;
ud = get(hFig,'UserData');
s = get(ud.hLen,'string');
ud.len = fix(str2num(s));
%Check Filter length
minlen = fix(1/ud.fc/ud.dt);
if ud.len<minlen && ud.len>0
%    ud.len = 3*minlen;    
    h=warndlg('Filter length is may be too small.','filter1d');
    uiwait(h);  
elseif ud.len==0
    h=warndlg('Filter length equal 0. No filtering applied.','filter1d');
    uiwait(h);   
end
maxlen = 20*minlen;
if ud.len>maxlen
    h=warndlg('Filter length may be too large. Reducing length.','filter1d');
    uiwait(h);    
    ud.len = 4*minlen;
end
lenmax = length(ud.x)-1;
if ud.len>lenmax
    h=warndlg('Filter length larger than time series. Reducing length.','filter1d');
    uiwait(h);    
    ud.len = lenmax;
end
set(ud.hLen,'string',num2str(ud.len));
set(hFig,'Userdata',ud)
drawnow
Analize
return
%-----------------------------------------------------
function SetType
hFig = gcbf;
ud = get(hFig,'UserData');
s = get(ud.htype,'string');
ind = get(ud.htype,'value');
ud.type = s(ind,:);
if strcmpi(ud.type(1:3),'Ban')
   set(ud.hfo,'Enable','On','BackgroundColor','w')
else
    set(ud.hfo,'Enable','Off','BackgroundColor',0.8*[1,1,1])
end
set(hFig,'Userdata',ud)
return
%-----------------------------------------------------
function SetWindow
hFig = gcbf;
ud = get(hFig,'UserData');
s = get(ud.hwin,'string');
ind = get(ud.hwin,'value');
ud.win = s(ind,:);
set(hFig,'Userdata',ud)
return
%-----------------------------------------------------
function setxlabel
hFig = gcbf;
ud = get(hFig,'UserData');
options.Resize = 'On';
txt = inputdlg('Enter X-axes label','X Label',1,{''},options);
if ~isempty(txt)
    xlabel(ud.ax(1),txt{1},'Fontsize',12)
end
return
%-----------------------------------------------------
function setylabel
hFig = gcbf;
ud = get(hFig,'UserData');
options.Resize = 'On';
txt = inputdlg('Enter Y-axes label','Y Label',1,{''},options);
if ~isempty(txt)
    ylabel(ud.ax(1),txt{1},'Fontsize',12)
end
return
%------------------------------------------------------
function xm = mirrordataends(x,N,nf)
xm = zeros(2*nf+N,1);
xm(1:nf) = -(x(nf+1:-1:2)-x(1))+x(1);
xm(nf+1:N+nf) = x;
xm(N+nf+1:N+2*nf) = -(x(N-1:-1:N-nf)-x(N))+x(N);
return
%-----------------------------------------------------
function Analize(hFig)
if nargin==0
    hFig = gcbf;
end
ud = get(hFig,'UserData');
m = length(ud.x) + ud.len;
iwb = (m>10000);
if iwb
    hwb = waitbar(0.1,'Filtering Data');
end

ud.tf = ud.t; %Important
N = length(ud.x); 
fn = 1/(2*ud.dt); %Nyquist frequency
if (ud.fc<=fn && ud.len>0)  
    if 1 %Data Mirrow Method
        nf = ud.len; %filter length
        x = ud.x;
        N = length(x);
        xm = mirrordataends(x,N,nf);
        [xmf,ud.Rh,ud.Fh] = winsinc(xm,ud.dt,...
            ud.fc,ud.len,ud.win,ud.type,ud.fo);
        ud.xf = xmf(nf+1:N+nf);
    else
        [ud.xf,ud.Rh,ud.Fh] = winsinc(ud.x,ud.dt,ud.fc,...
            ud.len,ud.win,ud.type,ud.fo);
    end
else
    ud.tf = ud.t;
    ud.xf = ud.x;
    df = 2*fn/(N-1); %delta frec
    ud.Fh = [df,fn]; 
    ud.Rh = [1,1]; 
end

%Response function
set(ud.hRes,'xdata',ud.Fh,'ydata',ud.Rh)
set(ud.hResInt,'Enable','On')
%Set limits
ylim = get(ud.ax(2),'ylim');
if ylim(1)<-0.5
    ylim(1) = -0.5;
end
if ylim(2)>=1
    ylim(2) = 1.5;
end
set(ud.ax(2),'ylim',ylim)
set(ud.hCut,'xdata',ud.fc*[1,1],'ydata',ylim)
%legend(ud.ax(3),'Original','Filtered',0)

%Alternate Method for Fixing ends
if 0
    ni = fix(ud.len/4);
    ii = linspace(2*pi,0,ni)';
    coef = sech(ii);
    rcoef = flipdim(coef,1);
    ud.xf(1:ni) = (1-coef).*ud.x(1:ni) + coef.*ud.xf(1:ni);
    ud.xf(end-ni+1:end) = (1-rcoef).*ud.x(end-ni+1:end) + rcoef.*ud.xf(end-ni+1:end);
end

if isempty(ud.dtf)
    ud.dtf = ud.dt;
    ud.tf = ud.t;
end

if ud.dtf~=ud.dt
    if(iwb)
        waitbar(0.5,hwb,'Interpolating');
    end
    ud.tf = (ud.t(1):ud.dtf:ud.t(end))';
    ud.xf = interp1(ud.t,ud.xf,ud.tf,'cubic');
    set(ud.hfilt,'xdata',ud.tf,'ydata',ud.xf)
else
    set(ud.hfilt,'xdata',ud.t,'ydata',ud.xf)    
end

if iwb
    waitbar(0.8,hwb,'Calculating Spectra')
end
%Input Spectral
[f,a] = ampspectra(ud.x,ud.dt);
set(ud.hSpec,'xdata',f,'ydata',a)
set(ud.horig,'xdata',ud.t,'ydata',ud.x)
%Output Spectral
[f,a] = ampspectra(ud.xf,ud.dtf);
set(ud.hSpecf,'xdata',f,'ydata',a)

%xlim = get(ud.ax(2),'xlim');
set(ud.ax(2:3),'xlim',[ud.Fh(1),ud.Fh(end)])
if isempty(ud.name)
    legend(ud.ax(1),'Original','Filtered',0)
else
    legend(ud.ax(1),ud.name,[ud.name,' filtered'],0)
end
% if strmatch(ud.mode,{'Excel','TXT','COOPS'})
%     datetick(ud.ax(1),'x',2)
% else %*.xys or *.*
%     set(ud.ax(1),'XLimMode','auto',...
%         'XTickLabelMode','auto','XTickMode','auto')    
% end
if iwb
    waitbar(1,hwb,'Finsished')
    pause(0.1)
    close(hwb)
end
set(ud.ax(1),'XMinorTick','On')
set(hFig,'Userdata',ud)
return
%-------------------------------------------------------------------
function saveTAPxlsfile(ud,fname)
%Saves data to a Microsoft Excel Workbook in TAP format
%Note: Currently only works with water level data
sd = datestr(ud.t(:),'mm/dd/yy HH:MM');
% data = cell(n+1,3);
data(1,:) = {'Record','Date and Time','Water level'};
for k=1:n
    data{k+1,1} = k;
    data{k+1,2} = sd(k,:);
    data{k+1,3} = ud.xf(k);
end
[type, sheets] = xlsfinfo(fname);
ns = length(sheets);
f{1} = 'Select output sheet name or enter new name';
f{2} = ' ';
f(3:ns+2) = sheets;
options.Resize = 'On';
txt = inputdlg(char(f),'Output Sheet',1,{'Filt'},options);
if isempty(txt)
    return
end
%    warning off
status = xlswrite(fname,data,txt{1},'A1');
%    warning on
if (status==0)
    errordlg('Could not save file. Make sure file is not open.','File Error');
    return
end
return
%-------------------------------------------------------------------
function savexysfile(ud,fname)
%saves an xys file
tf = ud.tf(:);
if strmatch(ud.mode,{'Excel','ASCII'}) %excel or prediction file
    tf = (tf-tf(1))*24;
end
tf = round(tf*60)/60; %round to the nearest minute
data = [tf, ud.xf(:)];
n = length(ud.tf);
fid = fopen(fname,'w');
fprintf(fid,'XYS 1 %d "Curve"\n',n);
fprintf(fid,'%f %f\n',data(1:n-1,:)');
fprintf(fid,'%f %f',data(n,:)');
fclose(fid);
return
%-------------------------------------------------------------------
function savetsdfile(ud,fname)
%Saves a Time Series Data file 
%
if (strcmpi(ud.mode,'OTHER') || strcmpi(ud.mode,'XYS')) %hours
    t = (ud.tf(:)-ud.tf(1))*3600; %convert from hours to seconds
    s = cell(4,1);
    s{1} = '';
    s{2} = 'Input Starting Time:';
    s{3} = '';
    s{4} = 'Format: mm/dd/yyyy HH:MM:SS';
    f = char(s);
    options.Resize = 'On';
    txt = inputdlg(f,'Input',1,{'01/01/2000 00:00:00'},options);  
    str = txt{1};
else % strmatch(ud.mode,{'TAPExcel','TAPpred','COOPS','TSD','ASCII'}) %days   
    t0 = ud.tf(1);
    t = (ud.tf(:)-t0)*86400; %convert to seconds and remove start    
    str = datestr(t0,'mm/dd/yyyy HH:MM:SS');
end
t = rounddec(t,3); %round to avoid precision errors
data = [t,ud.xf(:)];
nrow = length(ud.xf);
if isempty(ud.name)
%     name = 'Sample'; 
    name = 'Water Level'; 
else
    name = ud.name;
end
if isempty(ud.curvetype)
    curvetype = 'Unassigned';
else
    curvetype = ud.curvetype;
end
fid = fopen(fname,'w');
fprintf(fid,'TIME_SERIES\n');
fprintf(fid,'"%s"      "%s"      2      %d     "%s"\n',name,curvetype,nrow,str);
fprintf(fid,' %16.2f %10.3f\n',data');
fclose(fid);
return
%-------------------------------------------------------------------
function savedata
hFig = gcbf;
ud = get(hFig,'UserData');
nofilt = 0;
if isempty(ud.xf)
%     errordlg('Click on Analyize button before saving output file.')
%     return
    nofilt = 1;
    ud.xf = ud.x;
    ud.tf = ud.t;
    ud.dtf = ud.dt;    
end
if ~isempty(ud.pathout)
    cd(ud.pathout)
elseif ~isempty(ud.pathin)
    cd(ud.pathin)
end
% if strcmpi(ext,'.xls') %to avoid overwrite message
filterspec = {'*.tsd','SMS 10.1+ Time Series Data File (*.tsd)';...
            '*.xys','SMS 10.0+ XY Series File(*.xys)';...
            '*.xls','TAP MS Excel Workbook (*.xls)';...
            '*.*','Generic ASCII File (*.*)';...
            };   
[oname,opath,filterindex] = uiputfile(filterspec,'Save Output File:',ud.fout);
if isnumeric(oname) || isempty(oname)
    return
end
ud.pathout = opath;
cd(ud.pathfilter1d)
[ipath,aname,ext] = fileparts(oname);
%Make sure extension matches file type
if strcmpi(ext,'.tsd')
    filterindex = 1;
elseif strcmpi(ext,'.xys')
    filterindex = 2;
elseif strcmpi(ext,'.xls')
    filterindex = 3;
else
    filterindex = 4;
end
fname = [opath,oname];
ud.fout = oname;
% n = length(ud.xf);
switch (filterindex)
    case 1 %TSD file
        savetsdfile(ud,fname)
    case 2 %XYS file
        savexysfile(ud,fname)
    case 3 %TAP xls file
        saveTAPxlsfile(ud,fname)
    case 4 %Standard ASCII
        ud = saveASCIIfile(ud,fname);
end
%Save output file name
fout = 'filter1d_outputfiles.txt';
fid = fopen(fout,'a');
fprintf(fid,'%s\n',fname);
fclose(fid);

if (nofilt)
    ud.xf = [];
    ud.tf = [];
    ud.dtf = [];
end

set(ud.hfout,'string',oname)
set(hFig,'Userdata',ud)
return
%------------------------------------------------------------------
function ud = loadTAPxlsfile(fname,ud)
%Reads a Microsoft Excel File in TAP format
%
[type, sheets] = xlsfinfo(fname);
nsh = length(sheets);
s = cell(nsh,1);
s{1} = 'Select Worksheet: (e.g. 3)';
for k=1:nsh
    s{k+1} = [num2str(k),'. ',sheets{k}];
end
f = char(s);
options.Resize = 'On';
txt = inputdlg(f,'Input Sheet',1,{'1'},options);
if isempty(txt)
    set(ud.hfin,'string',ud.fin,'ForegroundColor','k')
    return
end
sh = sheets{str2double(txt{1})};
[data,txt,raw] = xlsread(fname,sh);
N = size(data,1);
%Check for header
if isnan(raw{1,1}) || ischar(raw{1,1})
    istart = 2;
else
    istart = 1;
end
if length(raw{5,2})==9 %size(data,2)>=4 && ~isempty(data(1,4))
    x = data(:,4);
    t = datenum(raw(istart+(0:N-1),2)) + data(:,3);
else
    x = data(:,3);
    t = datenum(raw(istart+(0:N-1),2));
end
ud.t = t(:);
ud.x = x(:);
%    set([ud.hStrHead,ud.hHead,ud.hStrCol,ud.hCol],'Visible','Off')
ud.dt = (t(2) - t(1)); %days
%    set(ud.hdt,'string',ud.dt);
%    set(ud.hStrdt,'String','dt, min')
xlabel(ud.ax(2),'Frequency, cycles/day')
ud.mode = 'TAPExcel';
return
%-------------------------------------------
function ud = loadTAPpredfile(fname,ud)
%Reads a TAP tidal prediction file
%
fid = fopen(fname,'r');
%read header lines
c = 1; %counter
while 1
    row = fgetl(fid);
    if ~ischar(row)
        %bad file format
        set(ud.hfin,'string',ud.fin,'ForegroundColor','k')
        return
    end
    nr = length(row);
    if strcmp(row(1:min(4,nr)),'Date')
        fclose(fid);
        break
    end
    c = c + 1;
end
ud.headerlines = c;
try
    [yyyy,mm,dd, HH,MM,SS, ud.x] = textread(fname,...
        '%d-%d-%d %d:%d:%d %f','headerlines',ud.headerlines);
catch
    errordlg('Could not read input file. Check input file and try again')
end
ud.t = datenum(yyyy,mm,dd,HH,MM,SS);
ud.dt = ud.t(2) - ud.t(1); %days
ud.mode = 'TAPpred';
xlabel(ud.ax(2),'Frequency, cycles/day')
return
%----------------------------------------------------------
function ud = loadCOOPSfile(fname,ud)
%Loads a NOAA COOPS file
%
ud.headerlines = 7;
try
    [yyyy,mm,dd,HH,MM,ud.x] = textread(fname,...
        '%*s %4d%2d%2d %d:%d %*f %f','headerlines',ud.headerlines);
catch
    errordlg('Could not read input file. Check input file and try again')
end
ud.t = datenum(yyyy,mm,dd,HH,MM,0);
ud.dt = ud.t(2) - ud.t(1); %days
ud.mode = 'COOPS';
xlabel(ud.ax(2),'Frequency, cycles/day')
return
%----------------------------------------------------------------
function ud = loadASCIIfile(fname,ud)
% Loads a Generic ASCII File
% assuming that all the lines haves the same format 
% written by Alex Sanchez, Last modified May 31, 2012
fid = fopen(fname,'r');
f = cell(1,14);
f{1} = 'First 10 lines of file:             ';
f{2} = ' ';
for k=1:10
    f{k+2} = fgetl(fid);
end
f{13} = ' ';
f{14} = 'Number of header lines';
ff = char(f);
ff(double(ff)==9) = char(32); %Replace tabs with spaces
fclose(fid);
options.Resize = 'On';
txt = inputdlg({ff;'Format (e.g. ''yyyymmdd HH:MM * z'')';...
    'Time inteval or Spacing (leave empty to determine from data)'},...
    'Set Parameters',1,{num2str(ud.headerlines);...
    num2str(ud.fmt); ' '},options);
if isempty(txt)
    if isfield(ud,'fin')
        set(ud.hfin,'string',ud.fin,'ForegroundColor','k')
    end
    return
end
%Allows the use of expression in the input dialog
eval(['ud.headerlines = ',txt{1},';'])
eval(['fmt = ''',txt{2},''';'])
ud.fmt = fmt;
%Initialize variables
year = 0;
month = 0;
day = 0;
hr = 0;
min = 0;
sec = 0;
outvar = ['[',fmt,']'];

outvar = regexprep(outvar,'yyyy','year ');
outvar = regexprep(outvar,'yy','year ');
ind = findstr(outvar,'mm');
if isempty(ind)
    outvar = regexprep(outvar,'m','month ');
else
    outvar = regexprep(outvar,'mm','month ');    
end
outvar = regexprep(outvar,'HH','hr ');
outvar = regexprep(outvar,'H','hr ');
outvar = regexprep(outvar,'MM','min ');
outvar = regexprep(outvar,'M','min ');
outvar = regexprep(outvar,'SS','sec ');
outvar = regexprep(outvar,'S','sec ');
outvar = regexprep(outvar,'*','');
outvar = regexprep(outvar,'/','');
outvar = regexprep(outvar,':','');
outvar = regexprep(outvar,'-','');

ind = findstr(fmt,'x');
if ~isempty(ind)
    ud.mode = 'OTHER';
else
    ud.mode = 'ASCII';
end

ind = findstr(fmt,'dd');
if isempty(ind)
    fmt = regexprep(fmt,'d','%2d');
    outvar = regexprep(outvar,'d','day ');
else
    fmt = regexprep(fmt,'dd','%2d');
    outvar = regexprep(outvar,'dd','day ');
end

fmt = regexprep(fmt,'yyyy','%4d');
fmt = regexprep(fmt,'yy','%2d');
fmt = regexprep(fmt,'mm','%2d');
fmt = regexprep(fmt,'m','%2d');
fmt = regexprep(fmt,'HH','%2d');
fmt = regexprep(fmt,'H','%2d');
fmt = regexprep(fmt,'MM','%2d');
fmt = regexprep(fmt,'M','%2d');
fmt = regexprep(fmt,'SS','%2d');
fmt = regexprep(fmt,'S','%2d');
fmt = regexprep(fmt,'x','%f');
fmt = regexprep(fmt,'z','%f');
fmt = regexprep(fmt,'t','%f');
fmt = regexprep(fmt,'*','%*s');

%Read first line to get terminator
fid = fopen(fname,'r');
aline = fgets(fid);
fclose(fid);
if ~isempty(regexp(aline,'\n','match'))
    fmt = [fmt,'%*[^\n]'];
elseif ~isempty(regexp(aline,'\r','match'))
    fmt = [fmt,'%*[^\r]'];
end
ud.format = fmt;
% outfmt = ['[',ud.fmt,']'];
% outfmt = regexprep(outfmt,'/',' ');
% outfmt = regexprep(outfmt,':',' ');

str = [outvar,' = textread(fname,fmt,',...
    '''headerlines'',ud.headerlines,''bufsize'',1e5);'];

try
    eval(str)
    if strcmp(ud.mode,'OTHER')
        ud.t = x;
        ud.x = z;
    else
        if ~exist('t','var')
            t = datenum(year,month,day,hr,min,sec);
        end
        ud.t = t;
        ud.x = z;
    end
catch
    errordlg('Could not read input file. Check input file and try again')
end
if ~isempty(deblank(txt{3}))
    eval(['ud.dt = ',txt{3},';'])
else
    ud.dt = ud.t(2) - ud.t(1);
end

% try
%     %          [ud.t,ud.x]= textread(fname,fmt,'headerlines',ud.headerlines,'bufsize',1e5);
%     fid = fopen(fname,'r');
%     for k=1:ud.headerlines
%         fgetl(fid);
%     end
%     data = fscanf(fid,fmt,[2,inf])';
%     ud.t = data(:,1);
%     ud.x = data(:,2);
%     %        str = [outfmt,'= textread(fname,fmt,''headerlines'',ud.headerlines);'];
%     %        eval(str)
% catch
%     errordlg('Could not read input file. Check input file and try again')
% end

%    set([ud.hStrHead,ud.hHead,ud.hStrCol,ud.hCol],'Visible','On')
return
%------------------------------------------------------------------------
function ud = saveASCIIfile(ud,fname)
%write a Generic ASCII file
%
options.Resize = 'On';
txt = inputdlg({'Text for Header (optional)',...
    'Format (e.g. ''yyyymmdd HH:MM z'')'},...
    'Set Parameters',1,{'';num2str(ud.fmt)},options);
if isempty(txt)
    return
end
t = ud.tf(:);
z = ud.xf(:);
[year,mon,day,hr,min,sec] = datevec(t);

header = deblank(txt{1});
fmt = txt{2};
fmtdat = fmt;

fmt = regexprep(fmt,'dd','%2d');
fmt = regexprep(fmt,'yyyy','%4d');
fmt = regexprep(fmt,'yy','%02d');
fmt = regexprep(fmt,'mm','%02d');
fmt = regexprep(fmt,'HH','%02d');
fmt = regexprep(fmt,'MM','%02d');
fmt = regexprep(fmt,'SS','%02d');
fmt = regexprep(fmt,'x','%f');
fmt = regexprep(fmt,'z','%f');
fmt = regexprep(fmt,'t','%f');
fmt = [fmt,'\n'];
 
fmtdat = regexprep(fmtdat,'dd','day,');
fmtdat = regexprep(fmtdat,'yyyy','year,');
fmtdat = regexprep(fmtdat,'yy','year,');
fmtdat = regexprep(fmtdat,'mm','mon,');
fmtdat = regexprep(fmtdat,'HH','hr,');
fmtdat = regexprep(fmtdat,'MM','min,');
fmtdat = regexprep(fmtdat,'SS','sec ');
fmtdat = regexprep(fmtdat,'x','t');
fmtdat = regexprep(fmtdat,':','');
fmtdat = regexprep(fmtdat,'\','');
fmtdat = regexprep(fmtdat,'/','');

eval(['data = [',fmtdat,'];'])

fid = fopen(fname,'w');
if ~isempty(header)
    fprintf(fid,'%s\n',header);
end
fprintf(fid,fmt,data');
fclose(fid);
ud.format = fmtdat;
ud.fmt = fmt;
return
%-------------------------------------------------------------------
function ud = loadTSDfile(fname,ud)
%Reads an SMS time series file
%
fid = fopen(fname,'r');
fgets(fid); %skip first line
C = textscan(fid,'%q %q %f %f %q',1);
ud.name = C{1}{1};
ud.curvetype = C{2}{1};
num = [C{3}, C{4}];
t0 = datenum(C{5},'mm/dd/yyyy HH:MM:SS');
curvetypes = {'Unassigned',...
            'Velocity - Components',...
            'Velocity - Mag. & Dir.',...
            'XY Series'};
fmt = '%f';
for k=2:num(1)
    fmt = [fmt,' %f'];
end
data = fscanf(fid,fmt,num)';
fclose(fid);
switch ud.curvetype
    case {curvetypes{1}}
        if num(1)>2
            f{1} = ['File has ',num2str(num(1)),' columns'];
            f{2} = 'Enter column number for input data';
            s = char(f);
            ans = inputdlg(s,'Input Column?');
            nth = str2num(ans{1});
        else
            nth = num(1);
        end        
        if nth>num(1)
            nth = num(1);
        end
        ud.x = data(:,nth);
    case {curvetypes{2},curvetypes{3}}
        ButtonName = questdlg('Select input data for filtering',...
            'Current Data','Magnitude',...
            'X-component','Y-Component','Magnitude');
        switch ud.curvetype
            case curvetypes{2}
                u = data(:,2);
                v = data(:,2);
                m = sqrt(u.^2 + v.^2);
            case curvetypes{3}
                m = data(:,2);
                u = m.*sin(data(:,3)*pi/180);
                v = m.*cos(data(:,3)*pi/180);
        end %switch
        switch ButtonName
            case ('Magnitude')
                ud.x = m;
            case ('X-component')
                ud.x = u;
            case ('Y-Component')
                ud.x = v;
        end %switch
        
end %switch
ud.t = t0  + data(:,1)/86400; %note conversion to days
ud.dt = ud.t(2) - ud.t(1);
ud.mode = 'TSD';
return
%------------------------------------------------------------------
function ud = getdata(fname,ud)
%Loads an input file specified as input or from a dialog
if nargin==0
    hFig = gcbf;
    ud = get(hFig,'UserData');
    if ~isempty(ud.pathin)
        cd(ud.pathin) 
    end
    filterspec = {'*.tsd','SMS 10.1+ Time Series Data file (*.tsd)';...
        '*.xys','SMS 10.0+ XY Series File (*.xys)';...
        '*.xls','TAP Excel Workbook (*.xls)';...
        '*.txt','TAP Prediction File (*.txt)';...
        '*.dat','COOPS ASCII format (*.dat)';...
        '*.*','Generic ASCII File (*.*)';...
        };
    [iname,ipath,filterindex] = uigetfile(filterspec,'Locate Input-File:');
    %'*.txt, NOAA COOPS format';
    if isnumeric(iname) || isempty(iname)
        return
    end
    ud.pathin = ipath;
    cd(ud.pathfilter1d)
    %set(ud.hfin,'string','Loading ...','ForegroundColor','r')
    %Make sure name has an appropriate extension
    [apath,aname,ext] = fileparts(iname);
    if isempty(ext)
        switch filterspec
            case 1
                ext = '.tsd';
            case 2
                ext = '.xys';
            case 3
                ext = '.xls';
            case 4
                ext = '.txt';
            case 5
                ext = '.dat';
            case 6
                ext = '.txt';
        end
        iname = [iname,ext];
    end
    fname = [ipath,iname];
else
    [ipath,iname,ext] = fileparts(fname);
    iname = [iname,ext];
end

if strcmpi(ext,'.xls') && (filterindex~=3)
    h = errordlg('To open an Excel choose the correct file type');
    uiwait(h)
    getdata
    return
end
if strcmpi(ext,'.tsd')
    filterindex = 1;
end
if strcmpi(ext,'.xys')
    filterindex = 2;
end
if filterindex>1
   ud.name = '';
   ud.curvetype = '';
end
switch filterindex
    case (1) %SMS Time Series Data File
        ud = loadTSDfile(fname,ud);
    case (2) % SMS 10 XY Series
        [ud.t,ud.x] = textread(fname,'%f %f','headerlines',1);
        ud.dt = ud.t(2) - ud.t(1);
        ud.mode = 'XYS';
%    xlabel(ud.ax(2),'Freq., cycles/hr')
    case (3) %TAP Excel File (*.xls)
        ud = loadTAPxlsfile(fname,ud);  
    case (4) % TAP Prediction File (*.txt)
        ud = loadTAPpredfile(fname,ud);
    case (5) % COOPS format (*.txt)
        ud = loadCOOPSfile(fname,ud);    
    case (6) %*.* Generic ASCII file
        ud = loadASCIIfile(fname,ud);            
end %select
ud.dtf = [];
ud.fin = iname;
ud.pathin = ipath;
ud.xo = ud.x;
ud.to = ud.t;
ud.fc = 1./(ud.dt*10); %initial value
ud.fc = rounddec(ud.fc,4);
minlen = 4*fix(1/ud.fc/ud.dt);
ud.len = 3*minlen;
N = length(ud.x);
if(ud.len>N)
    ud.len = N-1; 
end

%Detect NaNs
isfine = isfinite(ud.x);
ud.x = ud.x(isfine);
ud.t = ud.t(isfine);

%Sort data
[temp,ind] = sort(ud.t);
if ~isequal(temp,ud.t)
    h = warndlg('Times steps found out of order. Data will be sorted.','filter1d');
    uiwait(h);
    ud.t = temp;
    ud.x = ud.x(ind);
end

%Look at time steps
d = diff(ud.t); %time intervals in minutes
maxdt = max(d);
ivardt = (d>1.00001*ud.dt | d<0.999999*ud.dt); 

%Eliminate repeated times steps
ind = find(d==0); %repeated times
if ~isempty(ind)
    h = warndlg('Repeated time steps found. Repeated times will be deleted.','filter1d');
    uiwait(h);
    ud.t(ind+1) = []; %delete repeated points,
    ud.x(ind+1) = []; %delete repeated points
end

%Check time-series length
nc = round((ud.t(end)-ud.t(1))/ud.dt) + 1;
N = length(ud.t); %might have changed
ivarn = (N~=nc);

%Inerpolate to regulare time step
if  any(ivardt) || ivarn
%     s{1} = 'Variable time interval or missing data identified.';
%     s{2} = 'Data will be linearly interpolated to regular time interval';
%     f = char(s);
%     h = warndlg(f,'filter1d');
    if strmatch(ud.mode,{'TAPExcel','TAPpred','COOPS','TSD','ASCII'})
        s{1} = '                                 TIME INTERVAL';
        s{2} = 'Variable time intervals or missing data identified.';
        s{3} = 'Regular time intervals are required for filtering.';
        s{4} = ['The calculated time step is ', num2str(ud.dt*1440),' minutes.'];
        s{5} = 'Data will be interpolated to this time step.';
        s{6} = '';
        s{7} = '                                  DATA GAPS';
        s{8} = 'Specify the maximum time interval in minutes for interpolation.';
        s{9} = 'The time series will be cropped to the first gap exceeding this value.';
        s{10} = 'Clicking on Cancel will load entire time series without cropping.';
        s{11} = '';
        f = char(s);
        options.Resize = 'On';
        txt = inputdlg(f,'Input',1,{num2str(maxdt*1440)},options);
        if ~isempty(txt)
            maxdt = str2double(txt{1}/1440); %convert from minutes to days      
        else
            maxdt = ud.t(end)-ud.t(1);
        end
    else
        s{1} = '                    INTERVAL';
        s{2} = 'Variable intervals or missing data identified. ';
        s{3} = 'Regular intervals are required for filtering.';
        s{4} = '';
        s{5} = '                   DATA GAPS';
        s{6} = 'Specify the maximum interval for interpolation.';
        s{7} = '';
        s{8} = 'The data will be cropped to the first gap exceeding this value.';
        s{9} = 'Clicking on Cancel will load entire time series without cropping.';
        s{10} = '';
        f = char(s);
        options.Resize = 'On';
        txt = inputdlg(f,'Data Gaps',1,{num2str(maxdt)},options);
        if ~isempty(txt)
            maxdt = str2double(txt{1}/1440); %convert from minutes to days      
        else
            maxdt = ud.t(end)-ud.t(1);
        end
    end        
    igap = find(d>maxdt,1,'first');
    if isempty(igap)
        igap = N;
    end
    ud.maxdt = maxdt;
    nc = round((ud.t(igap)-ud.t(1))/ud.dt) + 1;
    t = ud.t(1) + ud.dt*(0:nc-1)';
    ud.x = interp1(ud.t,ud.x,t,'linear','extrap');
    ud.t = t; %replace old time with regular vector
end
ud.xf = [];
ud.tf = [];

%Input Amplitude Spectra
[f,a] = ampspectra(ud.x,ud.dt);
set(ud.hSpec,'Xdata',f,'Ydata',a)
fn = 1/(2*ud.dt); %Nyquist frequency
df = 2*fn/(N-1); %delta frec
set(ud.ax(2:3),'xlim',[df,fn])
set(ud.hCut,'xdata',NaN,'ydata',NaN)
set(ud.hfc,'String',num2str(ud.fc))
set(ud.hLen,'String',num2str(ud.len))
set([ud.hRes,ud.hSpecf,ud.hfilt],'xdata',NaN,'ydata',NaN)
set(ud.horig,'xdata',ud.t,'ydata',ud.x)
set(ud.ax(1),'XLimMode','auto','XTickLabelMode','auto','XTickMode','auto',...
        'YLimMode','auto','YTickLabelMode','auto','YTickMode','auto') 
set([ud.hRestore,ud.hResInt] ,'Enable','off')

if strmatch(ud.mode,{'TAPExcel','TAPpred','COOPS','TSD','ASCII'})
    datetick(ud.ax(1),'x',2)
    xlabel(ud.ax(1),'')
% else %*.xys or *.*
%     set(ud.ax(1),'XLimMode','auto',...
%         'XTickLabelMode','auto','XTickMode','auto')    
end

if any(isfinite(ud.x))
    if isempty(ud.name)
        legend(ud.ax(1),'Original',0)
    else
        legend(ud.ax(1),ud.name,0)
    end
end

%Reset Output File
set(ud.hfout,'string','')
set(ud.hfin,'string',iname,'ForegroundColor','k')
if nargin==0
    set(hFig,'Userdata',ud)
end
return
%---------------------------------------------------------------------
function [xf,Rh,Fh] = winsinc(x,dt,fc,ns,win,ftype,fo)
% Applies a windowed sinc filter
%
% USAGE:
%   [xf,Rh,Fh] = winsinc(x,dt,fc,ns,win,ftype,fo)
%
%INPUT VARIABLES:
%   x = vector
%   dt = Sampling Interval
%   fc = Cutoff frequency
%   ns = Width of filter 
%       (>>ns produces a steeper cutoff
%       but is more computationally expensive)
%       2*ns + 1 = Order of Filter, (Number of Filter Coef.)
%   win = Window Type:
%       'welch','parzen','hanning','hamming'
%       'blackman','lanczos','kaiser','gauss'
%   
%   ftype = 'low', 'high' o 'band'
%   fo = Central Frequency (for bandpass only) 
%
%OUTPUT VARIABLES:
%   xf = filtered vector
%   Rh = Filter Response
%   Fh = Frequencies for Filter Response
%
%Alejandro Sanchez

N = length(x); 
fn = 1/(2*dt); %Nyquist frequency
df = 2*fn/(N-1); %delta frec
ftype = lower(ftype);
if strcmpi(ftype(1:3),'ban')
    fcc = fo;
    fo = fc;
    fc = fcc;
end
%Frecuencia de corte para pasa alta
if strcmpi(ftype(1:3),'hig')
    fc = fn - fc;   
    if mod(ns,2)>0
       ns = ns + 1;
    end
end
n = -ns:ns;
w = 0*n; %initialize
c = 0*n;
an = abs(n);
ic = (an<ns & n~=0);
switch lower(win(1:3))
    case 'wel' %welch
        w(ic) = 1 - (n(ic)/ns).^2;
        w(~ic) = 0;
    case 'han' %hanning
        alpha2 = 1/2;
        w(ic) = alpha2 + (1 - alpha2)*cos(pi*n(ic)/ns);
        w(~ic) = 0;
    case 'ham' %hamming
        alpha2 = 0.54;
        w(ic) = alpha2 + (1 - alpha2)*cos(pi*n(ic)/ns);
        w(~ic) = 0;
    case 'bla' %blackman
        w(ic) = 0.42 + 0.5*cos(pi*n(ic)/ns) ...
            + 0.08*cos(2*pi*n(ic)/ns);
        w(~ic) = 0;
    case 'lan' %lanczos
        w(ic) = sin(pi*n(ic)/ns)./(pi*n(ic)/ns);
        w(~ic) = 0;
    case 'gau' %gauss
        sigma = ns/2;
        w(ic) = exp(-(n(ic)/sigma).^2);
        w(~ic) = 0;
    case 'box' %boxcar
        w = ones(size(n));
    otherwise
        error('Invalid Window')
        %w = ones(ns,1)/ns;
end %switch
%Basic Filter (sinc)
ho = fc/fn; %0<ho<1
iz = (n==0);
c(~iz) = sin(pi*n(~iz)*ho)./(pi*n(~iz));
%Obtain filter Coefficients
%Multiplication is time domain
%is a convolution in frequency domain
h = w.*c;
h(iz) = ho; 
%Make sure sum is equal to 1
h = h/sum(h);
%High pass, clever trick
if strcmpi(ftype(1:3),'hig')
    h(2:2:length(h)) = -1*h(2:2:length(h));
end
%Filter Frequencies
Fh = df:df:fn; 
Rh = 0*Fh; %initialize
switch ftype(1:3)
case {'low','hig'}
    %Filter Response
    for k=1:length(Fh)
        Rh(k) = sum(h.*cos(2*pi*n*Fh(k)/(2*fn)));
    end
    try 
        xf = conv2(h,1,x,'same'); %note that conv2 is used to use 'same' option       
    catch %in case conv2 is not available
        xf = calcconv(h,x);
    end
case 'ban'
    %Filter Response
    for k=1:length(Fh);
        Rh(k) = sum(h.*cos(2*pi*n*(Fh(k)-fo)/(2*fn)));
    end
    %Corrimiento de los pesos del filtro
    th = n*dt;
    i = sqrt(-1); %Redundant but instructive
    hpri = h.*exp(2.0*i*pi*fo.*th);
    try 
        xfc = conv2(hpri,1,x,'same'); %note that conv2 is used to use 'same' option 
    catch %in case conv2 is not available
        xfc = calcconv(hpri,x);
    end
    xf = 2*real(xfc);         
end
xf = xf(:);
if strcmpi(ftype,'band-stop')
    xf = x - xf;
    Rh = 1 - Rh;
end
return
%---------------------------------------------------------------
function xc = calcconv(w,x)
%Numerical computation of convolution
m = length(w);
n = length(x);
len = m+n-1;
y = zeros(len,1); %Initialize
for i=1:len
    y(i) = 0;
    s1 = max(1,i+1-n);
    s2 = min(m,i);
    for j=s1:s2
       y(i) = y(i) + w(j)*x(i+1-j);
    end
end
if mod(m,2)==1 %odd
    ns = (m-1)/2;
    ind = ns+1:m+n-ns-1; 
else %even  
    ns = m/2;
    ind = ns+1:m+n-ns;
end
xc = y(ind);
%---------------------------------------------------------------
function [f,a] = ampspectra(x,dt)
% Calculates the amplitude spectra
% using the windowed periodogram method
x = x(:);
isf = isfinite(x);
x = detrend(x(isf)); %also removes mean
N = sum(isf);
fn = 1/(2*dt); %Nyquist frequency
nfft = max(256,2^(nextpow2(N)));
f = (1:nfft/2)'/(nfft/2)*fn; %frequencies, Note: starts from 1
w = 0.54 - 0.46*cos(2*pi*(0:N-1)'/(N-1)); %Compute hamming window
x = w.*x; %apply window
R = w'*w;  %compute sum of squares of window
c = fft(x,nfft); %Compute Fourier Coefficients, Pad with zeros to length nfft
K = sqrt(5/2); %coefficient for hamming window
a = K*abs(c(2:nfft/2+1))/R; %amplitudes, note skip of mean
return
%---------------------------------------------------------------
function y = rounddec(x,n)
f = 10.^n;
y = round(x.*f)./f;
return
%----------------------------------------------------------------
function x = ramp(t,x,T)
r = 1-sech(2.5*(t-t(1))./T).^4;
x = x.*r;
return
%-----------------------------------------------------------
function AboutFilter1d
f{1} = 'Program filter1d ';
f{2} = 'Version 1.0';
f{3} = 'Programmed in Matlab 7.9 (R2009b)';
f{4} = 'Last Modified Dec 17, 2010 ';
f{5} = 'Funded by the U.S. Army Corps of Engineers';
f{6} = 'Coastal Inlets Research Program';
f{7} = 'Author: Alejandro Sánchez        ';
msg = char(f);
msgbox(msg,'About filter1d')
return
%-----------------------------------------------------------
function Documentation
% f{1} = '';
% f{2} = '     See filter1d User''s manual at         ';
% f{3} = '';
% f{4} = '  http://cirp.usace.army.mil/wiki/Filter1d  ';
% f{5} = '';
% msg = char(f);
% msgbox(msg,'Documentation')
web http://cirp.usace.army.mil/wiki/Filter1d -browser 
return
%-----------------------------------------------------------
function Disclaimer
f{1} = '';
f{2} = 'PLEASE READ CAREFULLY';
f{3} = '';
f{4} = 'This program is provided ''as is'' without any guaranty or warranty.';
f{5} = '';
f{6} = 'The user assumes all risk involved in the application of the present software.';
f{7} = '';
msg = char(f);
msgbox(msg,'Disclaimer')
return



