From c6bd6abffc270004b2ac647505e1f83887c5abe3 Mon Sep 17 00:00:00 2001
From: Prousselot <pierre.rousselot@ird.fr>
Date: Tue, 24 Dec 2019 16:08:03 +0100
Subject: [PATCH] Cleaned version

---
 .gitignore                            |    4 +
 filter/clean_median.m                 |  201 +++++
 filter/diurne.m                       |  176 ++++
 filter/fft_filter.m                   |   80 ++
 filter/filterfft.m                    |   71 ++
 filter/filtfilt.m                     |  310 +++++++
 filter/hamm.m                         |   14 +
 filter/hammfilter_nodec.m             |   29 +
 filter/hamming.m                      |   36 +
 filter/lisse.m                        |  131 +++
 filter/medfilt1_perso.m               |   22 +
 filter/median_filter.m                |  106 +++
 filter/mfilter.m                      | 1075 +++++++++++++++++++++++++
 filter/naninterp.m                    |    5 +
 filter/nanmedfilt2.m                  |   44 +
 filter/ndnanfilter.m                  |  484 +++++++++++
 filter/spikeRemoval.m                 |  601 ++++++++++++++
 moored_adcp_proc/adcp_check_surface.m |   12 +-
 moored_adcp_proc/adcp_filt_sub.m      |  129 +--
 moored_adcp_proc/adcp_shadowzone.m    |    6 +-
 moored_adcp_proc/adcp_surface_fit.m   |   88 +-
 moored_adcp_proc/filtfilt.m           |    2 +-
 template_get_adcp_data.m              |  368 +++++----
 23 files changed, 3735 insertions(+), 259 deletions(-)
 create mode 100644 .gitignore
 create mode 100644 filter/clean_median.m
 create mode 100644 filter/diurne.m
 create mode 100644 filter/fft_filter.m
 create mode 100644 filter/filterfft.m
 create mode 100644 filter/filtfilt.m
 create mode 100644 filter/hamm.m
 create mode 100644 filter/hammfilter_nodec.m
 create mode 100644 filter/hamming.m
 create mode 100644 filter/lisse.m
 create mode 100644 filter/medfilt1_perso.m
 create mode 100644 filter/median_filter.m
 create mode 100644 filter/mfilter.m
 create mode 100644 filter/naninterp.m
 create mode 100644 filter/nanmedfilt2.m
 create mode 100644 filter/ndnanfilter.m
 create mode 100644 filter/spikeRemoval.m

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..8a8a543
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+*.000
+*.pdf
+*.csv
+*.mat
diff --git a/filter/clean_median.m b/filter/clean_median.m
new file mode 100644
index 0000000..49b5f7a
--- /dev/null
+++ b/filter/clean_median.m
@@ -0,0 +1,201 @@
+function [par,info_rejet_med] = clean_median(param,window_lgth,nb_std,delta_lim,iter,fillval)
+% par = clean_median(param,window_length,nb_std,delta_lim,iter,fillval);
+% param = vector of input values
+% window_lgth: number of points in the windows used to define a mean and an std_dev
+% nb_std: number of std dev beyond which data are rejected
+% delta_lim = [delta_min delta_max] : interval of acceptable std_dev in window_lgth 
+%             (NaN possible if any value is accepted)
+% iter = number of iterations of the cleaning
+% fillval = value that will replace the value of the rejected data
+%
+% remarks: 
+%  - nb_std can be a scalar or a vector with the size of 1 x iter
+%  - delta_min is important in case the signal is constant over window_lgth (when std_dev=0)
+%  - delta_max id also useful when the signal is very bad
+%
+% examples:
+%    pressure = clean_median(pressure,20,[3 2.8],[1.5 10],2,-9.990e-29);
+%    temp1 = clean_median(temp1,12,3,[0.05 0.4],2,-9.990e-29);
+%    cond1 = clean_median(cond1,10,2.8,[0.01 0.4],3,-9.990e-29);
+%    oxy1 = clean_median(oxy1,10,2.8,[0.01 0.4],3,-9.990e-29);
+%
+% EX: P = [1:6 40 8:20]
+%  clean_median(P,5,3,[1.5 10],2,-100) : On travaille sur les donnees par paquet de 5 (=longueur de fenetre)
+%  ie les donnees [1:5], puis [6 40 8 9 10], puis [11:15] puis [16:20].
+% Pour chaque groupe de donnnees, on calcule la mediane et l'ecart-type sur
+% ce paquet de donnees auquel on ajoute npts donnees de part et autre; npts
+% etant la longueur de la demie-fenetre, soit ici npts = round(5/2) = 3.
+% Ainsi, dans ce cas, la mediane et l'ecart-type sont calcules sur 11 pts
+% au total (5(longueur de fenetre) + 2*npts = 5 + 2*3 = 11).
+% Chacun des 5 points de la fenetre est ensuite compare avec cette mediane et cet ecart-type
+% et est garde ou non.
+%
+% ??/??/2010 : P. Lherminier : Creation
+%       2017 : P. Rousselot  : Modification valeurs absolues -> par(abs(par)<abs(med-dmax) | abs(par)>abs(med+dmax))
+
+if length(nb_std) == 1,
+    nb_std = repmat(nb_std,[iter,1]);
+end
+perm_param = 0;
+
+size_param = size(param);
+if size_param(1) == 1,
+    param = param(:);
+    perm_param = 1;
+end
+deja_fillval = find(param==fillval); % On repere les donnees deja a badval.
+param(param==fillval) = NaN;
+nparam = length(param);
+nwin = floor(nparam/window_lgth)+1;
+nlast = nwin*window_lgth-nparam;
+
+par=[param;medianoutnan(param(end-nlast+1:end))*ones(nlast,1)];
+par = reshape(par,[window_lgth nwin]);
+% npts data are added before and after the window to calculate the std_dev (so that 
+% the extremes points of the segments are not eliminated in strong gradients)
+npts = round(window_lgth/2);
+
+for iiter = 1:iter,
+    parstd = [[nan(npts,1) par(end-npts+1:end,1:end-1)];par;[par(1:npts,2:end) nan(npts,1)]];
+    [stdm,med] = stdmedian(parstd,0,1);
+    med = med(npts+1:end-npts,:);
+    dmax = repmat(min(max(nb_std(iiter)*stdm,delta_lim(1)),delta_lim(2)),[window_lgth,1]);
+    par(abs(med-par)>=dmax) = NaN;            %%%%%%%%%%%%%%%%%%%%%%%
+end
+par = par(:); med = med(:); dmax = dmax(:);
+par = par(1:nparam); med = med(1:nparam); dmax = dmax(1:nparam);
+
+% Calcule le % de rejet apres le test d'ecart a le mediane.
+nouveau_fillval = find(isnan(par));
+info_rejet_med=(length(nouveau_fillval)-length(deja_fillval))*100/nparam;
+
+par(isnan(par)) = fillval;
+
+
+if perm_param == 1,
+    par = par';
+end
+
+%clf
+%plot(param,'b.-');
+%hold on;grid on
+%plot(par,'r.');
+%plot([med-dmax med+dmax],'c-');
+
+function [y,med] = stdmedian(x,flag,dim)
+%STD    Standard deviation from the median and not the mean.
+%   For vectors, STD(X) returns a pseudo standard deviation. For matrices,
+%   STD(X) is a row vector containing the standard deviation of each
+%   column.  For N-D arrays, STD(X) is the standard deviation of the
+%   elements along the first non-singleton dimension of X.
+%
+%   STD(X) normalizes by (N-1) where N is the sequence length.  This
+%   makes STD(X).^2 the best unbiased estimate of the variance if X
+%   is a sample from a normal distribution.
+%
+%   STD(X,1) normalizes by N and produces the second moment of the
+%   sample about its mean.  STD(X,0) is the same as STD(X).
+%
+%   STD(X,FLAG,DIM) takes the standard deviation along the dimension
+%   DIM of X.  When FLAG=0 STD normalizes by (N-1), otherwise STD
+%   normalizes by N.
+%
+%   [Y,M] = STD(X) returns the pseudo std Y and the median M
+%
+%   Example: If X = [4 -2 1
+%                    9  5 7]
+%     then std(X,0,1) is [ 3.5355 4.9497 4.2426] and std(X,0,2) is [3.0
+%                                                                   2.0]
+%   See also COV, MEAN, MEDIAN, CORRCOEF.
+
+%   J.N. Little 4-21-85
+%   Revised 5-9-88 JNL, 3-11-94 BAJ, 5-26-95 dlc, 5-29-96 CMT.
+%   Copyright (c) 1984-98 by The MathWorks, Inc.
+%   $Revision: 5.17 $  $Date: 1997/11/21 23:24:08 $
+
+if nargin<2, flag = 0; end
+if nargin<3,
+  dim = find(size(x)~=1, 1 );
+  if isempty(dim), dim = 1; end
+end
+
+% Avoid divide by zero.
+if size(x,dim)==1, y = zeros(size(x)); return, end
+
+tile = ones(1,max(ndims(x),dim));
+tile(dim) = size(x,dim);
+
+%xc = x - repmat(sum(x,dim)/size(x,dim),tile);  % Remove mean
+med = repmat(medianoutnan(x,dim),tile);
+xc = x - med;  % Remove median
+xcnan=(isnan(xc));
+coef=sum(~xcnan,dim);
+coef(coef==0)=NaN;
+xc(xcnan)=0;
+if flag,
+  y = sqrt(sum(conj(xc).*xc,dim)./coef);
+else
+  coef(coef==1)=2; %then y=0
+  y = sqrt(sum(conj(xc).*xc,dim)./(coef-1));
+end
+
+function y = medianoutnan(x,dim)
+%MEDIAN Median value.
+%   For vectors, MEDIANOUTNAN(X) is the median value of the finite
+%   elements in X. For matrices, MEDIANOUTNAN(X) is a row vector
+%   containing the median value of each column.  For N-D arrays,
+%   MEDIANOUTNAN(X) is the median value of the elements along the
+%   first non-singleton dimension of X.
+%
+%   MEDIANOUTNAN(X,DIM) takes the median along the dimension DIM of X.
+%
+%   Example: If X = [0 1 2   NaN
+%                    3 4 NaN NaN]
+%
+%   then medianoutnan(X,1) is [1.5 2.5 2 NaN]
+%   and medianoutnan(X,2) is [1
+%                             3.5]
+%
+%   See also MEANOUTNAN, STDOUTNAN, MIN, MAX, COV.
+
+%   P. Lherminier, 21/03/2002. From modifications of median.m
+
+if nargin==1,
+  dim = min(find(size(x)~=1));
+  if isempty(dim), dim = 1; end
+end
+if isempty(x), y = []; return, end
+
+siz = [size(x) ones(1,dim-ndims(x))];
+n = size(x,dim);
+
+% Permute and reshape so that DIM becomes the row dimension of a 2-D array
+perm = [dim:max(length(size(x)),dim) 1:dim-1];
+x = reshape(permute(x,perm),n,prod(siz)/n);
+
+% Sort along first dimension
+x = sort(x,1);
+sizx=size(x);
+
+[nnan,ncol]=find(diff(isnan(x)));
+nn=nan*ones(prod(siz)/n,1);
+y=nn';
+nn(ncol)=nnan;
+nn(~isnan(x(end,:)))=n;
+
+ieven=(rem(nn,2)==0);           % Even number of elements along DIM
+iodd=find(~ieven & ~isnan(nn)); % Odd number of elements along DIM
+ieven=find(ieven==1);
+x=x(:);
+if ~isempty(iodd),
+    y(iodd)=x(sub2ind(sizx,(nn(iodd)+1)/2,iodd));
+end;
+if ~isempty(ieven),
+    y(ieven)=(x(sub2ind(sizx,nn(ieven)/2,ieven))+ x(sub2ind(sizx,nn(ieven)/2+1,ieven)))/2;
+end;
+
+% Permute and reshape back
+siz(dim) = 1;
+y = ipermute(reshape(y,siz(perm)),perm);
+
+
diff --git a/filter/diurne.m b/filter/diurne.m
new file mode 100644
index 0000000..2a4beda
--- /dev/null
+++ b/filter/diurne.m
@@ -0,0 +1,176 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Biais capteurs - zoom cycle diurne                                      %
+% CNRS / Météo-France / Coriolis  --- AtlantOS                            %
+% Date: 01/08/2016  ---  Auteur: Pierre Rousselot                         %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%addpath('
+date = datestr(positions.timef(:,1));
+
+% Profil fin de nuit
+time =datestr(positions.timefinterp(:,1),'HH:MM');
+j=1;
+jj=1;
+for ii = 2:length(time(:,1))
+    %end_night(ii)=strcmp(time(ii,:),'03:00');
+    if strcmp(time(ii,:),'03:00')
+        end_night(j)=ii;
+        j=j+1;
+        if strcmp(time(ii+3,:),'06:00')
+            end_night_bis(jj)=ii+3;
+            jj=jj+1;
+        else
+            if strcmp(time(ii+2,:),'05:00')
+                end_night_bis(jj)=ii+2;
+                jj=jj+1;              
+            else
+                if strcmp(time(ii+1,:),'04:00')
+                    end_night_bis(jj)=ii+1;
+                    jj=jj+1;    
+                else
+                    if strcmp(time(ii,:),'03:00')
+                        end_night_bis(jj)=ii;
+                        jj=jj+1; 
+                    end
+                end
+            end
+        end
+    end
+
+end
+
+for ii= 1:331
+    mean_temp(ii,:) = mean(tempf(end_night(ii):end_night_bis(ii),:),1);
+    mean_depth(ii,:) = mean(depth(end_night(ii):end_night_bis(ii),:),1);
+    mean_time(ii,:) = mean(positions.timefinterp(end_night(ii):end_night_bis(ii),:),1);
+end
+
+% %Conv filtering
+% for ii=1:length(tempf(1,:))
+%     ttt(:,ii)=conv(tempf(end_night,ii),ones(1,3)/3,'valid');
+% end
+% 
+% %Median filtering
+% %tt=medfilt1(tempf(end_night,:));
+% hh=hann(24);
+% for ii=1:length(tempf(1,:))
+%     tt(:,ii)=conv(tempf(end_night,ii),hh,'valid');
+% end
+% 
+% 
+% %Filter
+% %tttt=filter(ones(1,3)/3,1,(tempf(end_night,:)));
+% 
+% %Zero-Phase filtering
+% ttttt=filtfilt(ones(1,3)/3,1,(tempf(end_night,:)));
+% 
+% %Savitky-Golay fitering
+% tttttt=sgolayfilt(tempf(:,:),3,7);
+% 
+% %LPO-ecart median
+% % for ii=1:length(tempf(1,:))
+% %     [par(:,ii),info_rejet_med]=clean_median(tempf(:,ii),6,2.8,[0.05 1],1,NaN);
+% % end
+% 
+% figure(1);
+% title('End-Night - Different filter')
+% subplot(5,1,1)
+% plot(tempf(end_night,:))
+% subplot(5,1,2)
+% plot(tt)
+% subplot(5,1,3)
+% plot(ttt)
+% subplot(5,1,4)
+% plot(ttttt)
+% subplot(5,1,5)
+% plot(tttttt)
+
+figure
+contourf(mean_time, -mean_depth, mean_temp);
+
+
+
+figure(2);
+datestr(positions.timefinterp(end_night,1));
+contourf(positions.timefinterp(end_night,:), -depth_ma(end_night,:), yyyy(:,end_night)')
+datetick('x','mm/yy','keepticks')
+cc = colorbar;
+xlabel('Time')
+ylabel('Depth [m]')
+ylabel(cc,'Temperature [{\circ}C]');
+caxis([min(min(tempf_interp)) max(max(tempf_interp))])
+grid on
+title('End-Night Mean (3h-6h) Time Series')
+
+% TT = 18;
+% %Zoom cycle diurne - biais capteur (fin de nuit)
+% figure(3);
+% plot(positions.timef(TT:TT+700,:),tempf(TT:TT+700,:))
+% datetick('x','HH dd/mm','keepticks')
+% grid on
+% j1m = mean(tempf(46:49,:));
+% j2m = mean(tempf(69:72,:));
+% diff_mean = j1m-j2m
+% title('Zoom 2 days')
+
+
+%Biais intercapteur
+% figure(4);
+% 
+% if options.marisonde
+%     %plot(positions.timef(500+24*3:500+24*6,:),tempf(500+24*3:500+24*6,:))
+%     plot(positions.timef(TT+24*3:TT+24*6,:),tempf(TT+24*3:TT+24*6,:))
+% else
+%     plot(positions.timef(5100+24*3:5100+24*6,:),tempf(5100+24*3:5100+24*6,:))
+% end
+% datetick('x','HH dd/mm','keepticks')
+% grid on
+% title('Zoom Homogeneous Period')
+
+
+% figure(5);
+% if options.marisonde
+%     %plot(tempf(500,:))
+%     plot(tempf(TT,:))
+%     hold on
+%     %plot(mean(tempf(500,:))*ones(1,17),'r--')
+%     plot(mean(tempf(TT,:))*ones(1,17),'r--')
+%     %max_ecart_intersensor = range(tempf(500,:))  
+%     max_ecart_intersensor = range(tempf(TT,:)) 
+% else
+%     plot(tempf(5178+24,:))
+%     hold on
+%     plot(mean(tempf(5178+24,:))*ones(1,17),'r--')
+%     max_ecart_intersensor = range(tempf(5178,:))
+% end
+% hold off
+% grid on
+% view([90 90])
+% title('Inter-sensor Bias')
+% 
+% %[maxbias_intersensor,ind_maxbias_intersensor]=max(std(tempf(4934:5280,:),[],2))
+% [maxbias_intersensor,ind_maxbias_intersensor]=max(std(tempf(TT:100,:),[],2))
+% 
+% % mm = mean(tempf(5100+24*3:5100+24*6,:));
+% % stdmm = std(mm);
+% mm = mean(tempf(TT+24*3:TT+24*6,:));
+% stdmm = std(mm);
+% 
+% %Derive dans le temps??
+% j=1;
+% %for ii=4950:7400 %1:length(time(:,1))
+% for ii=TT:TT+1100 %1:length(time(:,1))
+%     %end_night(ii)=strcmp(time(ii,:),'03:00');
+%     if strcmp(time(ii,:),'03:00')
+%         end_night2(j)=ii;
+%         j=j+1;
+%     end
+% end
+% mmean = nanmean(tempf(end_night2,:),2);
+% for ii = 1:length(tempf(1,:))
+% diff(:,ii) = tempf(end_night2,ii)-mmean;
+% end
+% mmmmm = nanmean(diff);
+% for ii=1:length(end_night2(1,:))
+% te(ii,:)=tempf(end_night2(1,ii),:)-mmmmm;
+% end
\ No newline at end of file
diff --git a/filter/fft_filter.m b/filter/fft_filter.m
new file mode 100644
index 0000000..1200681
--- /dev/null
+++ b/filter/fft_filter.m
@@ -0,0 +1,80 @@
+function [reconstructed] = fft_filter(data,dimension,bandlimits)
+%The filter is designed to filter gridded data in an array of up to 5 dimensions (DATA)
+%with the filtering done only in the dimension specified by DIMENSION.
+%The signal in data is first decomposed in the frequency domain (amplitude vs frequency).
+%Then the amplitudes associated with frequencies outside of the BANDLIMITS are set to zero.
+%Then an inverse Fourier transform is performed resulting in the reconstructed signal.
+ 
+% DATA : array to filter of up to 5 dimensions, the dimension to filter must be evenly spaced and be a multiple of 2.
+% if the dimension is spaced unevenly, interpolation to regular interval is advised.
+% DIMENSION : 1 integer indication the dimension to filter
+% BANDLIMITS : 2 numbers indicating the boundary of the selected frequency band. The bandlimit is actually
+% expressed as the period. Bandlimits have the same unit as the spacing between measurements
+% ex: if filtering with measurements every day, bandlimit of [5 10] means only oscillation with periods
+% of 5 to 10 days are kept in the signal. To construct a low-pass or high-pass, on can use [5 Inf] or [0 5]
+% RECONSTRUCTED : array of same size as data containing the filtered signal.
+ 
+% ATTENTION : Fourier transform are accurate for periodic signals. If the filtered values are not periodic,
+% it is advised to make them periodic. Ex signal= 1 2 5 coult be modified to 1 2 5 5 2 1 for periodicity.
+% Only the first half of reconstructed should be used in this case.
+ 
+%Bring dimension to filter to dimension 1
+si_data=size(data);
+order=1:length(si_data);
+otherdim=order(ismember(order,dimension)==0);
+neworder=[dimension,otherdim];
+data=permute(data,[dimension,otherdim]);
+ 
+%Fourier transform and associated frequencies
+amplitudes = fft(data,[],1);
+n1=size(data);
+n=size(amplitudes,1); %replaced dimension by 1
+nyquist=1/2;
+nyquist_location=n/2+1;
+ 
+%Generate a list of the frequencies
+if rem(size(data,1),2)==0;
+frequencies=[0,(1:n/2-1)/n,n/2/n,(n/2-1:-1:1)/n];
+end
+if rem(size(data,1),2)==1;
+frequencies=[0,(1:(n-1)/2)/n,((n-1)/2:-1:1)/n];
+end
+ 
+%Computes periods and decide of the ones to keep
+periods=1./frequencies;
+periodstokeep=(periods>=bandlimits(1) & periods<=bandlimits(2));
+periods(periodstokeep);
+filteredamps=amplitudes*0;
+ 
+%Only keep desired periods (amplitudes)
+if ndims(data)==1
+filteredamps(periodstokeep)=amplitudes(periodstokeep);
+end
+if ndims(data)==2
+filteredamps(periodstokeep,:)=amplitudes(periodstokeep,:);
+end
+if ndims(data)==3
+filteredamps(periodstokeep,:,:)=amplitudes(periodstokeep,:,:);
+end
+if ndims(data)==4
+filteredamps(periodstokeep,:,:,:)=amplitudes(periodstokeep,:,:,:);
+end
+if ndims(data)==5
+filteredamps(periodstokeep,:,:,:,:)=amplitudes(periodstokeep,:,:,:,:);
+end
+ 
+%Reconstruct the signal
+reconstructed=ifft(filteredamps,[],1);
+reconstructedr=real(reconstructed);
+reconstructedi=imag(reconstructed);
+ 
+%Return matrix in original order
+ordereturn=zeros([1,length(si_data)]);
+ordereturn(dimension)=1;
+iszero=find(ordereturn==0);
+for i=1:length(iszero)
+ordereturn(iszero(i))=i+1;
+end
+reconstructed=permute(reconstructed,ordereturn);
+ 
+end
\ No newline at end of file
diff --git a/filter/filterfft.m b/filter/filterfft.m
new file mode 100644
index 0000000..4eb546b
--- /dev/null
+++ b/filter/filterfft.m
@@ -0,0 +1,71 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Plot filtering data                                                     %
+% CNRS / Météo-France / Coriolis  --- AtlantOS                            %
+% Date: 14/09/2016  ---  Auteur: Pierre Rousselot                         %
+% fft_filter: Patrick Martineau                                           %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% tempf(end+1:length(tempf)*2,:)=tempf;
+% positions.timef(end+1:length(positions.timef)*2,:)=positions.timef;
+% depth(end+1:length(depth)*2,:)=depth;
+% positions.timef = linspace(1,15940,15940);
+% for i = 2:17
+%     positions.timef(i,:)=positions.timef(1,:);
+% end
+
+[reconstructed2]=fft_filter(tempf,1,[0 24]);
+
+figure
+plot(positions.timef,reconstructed)
+datetick('x','mm/yy','keepticks')
+xlabel('Time')
+ylabel('Temperature [degC]')
+%axis([min(positions.timef(:,1)) max(positions.timef(:,1)) 11 22])
+grid on
+title('Filtering process diff (<1week) - (<1day)')
+
+% figure;
+% contourf(positions.timef', -depth, reconstructed)
+
+figure
+subplot(4,1,1)
+contourf(vec_mean,-depth_mean,temp_mean)
+axis([7.34313e5 7.34318e5 -83 0])
+datetick('x','dd/mm','keeplimits')
+xlabel('Time')
+ylabel('Depth [m]')
+cc=colorbar;
+ylabel(cc,'Temperature [{\circ}C]')
+title('Daily mean')
+
+subplot(4,1,2)
+contourf(positions.timef,-depth,tempf)
+axis([7.34313e5 7.34318e5 -83 0])
+datetick('x','dd/mm','keeplimits')
+xlabel('Time')
+ylabel('Depth [m]')
+cc=colorbar;
+ylabel(cc,'Temperature [{\circ}C]')
+title('Raw Data')
+
+subplot(4,1,3)
+contourf(positions.timef,-depth,reconstructed)
+axis([7.34313e5 7.34318e5 -83 0])
+datetick('x','dd/mm','keeplimits')
+xlabel('Time')
+ylabel('Depth [m]')
+cc=colorbar;
+ylabel(cc,'Temperature [{\circ}C]')
+title('Filtered Data (>24h)')
+
+subplot(4,1,4)
+contourf(positions.timef,-depth,reconstructed2)
+axis([7.34313e5 7.34318e5 -83 0])
+xlabel('Time')
+ylabel('Depth [m]')
+cc=colorbar;
+ylabel(cc,'Temperature Anomaly [{\circ}C]')
+title('Diurnal Temperature Anomaly (<24h)')
+
+
+
diff --git a/filter/filtfilt.m b/filter/filtfilt.m
new file mode 100644
index 0000000..d4e2acd
--- /dev/null
+++ b/filter/filtfilt.m
@@ -0,0 +1,310 @@
+    function y = filtfilt(b,a,x)
+%FILTFILT Zero-phase forward and reverse digital IIR filtering.
+%   Y = FILTFILT(B, A, X) filters the data in vector X with the filter
+%   described by vectors A and B to create the filtered data Y.  The
+%   filter is described by the difference equation:
+%
+%     a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
+%                           - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
+%
+%   The length of the input X must be more than three times
+%   the filter order, defined as max(length(B)-1,length(A)-1). 
+%
+%   Y = FILTFILT(SOS, G, X) filters the data in vector X with the
+%   second-order section (SOS) filter described by the matrix SOS and the
+%   vector G.  The coefficients of the SOS matrix must be expressed using
+%   an Lx6 matrix where L is the number of second-order sections. The scale
+%   values of the filter must be expressed using the vector G. The length
+%   of G must be between 1 and L+1 and the length of input X must be more
+%   than 6 samples. The SOS matrix should have the following form:
+%
+%   SOS = [ b01 b11 b21 a01 a11 a21
+%           b02 b12 b22 a02 a12 a22
+%           ...
+%           b0L b1L b2L a0L a1L a2L ]
+%
+%   After filtering in the forward direction, the filtered sequence is then
+%   reversed and run back through the filter; Y is the time reverse of the
+%   output of the second filtering operation.  The result has precisely
+%   zero phase distortion, and magnitude modified by the square of the
+%   filter's magnitude response.  Startup and ending transients are
+%   minimized by matching initial conditions.
+%
+%   Note that FILTFILT should not be used when the intent of a filter is
+%   to modify signal phase, such as differentiators and Hilbert filters.
+%
+%   % Example:
+%   %   Zero-phase filter a noisy ECG waveform using an IIR filter.
+%  
+%   x = ecg(500)'+0.25*randn(500,1);        % noisy waveform
+%   [b,a] = butter(12,0.2,'low');           % IIR filter design
+%   y = filtfilt(b,a,x);                    % zero-phase filtering
+%   y2 = filter(b,a,x);                     % conventional filtering
+%   plot(x,'k-.'); grid on ; hold on
+%   plot([y y2],'LineWidth',1.5);
+%   legend('Noisy ECG','Zero-phase Filtering','Conventional Filtering');
+%
+%   See also FILTER, SOSFILT.
+
+%   References:
+%     [1] Sanjit K. Mitra, Digital Signal Processing, 2nd ed.,
+%         McGraw-Hill, 2001
+%     [2] Fredrik Gustafsson, Determining the initial states in forward-
+%         backward filtering, IEEE Transactions on Signal Processing,
+%         pp. 988-992, April 1996, Volume 44, Issue 4
+
+%   Copyright 1988-2011 The MathWorks, Inc.
+%   $Revision: 1.7.4.10 $  $Date: 2012/12/25 21:34:47 $
+
+error(nargchk(3,3,nargin,'struct'))
+
+% Only double precision is supported
+if ~isa(b,'double') || ~isa(a,'double') || ~isa(x,'double')
+    error(message('signal:filtfilt:NotSupported'));
+end
+if isempty(b) || isempty(a) || isempty(x)
+    y = [];
+    return
+end
+
+% If input data is a row vector, convert it to a column
+isRowVec = size(x,1)==1;
+if isRowVec
+    x = x(:);
+end
+[Npts,Nchans] = size(x);
+
+% Parse SOS matrix or coefficients vectors and determine initial conditions
+[b,a,zi,nfact,L] = getCoeffsAndInitialConditions(b,a,Npts);
+
+% Filter the data
+if Nchans==1
+    if Npts<10000
+        y = ffOneChanCat(b,a,x,zi,nfact,L);
+    else
+        y = ffOneChan(b,a,x,zi,nfact,L);
+    end
+else
+    y = ffMultiChan(b,a,x,zi,nfact,L);
+end
+
+if isRowVec
+    y = y.';   % convert back to row if necessary
+end
+
+%--------------------------------------------------------------------------
+function [b,a,zi,nfact,L] = getCoeffsAndInitialConditions(b,a,Npts)
+
+[L, ncols] = size(b);
+na = numel(a);
+
+% Rules for the first two inputs to represent an SOS filter:
+% b is an Lx6 matrix with L>1 or,
+% b is a 1x6 vector, its 4th element is equal to 1 and a has less than 2
+% elements. 
+if ncols==6 && L==1 && na<=2,
+    if b(4)==1,
+        warning(message('signal:filtfilt:ParseSOS', 'SOS', 'G'));
+    else
+        warning(message('signal:filtfilt:ParseB', 'a01', 'SOS'));
+    end
+end
+issos = ncols==6 && (L>1 || (b(4)==1 && na<=2));
+if issos,
+    %----------------------------------------------------------------------
+    % b is an SOS matrix, a is a vector of scale values
+    %----------------------------------------------------------------------
+    g = a(:);
+    ng = na;
+    if ng>L+1,
+        error(message('signal:filtfilt:InvalidDimensionsScaleValues', L + 1));
+    elseif ng==L+1,
+        % Include last scale value in the numerator part of the SOS Matrix
+        b(L,1:3) = g(L+1)*b(L,1:3);
+        ng = ng-1;
+    end
+    for ii=1:ng,
+        % Include scale values in the numerator part of the SOS Matrix
+        b(ii,1:3) = g(ii)*b(ii,1:3);
+    end
+    a = b(:,4:6).';
+    b = b(:,1:3).';
+    
+    nfact = 6;  % length of edge transients for each biquad
+    if Npts <= nfact      % input data too short
+        error(message('signal:filtfilt:InvalidDimensionsData'));
+    end
+    
+    % Compute initial conditions to remove DC offset at beginning and end of
+    % filtered sequence.  Use sparse matrix to solve linear system for initial
+    % conditions zi, which is the vector of states for the filter b(z)/a(z) in
+    % the state-space formulation of the filter.
+    zi = zeros(2,L);
+    for ii=1:L,
+        rhs  = (b(2:3,ii) - b(1,ii)*a(2:3,ii));
+        zi(:,ii) = ( eye(2) - [-a(2:3,ii),[1;0]] ) \ rhs;
+    end
+    
+else
+    %----------------------------------------------------------------------
+    % b and a are vectors that define the transfer function of the filter
+    %----------------------------------------------------------------------
+    L = 1;
+    % Check coefficients
+    b = b(:);
+    a = a(:);
+    nb = numel(b);
+    nfilt = max(nb,na);
+    nfact = 3*(nfilt-1);  % length of edge transients
+    if Npts <= nfact      % input data too short
+        error(message('signal:filtfilt:InvalidDimensionsDataShortForFiltOrder'));
+    end
+    % Zero pad shorter coefficient vector as needed
+    if nb < nfilt
+        b(nfilt,1)=0;
+    elseif na < nfilt
+        a(nfilt,1)=0;
+    end
+    
+    % Compute initial conditions to remove DC offset at beginning and end of
+    % filtered sequence.  Use sparse matrix to solve linear system for initial
+    % conditions zi, which is the vector of states for the filter b(z)/a(z) in
+    % the state-space formulation of the filter.
+    if nfilt>1
+        rows = [1:nfilt-1, 2:nfilt-1, 1:nfilt-2];
+        cols = [ones(1,nfilt-1), 2:nfilt-1, 2:nfilt-1];
+        vals = [1+a(2), a(3:nfilt).', ones(1,nfilt-2), -ones(1,nfilt-2)];
+        rhs  = b(2:nfilt) - b(1)*a(2:nfilt);
+        zi   = sparse(rows,cols,vals) \ rhs;
+        % The non-sparse solution to zi may be computed using:
+        %      zi = ( eye(nfilt-1) - [-a(2:nfilt), [eye(nfilt-2); ...
+        %                                           zeros(1,nfilt-2)]] ) \ ...
+        %          ( b(2:nfilt) - b(1)*a(2:nfilt) );
+    else
+        zi = zeros(0,1);
+    end
+end
+%--------------------------------------------------------------------------
+function y = ffOneChanCat(b,a,y,zi,nfact,L)
+
+for ii=1:L
+    % Single channel, data explicitly concatenated into one vector
+    y = [2*y(1)-y(nfact+1:-1:2); y; 2*y(end)-y(end-1:-1:end-nfact)]; %#ok<AGROW>
+    
+    % filter, reverse data, filter again, and reverse data again
+    y = filter(b(:,ii),a(:,ii),y,zi(:,ii)*y(1));
+    y = y(end:-1:1);
+    y = filter(b(:,ii),a(:,ii),y,zi(:,ii)*y(1));
+    
+    % retain reversed central section of y
+    y = y(end-nfact:-1:nfact+1);
+end
+
+%--------------------------------------------------------------------------
+function y = ffOneChan(b,a,xc,zi,nfact,L)
+% Perform filtering of input data with no phase distortion
+%
+% xc: one column of input data
+% yc: one column of output, same dimensions as xc
+% a,b: IIR coefficients, both of same order/length
+% zi: initial states
+% nfact: scalar
+% L: scalar
+
+% Extrapolate beginning and end of data sequence using reflection.
+% Slopes of original and extrapolated sequences match at end points,
+% reducing transients.
+%
+% yc is length numel(x)+2*nfact
+%
+% We wish to filter the appended sequence:
+% yc = [2*xc(1)-xc(nfact+1:-1:2); xc; 2*xc(end)-xc(end-1:-1:end-nfact)];
+%
+% We would use the following code:
+% Filter the input forwards then again in the backwards direction
+%   yc = filter(b,a,yc,zi*yc(1));
+%   yc = yc(length(yc):-1:1); % reverse the sequence
+%
+% Instead of reallocating and copying, just do filtering in pieces
+% Filter the first part, then xc, then the last part
+% Retain each piece, feeding final states as next initial states
+% Output is yc = [yc1 yc2 yc3]
+%
+% yc2 can be long (matching user input length)
+% yc3 is short (nfilt samples)
+%
+for ii=1:L
+    
+    xt = -xc(nfact+1:-1:2) + 2*xc(1);
+    [~,zo] = filter(b(:,ii),a(:,ii), xt, zi(:,ii)*xt(1)); % yc1 not needed
+    [yc2,zo] = filter(b(:,ii),a(:,ii), xc, zo);
+    xt = -xc(end-1:-1:end-nfact) + 2*xc(end);
+    yc3 = filter(b(:,ii),a(:,ii), xt, zo);
+    
+    % Reverse the sequence
+    %   yc = [flipud(yc3) flipud(yc2) flipud(yc1)]
+    % which we can do as
+    %   yc = [yc3(end:-1:1); yc2(end:-1:1); yc1(end:-1:1)];
+    %
+    % But we don't do that explicitly.  Instead, we wish to filter this
+    % reversed result as in:
+    %   yc = filter(b,a,yc,zi*yc(1));
+    % where yc(1) is now the last sample of the (unreversed) yc3
+    %
+    % Once again, we do this in pieces:
+    [~,zo] = filter(b(:,ii),a(:,ii), yc3(end:-1:1), zi(:,ii)*yc3(end));
+    yc5 = filter(b(:,ii),a(:,ii), yc2(end:-1:1), zo);
+    
+    % Conceptually restore the sequence by reversing the last pieces
+    %    yc = yc(length(yc):-1:1); % restore the sequence
+    % which is done by
+    %    yc = [yc6(end:-1:1); yc5(end:-1:1); yc4(end:-1:1)];
+    %
+    % However, we only need to retain the reversed central samples of filtered
+    % output for the final result:
+    %    y = yc(nfact+1:end-nfact);
+    %
+    % which is the reversed yc5 piece only.
+    %
+    % This means we don't need yc4 or yc6.  We need to compute yc4 only because
+    % the final states are needed for yc5 computation.  However, we can omit
+    % yc6 entirely in the above filtering step.
+    xc = yc5(end:-1:1);
+end
+y = xc;
+
+%--------------------------------------------------------------------------
+function y = ffMultiChan(b,a,xc,zi,nfact,L)
+% Perform filtering of input data with no phase distortion
+%
+% xc: matrix of input data
+% yc: matrix of output data, same dimensions as xc
+% a,b: IIR coefficients, both of same order/length
+% zi: initial states
+% nfact: scalar
+% L: scalar
+
+% Same comments as in ffOneChan, except here we need to use bsxfun.
+% Instead of doing scalar subtraction with a vector, we are doing
+% vector addition with a matrix.  bsxfun replicates the vector
+% for us.
+%
+% We also take care to preserve column dimensions
+%
+sz = size(xc);
+xc = reshape(xc,sz(1),[]);%converting N-D matrix to 2-D.
+
+for ii=1:L
+    xt = bsxfun(@minus, 2*xc(1,:), xc(nfact+1:-1:2,:));
+    [~,zo] = filter(b(:,ii),a(:,ii),xt,zi(:,ii)*xt(1,:)); % outer product
+    [yc2,zo] = filter(b(:,ii),a(:,ii),xc,zo);
+    xt = bsxfun(@minus, 2*xc(end,:), xc(end-1:-1:end-nfact,:));
+    yc3 = filter(b(:,ii),a(:,ii),xt,zo);
+    
+    [~,zo] = filter(b(:,ii),a(:,ii),yc3(end:-1:1,:),zi(:,ii)*yc3(end,:)); % outer product
+    yc5 = filter(b(:,ii),a(:,ii),yc2(end:-1:1,:), zo);
+    
+    xc = yc5(end:-1:1,:);
+end
+y = xc;
+y = reshape(y,sz);% To match the input size.
diff --git a/filter/hamm.m b/filter/hamm.m
new file mode 100644
index 0000000..3a8c5ab
--- /dev/null
+++ b/filter/hamm.m
@@ -0,0 +1,14 @@
+function w = hamm(n)
+%HAMM HAMM(N) returns the N-point Hamming window.
+
+%   Copyright (c) 1988-96 by The MathWorks, Inc.
+%       $Revision: 1.7 $  $Date: 1996/07/25 16:36:49 $
+
+if n > 1
+    w = .54 - .46*cos(2*pi*(0:n-1)'/(n-1));
+elseif n == 1
+    w = .08;
+else
+    error('N must be greater than 0.')
+end
+
diff --git a/filter/hammfilter_nodec.m b/filter/hammfilter_nodec.m
new file mode 100644
index 0000000..d3d616d
--- /dev/null
+++ b/filter/hammfilter_nodec.m
@@ -0,0 +1,29 @@
+function uf = hammfilter_nodec(u,aver)
+
+% run a box filter of length aver (points) over time series u
+% replaces lowpass filter if toolbox license is missing
+% no decimation (subsampling)
+%
+% -------- R. Zantopp  24 Nov 2006 
+
+nn = length(u);
+na = floor(aver/2);
+% [w]=hamming(aver);
+w  = hamm(aver);       % weights
+uu = reshape(u,1,length(u));
+
+for i = 1:nn-aver
+    u1  = uu(i:i+aver-1);
+    au  = u1'.*w;
+    nau = sum(isnan(au));
+    if nau<=floor(aver/2) 
+      au_good  = find(~isnan(au));
+      au2      = meanmiss(au)/meanmiss(w(au_good));
+      u2(i+na) = au2;
+    else
+      u2(i+na) = NaN;
+    end;
+end;
+u2(1:na)     = NaN;
+u2(nn-na:nn) = NaN;
+uf           = u2;
diff --git a/filter/hamming.m b/filter/hamming.m
new file mode 100644
index 0000000..b412338
--- /dev/null
+++ b/filter/hamming.m
@@ -0,0 +1,36 @@
+
+% HAMMING.M	- returns the N-point Hamming window.
+%
+% Input		: n = number
+%
+% Output	: w = vector
+%
+% Usage		: w = hamming (n)
+%
+% Comments	: allows also the call:  hamming(xx), taking only format from signal xx
+%
+% See also	: HAMMING2
+
+% modification of original MATLAB (3.5)  file
+% HGFei, 1990  
+
+function w = hamming(n)
+
+if nargin == 0
+ help hamming;
+ return;
+end;
+
+n = n(:); 
+[hig, wid] = size(n);
+
+if  wid > 1;  n = wid; end; 
+
+w = (.54 - .46*cos(2*pi*(0:n-1)'/(n-1))).';
+
+if  hig > 1;
+ww = w;  
+ for jj = 2 : hig;
+ w(jj,:) =  ww;
+ end; 
+end; 
diff --git a/filter/lisse.m b/filter/lisse.m
new file mode 100644
index 0000000..fadfd79
--- /dev/null
+++ b/filter/lisse.m
@@ -0,0 +1,131 @@
+function [yl,yy] = lisse(y0,tab,typ,dk,fillval)
+%
+% Yl = lisse(Y,TAB) is a 1-D running window filter. Y = Y(i) is the signal to be filtered, 
+% TAB is a vector containing the non-zero weights defining the shape of the window [i-n i+n],  
+% where n*2+1 = length(TAB), n being the half-width of the window, and length(TAB) must be odd.  
+% TAB allows to define manually the window size and shape. If the weights provided in TAB do not sum to 1, 
+% they are normalized so that the energy of the signal is conserved. 
+% 
+% Yl = lisse(Y,TAB,typ) is equivalent to Yl = lisse(Y,TAB) if typ = 'mean'. 
+% Using typ = 'min'('max') allows to select the minimum (maximum) value over the running window 
+% rather than apply the weights.  
+%
+% For the n first and last points, the signal is only partially defined over the window. 
+% Yl = lisse(Y,TAB,typ,dk,fillval) permits to provide a fill-value to pad the signal 
+% over the window for these first and last points. 
+% If fillval is the string 'extrapol', Y(1) or Y(end) are used for padding. 
+% If fillval is a scalar, it is used for padding anywhere. 
+% If fillval is a pair of scalars, the first value is used for padding at the beginning, 
+% the second one at the end.
+%
+% If dk differs from 0, the signal is further shifted forward or backward and padded 
+% with the corresponding fill-value, so that 
+%    Yl(k) = Yl(k+dk)
+% 
+y = y0(:)'; 
+if sum(tab)~=1, tab = tab/sum(tab); end
+if nargin<5, 
+  fillval_beg = NaN; 
+  fillval_end = NaN; 
+else
+  if isstr(fillval),
+    if strcmp(fillval,'extrapol'), 
+      fillval_beg = y(1); fillval_end = y(end); 
+    end
+  else
+    switch length(fillval)
+     case 1, fillval_beg = fillval; fillval_end = fillval; 
+     case 2, fillval_beg = fillval(1); fillval_end = fillval(2); 
+    end
+  end
+end
+if nargin<4, dk = 0; end
+if nargin<3 | isempty(typ), typ = 'mean'; end
+%if strcmp(typ,'max') | strcmp(typ,'min'), tab = ones(size(tab)); end
+%
+in = floor((length(tab)+1)/2); halfwidth = in-1; 
+
+if size(y0,1)==1, % y0 is a vector
+  %
+  yy = zeros(length(tab),length(y));
+  %
+  switch typ
+   case 'mean',
+    yy(in,:) = tab(in)*y; 
+    for k = 1:in-1
+      decneg = in-k; 
+      yy(k,:) = tab(k)*[fillval_beg*ones(1,decneg) y(1:end-decneg)]; 
+    end
+    for k = in+1:length(tab)
+      decpos = k-in; 
+      yy(k,:) = tab(k)*[y(1+decpos:end) fillval_end*ones(1,decpos)]; 
+    end
+    yl = nansum(yy); 
+    %
+   case 'mean2',
+    yy(in,:) = tab(in)*y; 
+    for k = 1:in-1
+      decneg = in-k; 
+      yy(k,:) = tab(k)*[fillval_beg*ones(1,decneg) y(1:end-decneg)].^2; 
+    end
+    for k = in+1:length(tab)
+      decpos = k-in; 
+      yy(k,:) = tab(k)*[y(1+decpos:end) fillval_end*ones(1,decpos)].^2; 
+    end
+    yl = sqrt(nansum(yy)); 
+    %
+   case 'median',
+    yy(in,:) = y; 
+    for k = 1:in-1
+      decneg = in-k; 
+      yy(k,:) = [fillval_beg*ones(1,decneg) y(1:end-decneg)]; 
+    end
+    for k = in+1:length(tab)
+      decpos = k-in; 
+      yy(k,:) = [y(1+decpos:end) fillval_end*ones(1,decpos)]; 
+    end
+    yl = nanmedian(yy); 
+    %
+   case 'max',
+    yy(in,:) = y; 
+    for k = 1:in-1
+      decneg = in-k; 
+      yy(k,:) = [fillval_beg*ones(1,decneg) y(1:end-decneg)]; 
+    end
+    for k = in+1:length(tab)
+      decpos = k-in; 
+      yy(k,:) = [y(1+decpos:end) fillval_end*ones(1,decpos)]; 
+    end
+    yl = nanmax(yy); 
+    %
+   case 'min',
+    yy(in,:) = y;
+    for k = 1:in-1
+      decneg = in-k;
+      yy(k,:) = [fillval_beg*ones(1,decneg) y(1:end-decneg)];
+    end
+    for k = in+1:length(tab)
+      decpos = k-in;
+      yy(k,:) = [y(1+decpos:end) fillval_end*ones(1,decpos)];
+    end
+    yl = nanmin(yy);
+    %
+  end
+else          % y0 is a matrice
+  'coucou'
+end
+%
+y0 = yl; 
+if dk == 0, yl = y0; 
+elseif dk < 0, 
+%keyboard
+%pause
+  yl = [fillval_beg*ones(1,-dk) y0(1:end+dk)]; 
+elseif dk > 0, 
+  yl = [y0(1+dk:end) fillval_end*ones(1,dk)]; 
+end
+
+%keyboard
+%pause
+
+return
diff --git a/filter/medfilt1_perso.m b/filter/medfilt1_perso.m
new file mode 100644
index 0000000..3dbe693
--- /dev/null
+++ b/filter/medfilt1_perso.m
@@ -0,0 +1,22 @@
+function y = medfilt1_perso(x,w)
+%MEDFILT1_PERSO 
+%Gere la non presence de signal processing toolbox
+%
+%odd w
+if(mod(w,1)==1)
+    w=w+1;
+end
+%
+n=length(x);
+for i=1:n
+    if(i<=w/2)
+       y(i)=median(x(1:i+w/2)); 
+    elseif(i>=n-w/2)
+       y(i)=median(x(i-w/2:end)); 
+    else
+       y(i)=median(x(i-w/2:i+w/2)); 
+    end
+end
+
+end
+
diff --git a/filter/median_filter.m b/filter/median_filter.m
new file mode 100644
index 0000000..619a3a6
--- /dev/null
+++ b/filter/median_filter.m
@@ -0,0 +1,106 @@
+function y=median_filter(x,n,blksz,DIM)
+%MEDFILT1  One dimensional median filter.
+%   Y = MEDFILT1(X,N) returns the output of the order N, one dimensional
+%   median filtering of X.  Y is the same size as X; for the edge points,
+%   zeros are assumed to the left and right of X.  If X is a matrix,
+%   then MEDFILT1 operates along the columns of X.
+%
+%   If you do not specify N, MEDFILT1 uses a default of N = 3.
+%   For N odd, Y(k) is the median of X( k-(N-1)/2 : k+(N-1)/2 ).
+%   For N even, Y(k) is the median of X( k-N/2 : k+N/2-1 ).
+%
+%   Y = MEDFILT1(X,N,BLKSZ) uses a for-loop to compute BLKSZ ("block size") 
+%   output samples at a time.  Use this option with BLKSZ << LENGTH(X) if 
+%   you are low on memory (MEDFILT1 uses a working matrix of size
+%   N x BLKSZ).  By default, BLKSZ == LENGTH(X); this is the fastest
+%   execution if you have the memory for it.
+%
+%   For matrices and N-D arrays, Y = MEDFILT1(X,N,[],DIM) or 
+%   Y = MEDFILT1(X,N,BLKSZ,DIM) operates along the dimension DIM.
+%
+%   See also MEDIAN, FILTER, SGOLAYFILT, and MEDFILT2 in the Image
+%   Processing Toolbox.
+
+%   Author(s): L. Shure and T. Krauss, 8-3-93
+%   Copyright 1988-2004 The MathWorks, Inc.
+%   $Revision: 1.8.4.2 $  $Date: 2004/12/26 22:16:21 $
+
+% Validate number of input arguments
+narginchk(1,4)
+if nargin < 2, n = []; end
+if nargin < 3, blksz = []; end
+if nargin < 4, DIM = []; end
+
+% Check if the input arguments are valid
+if isempty(n)
+  n = 3;
+end
+
+if ~isempty(DIM) && DIM > ndims(x)
+	error('Dimension specified exceeds the dimensions of X.')
+end
+
+% Reshape x into the right dimension.
+if isempty(DIM)
+	% Work along the first non-singleton dimension
+	[x, nshifts] = shiftdim(x);
+else
+	% Put DIM in the first (row) dimension (this matches the order 
+	% that the built-in filter function uses)
+	perm = [DIM,1:DIM-1,DIM+1:ndims(x)];
+	x = permute(x,perm);
+end
+
+% Verify that the block size is valid.
+siz = size(x);
+if isempty(blksz),
+	blksz = siz(1); % siz(1) is the number of rows of x (default)
+else
+	blksz = blksz(:);
+end
+
+% Initialize y with the correct dimension
+y = zeros(siz); 
+
+% Call medfilt1D (vector)
+for i = 1:prod(siz(2:end)),
+	y(:,i) = medfilt1D(x(:,i),n,blksz);
+end
+
+% Convert y to the original shape of x
+if isempty(DIM)
+	y = shiftdim(y, -nshifts);
+else
+	y = ipermute(y,perm);
+end
+
+
+%-------------------------------------------------------------------
+%                       Local Function
+%-------------------------------------------------------------------
+function y = medfilt1D(x,n,blksz)
+%MEDFILT1D  One dimensional median filter.
+%
+% Inputs:
+%   x     - vector
+%   n     - order of the filter
+%   blksz - block size
+
+nx = length(x);
+if rem(n,2)~=1    % n even
+    m = n/2;
+else
+    m = (n-1)/2;
+end
+X = [zeros(m,1); x; zeros(m,1)];
+y = zeros(nx,1);
+
+% Work in chunks to save memory
+indr = (0:n-1)';
+indc = 1:nx;
+for i=1:blksz:nx
+    ind = indc(ones(1,n),i:min(i+blksz-1,nx)) + ...
+          indr(:,ones(1,min(i+blksz-1,nx)-i+1));
+    xx = reshape(X(ind),n,min(i+blksz-1,nx)-i+1);
+    y(i:min(i+blksz-1,nx)) = median(xx,1);
+end
\ No newline at end of file
diff --git a/filter/mfilter.m b/filter/mfilter.m
new file mode 100644
index 0000000..654e050
--- /dev/null
+++ b/filter/mfilter.m
@@ -0,0 +1,1075 @@
+function [y,c]=mfilter(x,f0,f1,f2,n,m);
+
+% MFILTER  FIR digital filter
+%
+% [Y,C] = MFILTER(X,F0,F1,F2,N,M);
+%
+% filter data vector X by N'th order FIR1 filter 
+% Hamming window second order used X timeseries along 2. Dimension !!!
+%
+% F0 sample frequency
+%
+% F1 , F2 filter frequency
+%
+% F1=0       F2=highpass
+% F1=lowpass F2=0
+%
+% F1 < F2 bandpass
+% F1 > F2 bandstopp
+%
+% N Number of elements to filter, default: F0/max(F1,F2)
+%
+% M Increment for Subsample X after filter, default: 1
+%
+% C Filter-Coefficients by FIR1 with length  N+1
+%
+%
+% see also: BFILTER, MEANIND1, MEANIND2, NOISE
+%           FIR1, FILTFILT (SignalProcessingToolbox)
+%
+%-------------------------------------------------------------
+% Example for LowPass:
+%
+%  x = linspace(0,50,1000);
+%  y = sin(x);               % RowVector !!!
+%
+%  dx  = median(diff(x));    % Inkrement
+%  int = 4*pi;
+%
+%  c   = round(int/dx);
+%  z   = mfilter( y , 1 , 1/c , 0 , c , 1 );
+%
+%  figure,plot(x,y),hold on,plot(x,z,'r')
+%
+
+%M. Visbeck 29.07.91
+
+%---------------------------------------------------------------------
+% SubFunctions from Signal and Matlab-Toolboxes
+%---------------------------------------------------------------------
+%
+% FILTFILT (FILT2) Zero-phase forward and reverse digital filtering.
+% y = filtfilt(b,a,x)
+%
+% FIR1   FIR filter design using the window method.
+% [b,a] = fir1(N,Wn,varargin)
+%
+% FIRCHK   Check if specified filter order is valid.
+% [n,msg1,msg2] = firchk(n,Fend,a,exception)
+%
+% FIRLS Linear-phase FIR filter design using least-squares error minimization.
+% [h,a]=firls(N,F,M,W,ftype);
+%
+% SINC Sin(pi*x)/(pi*x) function.
+% y=sinc(x)
+%
+% HAMMING   Hamming window.
+% w = hamming(varargin)
+%
+% GENCOSWIN   Returns one of the generalized cosine windows.
+% [w,msg] = gencoswin(varargin)
+%
+%---------------------------------------------------------------------
+
+if nargin < 5, n=ceil(f0/max(f1,f2)); end
+if nargin < 5, m=1; end
+
+%lp
+if     f2==0, Wn = f1;      typ=''; 
+%hp
+elseif f1==0, Wn = f2;      typ='high'; 
+%bp
+elseif f2>f1, Wn = [f1 f2]; typ=''; 
+%bs
+elseif f1>f2, Wn = [f2 f1]; typ='stop'; 
+end
+
+Wn = 2 * Wn / f0;
+
+b=fir1(n,Wn,typ);
+
+si=size(x);
+
+if any(si([1 2])) == 1
+   y=filt2(b,1,x);
+   y=y(1:m:end);
+else
+   y=x;
+   for ii=1:si(1)
+       y(ii,:)=filt2(b,1,x(ii,:));
+   end
+   y=y(:,1:m:si(2));
+end
+
+%**************************************************************
+%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+function y = filt2(b,a,x)
+%FILT2 Zero-phase forward and reverse digital filtering.
+%   Y = FILT2(B, A, X) filters the data in vector X with the filter described
+%   by vectors A and B to create the filtered data Y.  The filter is described 
+%   by the difference equation:
+%
+%     y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
+%                      - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
+%
+%
+%   After filtering in the forward direction, the filtered sequence is then 
+%   reversed and run back through the filter; Y is the time reverse of the 
+%   output of the second filtering operation.  The result has precisely zero 
+%   phase distortion and magnitude modified by the square of the filter's 
+%   magnitude response.  Care is taken to minimize startup and ending 
+%   transients by matching initial conditions.
+%
+%   The length of the input x must be more than three times
+%   the filter order, defined as max(length(b)-1,length(a)-1).
+%
+%   Note that FILT2 should not be used with differentiator and Hilbert FIR
+%   filters, since the operation of these filters depends heavily on their
+%   phase response.
+%
+%   See also FILTER.
+
+%   References: 
+%     [1] Sanjit K. Mitra, Digital Signal Processing, 2nd ed., McGraw-Hill, 2001
+%     [2] Fredrik Gustafsson, Determining the initial states in forward-backward 
+%         filtering, IEEE Transactions on Signal Processing, pp. 988--992, April 1996, 
+%         Volume 44, Issue 4
+
+%   Author(s): L. Shure, 5-17-88
+%   revised by T. Krauss, 1-21-94
+%   Initial Conditions: Fredrik Gustafsson
+%   Copyright 1988-2004 The MathWorks, Inc.
+%   $Revision: 1.7.4.2 $  $Date: 2004/12/26 22:15:49 $
+
+    narginchk(3,3)
+    if (isempty(b) || isempty(a) || isempty(x))
+        y = [];
+        return
+    end
+
+    [m,n] = size(x);
+    if (n>1) && (m>1)
+        y = x;
+        for i=1:n  % loop over columns
+            y(:,i) = filt2(b,a,x(:,i));
+        end
+        return
+        % error('Only works for vector input.')
+    end
+    if m==1
+        x = x(:);   % convert row to column
+    end
+    len = size(x,1);   % length of input
+    b = b(:).';
+    a = a(:).';
+    nb = length(b);
+    na = length(a);
+    nfilt = max(nb,na);
+
+    nfact = 3*(nfilt-1);  % length of edge transients
+
+    if (len<=nfact),    % input data too short!
+        error('Data must have length more than 3 times filter order.');
+    end
+
+% set up filter's initial conditions to remove dc offset problems at the 
+% beginning and end of the sequence
+    if nb < nfilt, b(nfilt)=0; end   % zero-pad if necessary
+    if na < nfilt, a(nfilt)=0; end
+% use sparse matrix to solve system of linear equations for initial conditions
+% zi are the steady-state states of the filter b(z)/a(z) in the state-space 
+% implementation of the 'filter' command.
+    rows = [1:nfilt-1  2:nfilt-1  1:nfilt-2];
+    cols = [ones(1,nfilt-1) 2:nfilt-1  2:nfilt-1];
+    data = [1+a(2) a(3:nfilt) ones(1,nfilt-2)  -ones(1,nfilt-2)];
+    sp = sparse(rows,cols,data);
+    zi = sp \ ( b(2:nfilt).' - a(2:nfilt).'*b(1) );
+% non-sparse:
+% zi = ( eye(nfilt-1) - [-a(2:nfilt).' [eye(nfilt-2); zeros(1,nfilt-2)]] ) \ ...
+%      ( b(2:nfilt).' - a(2:nfilt).'*b(1) );
+
+% Extrapolate beginning and end of data sequence using a "reflection
+% method".  Slopes of original and extrapolated sequences match at
+% the end points.
+% This reduces end effects.
+    y = [2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)];
+
+% filter, reverse data, filter again, and reverse data again
+    y = filter(b,a,y,zi*y(1));
+    y = y(length(y):-1:1);
+    y = filter(b,a,y,zi*y(1));
+    y = y(length(y):-1:1);
+
+% remove extrapolated pieces of y
+    y([1:nfact len+nfact+(1:nfact)]) = [];
+
+    if m == 1
+        y = y.';   % convert back to row if necessary
+    end
+
+
+%**************************************************************
+%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+function [b,a] = fir1(N,Wn,varargin)
+%FIR1   FIR filter design using the window method.
+%   B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter
+%   and returns the filter coefficients in length N+1 vector B.
+%   The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0 
+%   corresponding to half the sample rate.  The filter B is real and
+%   has linear phase.  The normalized gain of the filter at Wn is
+%   -6 dB.
+%
+%   B = FIR1(N,Wn,'high') designs an N'th order highpass filter.
+%   You can also use B = FIR1(N,Wn,'low') to design a lowpass filter.
+%
+%   If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an
+%   order N bandpass filter with passband  W1 < W < W2. You can
+%   also specify B = FIR1(N,Wn,'bandpass').  If Wn = [W1 W2], 
+%   B = FIR1(N,Wn,'stop') will design a bandstop filter.
+%
+%   If Wn is a multi-element vector, 
+%          Wn = [W1 W2 W3 W4 W5 ... WN],
+%   FIR1 returns an order N multiband filter with bands
+%    0 < W < W1, W1 < W < W2, ..., WN < W < 1.
+%   B = FIR1(N,Wn,'DC-1') makes the first band a passband.
+%   B = FIR1(N,Wn,'DC-0') makes the first band a stopband.
+%
+%   B = FIR1(N,Wn,WIN) designs an N-th order FIR filter using 
+%   the N+1 length vector WIN to window the impulse response.
+%   If empty or omitted, FIR1 uses a Hamming window of length N+1.
+%   For a complete list of available windows, see the help for the
+%   WINDOW function. KAISER and CHEBWIN can be specified with an 
+%   optional trailing argument.  For example, B = FIR1(N,Wn,kaiser(N+1,4)) 
+%   uses a Kaiser window with beta=4. B = FIR1(N,Wn,'high',chebwin(N+1,R)) 
+%   uses a Chebyshev window with R decibels of relative sidelobe 
+%   attenuation.
+%
+%   For filters with a gain other than zero at Fs/2, e.g., highpass
+%   and bandstop filters, N must be even.  Otherwise, N will be
+%   incremented by one.  In this case the window length should be
+%   specified as N+2.
+%   
+%   By default, the filter is scaled so the center of the first pass band 
+%   has magnitude exactly one after windowing. Use a trailing 'noscale' 
+%   argument to prevent this scaling, e.g. B = FIR1(N,Wn,'noscale'), 
+%   B = FIR1(N,Wn,'high','noscale'), B = FIR1(N,Wn,wind,'noscale').  You
+%   can also specify the scaling explicitly, e.g. FIR1(N,Wn,'scale'), etc.
+%
+%   See also KAISERORD, FIRCLS1, FIR2, FIRLS, FIRCLS, CFIRPM,
+%            FIRPM, FREQZ, FILTER, WINDOW.
+
+%   FIR1 is an M-file implementation of program 5.2 in the IEEE
+%   Programs for Digital Signal Processing tape. 
+
+%   Author(s): L. Shure
+%              L. Shure, 4-5-90, revised
+%              T. Krauss, 3-5-96, revised
+%   Copyright 1988-2004 The MathWorks, Inc.
+%   $Revision: 1.15.4.2 $  $Date: 2004/04/13 00:17:49 $
+
+%   Reference(s):
+%     [1] "Programs for Digital Signal Processing", IEEE Press
+%         John Wiley & Sons, 1979, pg. 5.2-1.
+
+narginchk(2,5);
+
+% Parse optional input arguments
+[Ftype,Wind,SCALING,msg] = parseoptargs(Wn,varargin{:});
+error(msg);
+
+% Compute the frequency vector
+[nbands,ff,Ftype,msg] = desiredfreq(Wn,Ftype);
+error(msg);
+
+% Compute the magnitude vector
+[aa,First_Band] = desiredmag(Ftype,nbands);
+
+% Check for appropriate filter order, increase when necessary
+[N,msg1,msg2] = firchk(N,ff(end),aa);
+error(msg1);
+warning(msg2);
+
+% Work with filter length (= order + 1)
+L = N + 1;
+
+% Check for valid window, or assign default if empty
+[Wind,msg] = chkwindow(Wind,L);
+error(msg);
+
+% Compute unwindowed impulse response
+hh = firls(L-1,ff,aa);
+
+% Window impulse response to get the filter
+b = hh.*Wind(:)'; 
+a = 1;
+
+if SCALING,   
+    % Scale so that passband is approx 1
+    b = scalefilter(b,First_Band,ff,L);   
+end
+
+%-------------------------------------------------------------------------
+function [Ftype,Wind,SCALING,msg] = parseoptargs(Wn,varargin)
+%PARSEOPTARGS   Parse optional input arguments.
+
+% Up to 3 optional input arguments, always in this order:
+%   1 - Filter type flag, can be 'low','high','bandpass','stop','DC-0','DC-1'
+%   2 - Window vector
+%   3 - 'noscale' flag
+
+% Initialize output args.
+msg = '';
+SCALING = [];
+
+[Ftype,Wind,Scale] = assignoptargs(Wn,varargin{:});
+
+[Ftype,Wind,Scale,msg] = validateargs(Wn,Ftype,Wind,Scale);
+if ~isempty(msg),
+    return
+end
+
+switch lower(Scale),
+case 'noscale';
+    SCALING = 0;
+case 'scale';
+    SCALING = 1;
+end
+
+%--------------------------------------------------------------------------
+function [Ftype,Wind,Scale] = assignoptargs(Wn,varargin)
+%ASSIGNOPTARGS  Assign optional input arguments to the appropriate variables.
+
+% default optional parameter values:
+Wind = [];
+Scale = 'scale';
+Ftype = defaultftype(Wn);
+
+
+switch length(varargin)
+case 1
+    if ischar(varargin{1}) && (length(varargin{1})>0),
+        s = upper(varargin{1});
+        switch upper(s)
+        case {'SCALE','NOSCALE'}
+            Scale = s;
+        otherwise
+            Ftype = s;
+        end
+    else
+        Wind = varargin{1};
+    end
+case 2
+    if ischar(varargin{1})
+        Ftype = varargin{1};
+    else
+        Wind = varargin{1};
+    end
+    if ischar(varargin{2})
+        Scale = varargin{2};
+    else
+        Wind = varargin{2};
+    end
+case 3
+    Ftype = varargin{1};
+    Wind = varargin{2};
+    Scale = varargin{3};
+end
+
+%--------------------------------------------------------------------------
+function [Ftype,Wind,Scale,msg] = validateargs(Wn,Ftype,Wind,Scale)
+%VALIDATEARGS  Test if arguments are valid.
+
+msg = '';
+
+% Assign a default Ftype when an empty is given. Backwards compatibility
+if isempty(Ftype),
+    Ftype = defaultftype(Wn);
+end
+
+
+Ftypeopts = {'LOW','HIGH','BANDPASS','STOP','DC-0','DC-1'};
+Scaleopts = {'NOSCALE','SCALE'};
+
+
+indx = strmatch(upper(Ftype),Ftypeopts);
+if isempty(indx),
+    msg = 'Unrecognized or ambiguous filter type specified.';
+    return
+else
+    Ftype = Ftypeopts{indx};
+end
+
+scaleindx = strmatch(upper(Scale),Scaleopts);
+if isempty(scaleindx),
+    msg = 'Scaling option must be ''noscale'' or ''scale''.';
+    return
+else
+    Scale = Scaleopts{scaleindx};
+end
+
+if ~any(size(Wind) <= 1),
+    msg = 'The window specified must be a vector';
+    return
+else
+    Wind = Wind(:).'; % Make it a row vector
+end
+
+%--------------------------------------------------------------------------
+function [nbands,ff,Ftype,msg] = desiredfreq(Wn,Ftype)
+%DESIREDFREQ  Compute the vector of frequencies to pass to FIRLS.
+%
+%   Inputs:
+%           Wn    - vector of cutoff frequencies.
+%           Ftype - string with desired response ('low','high',...)
+%
+%   Outputs:
+%           nbands - number of frequency bands.
+%           ff     - vector of frequencies to pass to FIRLS.
+%           Ftype  - converted filter type (if it's necessary to convert)
+
+% Initialize output args.
+nbands = [];
+ff     = [];
+msg    = '';
+
+
+if  any( Wn<0 | Wn>1 ),
+   msg = 'Frequencies must fall in range between 0 and 1.';
+   return
+end
+if  any(diff(Wn)<0),
+   msg = 'Frequencies must be increasing';
+   return
+end
+
+Wn = Wn(:)';
+
+nbands = length(Wn) + 1;
+
+if (nbands > 2) && strcmp(lower(Ftype),'bandpass'),
+    Ftype = 'DC-0';  % make sure default 3 band filter is bandpass
+end
+
+ff = [0,Wn(1:nbands-1); Wn(1:nbands-1),1];
+
+ff = ff(:);
+
+%-------------------------------------------------------------------------
+function [aa,First_Band] = desiredmag(Ftype,nbands)
+%DESIREDMAG  Compute the magnitude vector to pass to FIRLS.
+
+First_Band = isempty(findstr('DC-0',Ftype)) && isempty(findstr('HIGH',Ftype));
+mags = rem( First_Band + (0:nbands-1), 2);
+aa = [mags(:)'; mags(:)'];
+
+aa = aa(:);
+%--------------------------------------------------------------------------
+function [Wind,msg] = chkwindow(Wind,L)
+%CHKWINDOW   Check if specified window is valid, assign default if empty.
+
+msg = '';
+
+if isempty(Wind),
+   % Replace the following with the default window of your choice.
+   Wind = hamming(L);
+end
+    
+if length(Wind) ~= L
+   msg = 'The window length must be the same as the filter length.';
+end
+%
+% to use Kaiser window, beta must be supplied
+% att = 60; % dB of attenuation desired in sidelobe
+% beta = 0.1102*(att-8.7);
+% wind = kaiser(L,beta);
+
+%---------------------------------------------------------------------------
+function b = scalefilter(b,First_Band,ff,L)
+%SCALEFILTER   Scale fitler to have passband approx. equal to one.
+
+if First_Band
+    b = b / sum(b);  % unity gain at DC
+else
+    if ff(4)==1
+        % unity gain at Fs/2
+        f0 = 1;
+    else
+        % unity gain at center of first passband
+        f0 = mean(ff(3:4));
+    end
+    b = b / abs( exp(-j*2*pi*(0:L-1)*(f0/2))*(b.') );
+end
+
+%----------------------------------------------------------------------------
+function Ftype = defaultftype(Wn)
+%DEFAULTFTYPE  Assign default filter type depending on number of bands.
+
+if length(Wn) == 1,
+    Ftype = 'low';
+elseif length(Wn) == 2,
+    Ftype = 'bandpass';
+elseif length(Wn) >= 3,
+    Ftype = 'dc-0';
+end
+
+%**************************************************************
+%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+function [n,msg1,msg2] = firchk(n,Fend,a,exception)
+%FIRCHK   Check if specified filter order is valid.
+%   FIRCHK(N,Fend,A) checks if the specified order N is valid given the
+%   final frequency point Fend and the desired magnitude response vector A.
+%   Type 2 linear phase FIR filters (symmetric, odd order) must have a
+%   desired magnitude response vector that ends in zero if Fend = 1.  This
+%   is because type 2 filters necessarily have a zero at w = pi.
+%
+%   If the order is not valid, a warning is given and the order
+%   of the filter is incremented by one.
+%
+%   If A is a scalar (as when called from fircls1), A = 0 is
+%   interpreted as lowpass and A = 1 is interpreted as highpass.
+%
+%   FIRCHK(N,Fend,A,EXCEPTION) will not warn or increase the order
+%   if EXCEPTION = 1.  Examples of EXCEPTIONS are type 4 filters
+%   (such as differentiators or hilbert transformers) or non-linear
+%   phase filters (such as minimum and maximum phase filters).
+
+%   Author : R. Losada
+%   Copyright 1988-2004 The MathWorks, Inc.
+%   $Revision: 1.7.4.3 $  $Date: 2004/04/13 00:18:44 $
+
+narginchk(3,4);
+
+if nargin == 3,
+    exception = false;
+end
+
+msg1 = '';
+msg2 = '';
+oddord = false; % Flag, initially we assume even order
+
+if isempty(n) || length(n) > 1 || ~isnumeric(n) || ~isreal(n) || n~=round(n) || n<=0,
+    msg1 = 'Filter order must be a real, positive integer.';
+    return
+end
+
+if rem(n,2) == 1,
+    oddord = true; % Overwrite flag
+end
+ 
+if (a(end) ~= 0) && Fend == 1 && oddord && ~exception,
+    str = ['Odd order symmetric FIR filters must have a gain of zero \n'...
+     'at the Nyquist frequency. The order is being increased by one.'];
+    msg2 = sprintf(str);
+    n = n+1;
+end
+    
+%**************************************************************
+%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+function [h,a]=firls(N,F,M,W,ftype);
+% FIRLS Linear-phase FIR filter design using least-squares error minimization.
+%   B=FIRLS(N,F,A) returns a length N+1 linear phase (real, symmetric
+%   coefficients) FIR filter which has the best approximation to the
+%   desired frequency response described by F and A in the least squares
+%   sense. F is a vector of frequency band edges in pairs, in ascending
+%   order between 0 and 1. 1 corresponds to the Nyquist frequency or half
+%   the sampling frequency. A is a real vector the same size as F
+%   which specifies the desired amplitude of the frequency response of the
+%   resultant filter B. The desired response is the line connecting the
+%   points (F(k),A(k)) and (F(k+1),A(k+1)) for odd k; FIRLS treats the
+%   bands between F(k+1) and F(k+2) for odd k as "transition bands" or
+%   "don't care" regions. Thus the desired amplitude is piecewise linear
+%   with transition bands.  The integrated squared error is minimized.
+%
+%   For filters with a gain other than zero at Fs/2, e.g., highpass
+%   and bandstop filters, N must be even.  Otherwise, N will be
+%   incremented by one. Alternatively, you can use a trailing 'h' flag to
+%   design a type 4 linear phase filter and avoid incrementing N.
+%
+%   B=FIRLS(N,F,A,W) uses the weights in W to weight the error. W has one
+%   entry per band (so it is half the length of F and A) which tells
+%   FIRLS how much emphasis to put on minimizing the integral squared error
+%   in each band relative to the other bands.
+%
+%   B=FIRLS(N,F,A,'Hilbert') and B=FIRLS(N,F,A,W,'Hilbert') design filters
+%   that have odd symmetry, that is, B(k) = -B(N+2-k) for k = 1, ..., N+1.
+%   A special case is a Hilbert transformer which has an approx. amplitude
+%   of 1 across the entire band, e.g. B=FIRLS(30,[.1 .9],[1 1],'Hilbert').
+%
+%   B=FIRLS(N,F,A,'differentiator') and B=FIRLS(N,F,A,W,'differentiator')
+%   also design filters with odd symmetry, but with a special weighting
+%   scheme for non-zero amplitude bands. The weight is assumed to be equal
+%   to the inverse of frequency, squared, times the weight W. Thus the
+%   filter has a much better fit at low frequency than at high frequency.
+%   This designs FIR differentiators.
+%
+%   % Example of a length 31 lowpass filter:
+%   	h=firls(30,[0 .1 .2 .5]*2,[1 1 0 0]);
+%
+%   % Example of a low-pass differentiator:
+%   	h=firls(44,[0 .3 .4 1],[0 .2 0 0],'differentiator');
+%
+%   % Example of a type 4 highpass filter:
+%       h=firls(25,[0 .4 .5 1],[0 0 1 1],'h');
+%
+%   See also FIRPM, FIR1, FIR2, FREQZ and FILTER.
+
+%       Author(s): T. Krauss
+%   History: 10-18-91, original version
+%            3-30-93, updated
+%            9-1-95, optimize adjacent band case
+%   Copyright 1988-2004 The MathWorks, Inc.
+%   $Revision: 1.11.4.2 $  $Date: 2004/04/13 00:17:54 $
+
+% check number of arguments, set up defaults.
+narginchk(3,5);
+
+if (max(F)>1) || (min(F)<0)
+    error('Frequencies in F must be in range [0,1].')
+end
+if (rem(length(F),2)~=0)
+    error('F must have even length.');
+end
+if (length(F) ~= length(M))
+    error('F and A must be equal lengths.');
+end
+if (nargin==3),
+    W = ones(length(F)/2,1);
+    ftype = '';
+end
+if (nargin==4),
+    if isstr(W),
+        ftype = W; W = ones(length(F)/2,1);
+    else
+        ftype = '';
+    end
+end
+if (nargin==5),
+    if isempty(W),
+        W = ones(length(F)/2,1);
+    end
+end
+if isempty(ftype)
+    ftype = 0;  differ = 0;
+else
+    ftype = lower(ftype);
+    if strcmpi(ftype,'h') || strcmpi(ftype,'hilbert')
+        ftype = 1;  differ = 0;
+    elseif strcmpi(ftype,'d') || strcmpi(ftype,'differentiator')
+        ftype = 1;  differ = 1;
+    else
+        error('Requires symmetry to be ''Hilbert'' or ''differentiator''.')
+    end
+end
+
+% Check for valid filter length
+[N,msg1,msg2] = firchk(N,F(end),M,ftype);
+error(msg1);
+
+if ~isempty(msg2),
+    msg2 = sprintf([msg2,'\r',...
+        '\nAlternatively, you can pass a trailing ''h'' argument,\r',...
+        'as in firls(N,F,A,W,''h''), to design a type 4 linear phase filter.']);
+end
+warning(msg2);
+
+
+N = N+1;                   % filter length
+F=F(:)/2;  M=M(:);  W=sqrt(W(:));  % make these guys columns
+dF = diff(F);
+
+if (length(F) ~= length(W)*2)
+    error('There should be one weight per band.');
+end;
+if any(dF<0),
+    error('Frequencies in F must be nondecreasing.')
+end
+if all(dF(2:2:length(dF)-1)==0)
+    fullband = 1;
+else
+    fullband = 0;
+end
+if all((W-W(1))==0)
+    constant_weights = 1;
+else
+    constant_weights = 0;
+end
+
+L=(N-1)/2;
+
+Nodd = rem(N,2);
+
+if (ftype == 0),  % Type I and Type II linear phase FIR
+    % basis vectors are cos(2*pi*m*f) (see m below)
+    if ~Nodd
+        m=(0:L)+.5;   % type II
+    else
+        m=(0:L);      % type I
+    end
+    k=m';
+    need_matrix = (~fullband) || (~constant_weights);
+    if need_matrix
+        I1=k(:,ones(size(m)))+m(ones(size(k)),:);    % entries are m + k
+        I2=k(:,ones(size(m)))-m(ones(size(k)),:);    % entries are m - k
+        G=zeros(size(I1));
+    end
+
+    if Nodd
+        k=k(2:length(k));
+        b0=0;       %  first entry must be handled separately (where k(1)=0)
+    end;
+    b=zeros(size(k));
+    for s=1:2:length(F),
+        m=(M(s+1)-M(s))/(F(s+1)-F(s));    %  slope
+        b1=M(s)-m*F(s);                   %  y-intercept
+        if Nodd
+            b0 = b0 + (b1*(F(s+1)-F(s)) + m/2*(F(s+1)*F(s+1)-F(s)*F(s)))...
+                * abs(W((s+1)/2)^2) ;
+        end
+        b = b+(m/(4*pi*pi)*(cos(2*pi*k*F(s+1))-cos(2*pi*k*F(s)))./(k.*k))...
+            * abs(W((s+1)/2)^2);
+        b = b + (F(s+1)*(m*F(s+1)+b1)*sinc(2*k*F(s+1)) ...
+            - F(s)*(m*F(s)+b1)*sinc(2*k*F(s))) ...
+            * abs(W((s+1)/2)^2);
+        if need_matrix
+            G = G + (.5*F(s+1)*(sinc(2*I1*F(s+1))+sinc(2*I2*F(s+1))) ...
+                - .5*F(s)*(sinc(2*I1*F(s))+sinc(2*I2*F(s))) ) ...
+                * abs(W((s+1)/2)^2);
+        end
+    end;
+    if Nodd
+        b=[b0; b];
+    end;
+
+    if need_matrix
+        a=G\b;
+    else
+        a=(W(1)^2)*4*b;
+        if Nodd
+            a(1) = a(1)/2;
+        end
+    end
+    if Nodd
+        h=[a(L+1:-1:2)/2; a(1); a(2:L+1)/2].';
+    else
+        h=.5*[flipud(a); a].';
+    end;
+elseif (ftype == 1),  % Type III and Type IV linear phase FIR
+    %  basis vectors are sin(2*pi*m*f) (see m below)
+    if (differ),      % weight non-zero bands with 1/f^2
+        do_weight = ( abs(M(1:2:length(M))) +  abs(M(2:2:length(M))) ) > 0;
+    else
+        do_weight = zeros(size(F));
+    end
+
+    if Nodd
+        m=(1:L);      % type III
+    else
+        m=(0:L)+.5;   % type IV
+    end;
+    k=m';
+    b=zeros(size(k));
+
+    need_matrix = (~fullband) || (any(do_weight)) || (~constant_weights);
+    if need_matrix
+        I1=k(:,ones(size(m)))+m(ones(size(k)),:);    % entries are m + k
+        I2=k(:,ones(size(m)))-m(ones(size(k)),:);    % entries are m - k
+        G=zeros(size(I1));
+    end
+
+    i = sqrt(-1);
+    for s=1:2:length(F),
+        if (do_weight((s+1)/2)),      % weight bands with 1/f^2
+            if F(s) == 0, F(s) = 1e-5; end     % avoid singularities
+            m=(M(s+1)-M(s))/(F(s+1)-F(s));
+            b1=M(s)-m*F(s);
+            snint1 = sineint(2*pi*k*F(s+1)) - sineint(2*pi*k*F(s));
+            %snint1 = (-1/2/i)*(expint(i*2*pi*k*F(s+1)) ...
+            %    -expint(-i*2*pi*k*F(s+1)) -expint(i*2*pi*k*F(s)) ...
+            %    +expint(-i*2*pi*k*F(s)) );
+            % csint1 = cosint(2*pi*k*F(s+1)) - cosint(2*pi*k*F(s)) ;
+            csint1 = (-1/2)*(expint(i*2*pi*k*F(s+1))+expint(-i*2*pi*k*F(s+1))...
+                -expint(i*2*pi*k*F(s))  -expint(-i*2*pi*k*F(s)) );
+            b=b + ( m*snint1 ...
+                + b1*2*pi*k.*( -sinc(2*k*F(s+1)) + sinc(2*k*F(s)) + csint1 ))...
+                * abs(W((s+1)/2)^2);
+            snint1 = sineint(2*pi*F(s+1)*(-I2));
+            snint2 = sineint(2*pi*F(s+1)*I1);
+            snint3 = sineint(2*pi*F(s)*(-I2));
+            snint4 = sineint(2*pi*F(s)*I1);
+            G = G - ( ( -1/2*( cos(2*pi*F(s+1)*(-I2))/F(s+1)  ...
+                - 2*snint1*pi.*I2 ...
+                - cos(2*pi*F(s+1)*I1)/F(s+1) ...
+                - 2*snint2*pi.*I1 )) ...
+                - ( -1/2*( cos(2*pi*F(s)*(-I2))/F(s)  ...
+                - 2*snint3*pi.*I2 ...
+                - cos(2*pi*F(s)*I1)/F(s) ...
+                - 2*snint4*pi.*I1) ) ) ...
+                * abs(W((s+1)/2)^2);
+        else      % use usual weights
+            m=(M(s+1)-M(s))/(F(s+1)-F(s));
+            b1=M(s)-m*F(s);
+            b=b+(m/(4*pi*pi)*(sin(2*pi*k*F(s+1))-sin(2*pi*k*F(s)))./(k.*k))...
+                * abs(W((s+1)/2)^2) ;
+            b = b + (((m*F(s)+b1)*cos(2*pi*k*F(s)) - ...
+                (m*F(s+1)+b1)*cos(2*pi*k*F(s+1)))./(2*pi*k)) ...
+                * abs(W((s+1)/2)^2) ;
+            if need_matrix
+                G = G + (.5*F(s+1)*(sinc(2*I1*F(s+1))-sinc(2*I2*F(s+1))) ...
+                    - .5*F(s)*(sinc(2*I1*F(s))-sinc(2*I2*F(s)))) * ...
+                    abs(W((s+1)/2)^2);
+            end
+        end;
+    end
+
+    if need_matrix
+        a=G\b;
+    else
+        a=-4*b*(W(1)^2);
+    end
+    if Nodd
+        h=.5*[flipud(a); 0; -a].';
+    else
+        h=.5*[flipud(a); -a].';
+    end
+    if differ, h=-h; end
+end
+
+if nargout > 1
+    a = 1;
+end
+
+%----------------------------------------------------------------------------
+function y = sineint(x)
+% SINEINT (a.k.a. SININT)   Numerical Sine Integral
+%   Used by FIRLS in the Signal Processing Toolbox.
+%   Untested for complex or imaginary inputs.
+%
+%   See also SININT in the Symbolic Toolbox.
+
+%   Was Revision: 1.5, Date: 1996/03/15 20:55:51
+
+i1 = find(real(x)<0);   % this equation is not valid if x is in the
+% left-hand plane of the complex plane.
+% use relation Si(-z) = -Si(z) in this case (Eq 5.2.19, Abramowitz
+%  & Stegun).
+x(i1) = -x(i1);
+y = zeros(size(x));
+ind = find(x);
+% equation 5.2.21 Abramowitz & Stegun
+%  y(ind) = (1/(2*i))*(expint(i*x(ind)) - expint(-i*x(ind))) + pi/2;
+y(ind) = imag(expint(i*x(ind))) + pi/2;
+y(i1) = -y(i1);
+
+%**************************************************************
+%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+function y=sinc(x)
+%SINC Sin(pi*x)/(pi*x) function.
+%   SINC(X) returns a matrix whose elements are the sinc of the elements 
+%   of X, i.e.
+%        y = sin(pi*x)/(pi*x)    if x ~= 0
+%          = 1                   if x == 0
+%   where x is an element of the input matrix and y is the resultant
+%   output element.
+%
+%   See also SQUARE, SIN, COS, CHIRP, DIRIC, GAUSPULS, PULSTRAN, RECTPULS,
+%   and TRIPULS.
+
+%   Author(s): T. Krauss, 1-14-93
+%   Copyright 1988-2002 The MathWorks, Inc.
+%       $Revision: 1.7 $  $Date: 2002/04/15 01:13:58 $
+
+y=ones(size(x));
+i=find(x);
+y(i)=sin(pi*x(i))./(pi*x(i));
+
+%**************************************************************
+%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+function w = hamming(varargin)
+%HAMMING   Hamming window.
+%   HAMMING(N) returns the N-point symmetric Hamming window in a column vector.
+% 
+%   HAMMING(N,SFLAG) generates the N-point Hamming window using SFLAG window
+%   sampling. SFLAG may be either 'symmetric' or 'periodic'. By default, a 
+%   symmetric window is returned. 
+%
+%   See also BLACKMAN, HANN, WINDOW.
+
+%   Copyright 1988-2002 The MathWorks, Inc.
+%   $Revision: 1.14 $  $Date: 2002/11/21 15:46:43 $
+
+% Check number of inputs
+narginchk(1,2);
+
+[w,msg] = gencoswin('hamming',varargin{:});
+error(msg);
+
+
+% [EOF] hamming.m
+
+%**************************************************************
+%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+function [w,msg] = gencoswin(varargin)
+%GENCOSWIN   Returns one of the generalized cosine windows.
+%   GENCOSWIN returns the generalized cosine window specified by the 
+%   first string argument. Its inputs can be
+%     Window name    - a string, any of 'hamming', 'hann', 'blackman'.
+%     N              - length of the window desired.
+%     Sampling flag  - optional string, one of 'symmetric', 'periodic'. 
+
+%   Copyright 1988-2002 The MathWorks, Inc.
+%   $Revision: 1.7 $  $Date: 2002/04/15 01:09:14 $ 
+
+% Parse the inputs
+window = varargin{1};
+n = varargin{2};
+msg = '';
+
+% Check for trivial orders
+[n,w,trivialwin] = check_order(n);
+if trivialwin, return, end;
+
+% Select the sampling option
+if nargin == 2, % no sampling flag specified, use default. 
+    sflag = 'symmetric';
+else
+    sflag = lower(varargin{3});
+end
+
+% Allow partial strings for sampling options
+allsflags = {'symmetric','periodic'};
+sflagindex = strmatch(sflag, allsflags);
+if length(sflagindex)~=1         % catch 0 or 2 matches
+    msg = 'Sampling flag must be either ''symmetric'' or ''periodic''.';
+    return;
+else	
+    sflag = allsflags{sflagindex};
+end
+
+% Evaluate the window
+switch sflag
+case 'periodic'
+    w = sym_window(n+1,window);
+    w(end) = [];
+case 'symmetric'
+    w = sym_window(n,window);
+end
+
+%---------------------------------------------------------------------
+function w = sym_window(n,window)
+%SYM_WINDOW   Symmetric generalized cosine window.
+%   SYM_WINDOW Returns an exactly symmetric N point generalized cosine 
+%   window by evaluating the first half and then flipping the same samples
+%   over the other half.
+
+if ~rem(n,2)
+    % Even length window
+    half = n/2;
+    w = calc_window(half,n,window);
+    w = [w; w(end:-1:1)];
+else
+    % Odd length window
+    half = (n+1)/2;
+    w = calc_window(half,n,window);
+    w = [w; w(end-1:-1:1)];
+end
+
+%---------------------------------------------------------------------
+function w = calc_window(m,n,window)
+%CALC_WINDOW   Calculate the generalized cosine window samples.
+%   CALC_WINDOW Calculates and returns the first M points of an N point
+%   generalized cosine window determined by the 'window' string.
+
+% For the hamming and blackman windows we force rounding in order to achieve
+% better numerical properties.  For example, the end points of the hamming 
+% window should be exactly 0.08.
+
+switch window
+case 'hann'
+    % Hann window
+    %    w = 0.5 * (1 - cos(2*pi*(0:m-1)'/(n-1))); 
+    a0 = 0.5;
+    a1 = 0.5;
+    a2 = 0;
+    a3 = 0;
+    a4 = 0;
+case 'hamming'
+    % Hamming window
+    %    w = (54 - 46*cos(2*pi*(0:m-1)'/(n-1)))/100;
+    a0 = 0.54;
+    a1 = 0.46;
+    a2 = 0;
+    a3 = 0;
+    a4 = 0;
+case 'blackman'
+    % Blackman window
+    %    w = (42 - 50*cos(2*pi*(0:m-1)/(n-1)) + 8*cos(4*pi*(0:m-1)/(n-1)))'/100;
+    a0 = 0.42;
+    a1 = 0.5;
+    a2 = 0.08;
+    a3 = 0;
+    a4 = 0;
+case 'flattopwin'
+    % Flattop window
+    % Original coefficients as defined in the reference (see flattopwin.m);
+    % a0 = 1;
+    % a1 = 1.93;
+    % a2 = 1.29;
+    % a3 = 0.388;
+    % a4 = 0.032;
+    %
+    % Scaled by (a0+a1+a2+a3+a4)
+    a0 = 0.2156;
+    a1 = 0.4160;
+    a2 = 0.2781;
+    a3 = 0.0836;
+    a4 = 0.0069;
+end
+
+x = (0:m-1)'/(n-1);
+w = a0 - a1*cos(2*pi*x) + a2*cos(4*pi*x) - a3*cos(6*pi*x) + a4*cos(8*pi*x);
+
+% [EOF] gencoswin.m
+
+%----------------------------------------------------------------------------
+
+function [n_out, w, trivalwin] = check_order(n_in)
+%CHECK_ORDER Checks the order passed to the window functions.
+% [N,W,TRIVALWIN] = CHECK_ORDER(N_ESTIMATE) will round N_ESTIMATE to the
+% nearest integer if it is not alreay an integer. In special cases (N is [],
+% 0, or 1), TRIVALWIN will be set to flag that W has been modified.
+
+%   Copyright 1988-2002 The MathWorks, Inc.
+%   $Revision: 1.6 $  $Date: 2002/04/15 01:07:36 $
+
+w = [];
+trivalwin = 0;
+
+% Special case of negative orders:
+if n_in < 0,
+   error('Order cannot be less than zero.');
+end
+
+% Check if order is already an integer or empty
+% If not, round to nearest integer.
+if isempty(n_in) | n_in == floor(n_in),
+   n_out = n_in;
+else
+   n_out = round(n_in);
+   warning('Rounding order to nearest integer.');
+end
+
+% Special cases:
+if isempty(n_out) | n_out == 0,
+   w = zeros(0,1);               % Empty matrix: 0-by-1
+   trivalwin = 1; 
+elseif n_out == 1,
+   w = 1;
+   trivalwin = 1;   
+end
diff --git a/filter/naninterp.m b/filter/naninterp.m
new file mode 100644
index 0000000..8f5eebf
--- /dev/null
+++ b/filter/naninterp.m
@@ -0,0 +1,5 @@
+function X = naninterp(X)
+% Interpolate over NaNs
+% See INTERP1 for more info
+X(isnan(X)) = interp1(find(~isnan(X)), X(~isnan(X)), find(isnan(X)),'cubic');
+return
\ No newline at end of file
diff --git a/filter/nanmedfilt2.m b/filter/nanmedfilt2.m
new file mode 100644
index 0000000..512be4c
--- /dev/null
+++ b/filter/nanmedfilt2.m
@@ -0,0 +1,44 @@
+function [M] = nanmedfilt2(A, sz)
+%This MATLAB function performs median filtering of the matrix A in two dimensions
+%while *ignoring* NaNs (based on discussions here
+%http://www.mathworks.com/matlabcentral/newsreader/view_thread/251787)
+
+%A is a 2d matrix, sz is tile size, M is filtered A.
+
+if nargin<2
+    sz = 5;
+end
+if length(sz)==1
+    sz = [sz sz];
+end
+
+if any(mod(sz,2)==0)
+    error('kernel size SZ must be odd)')
+end
+margin=(sz-1)/2;
+AA = nan(size(A)+2*margin);
+AA(1+margin(1):end-margin(1),1+margin(2):end-margin(2))=A;
+[iB,jB]=ndgrid(1:sz(1),1:sz(2));
+is=sub2ind(size(AA),iB,jB);
+[iA, jA]=ndgrid(1:size(A,1),1:size(A,2));
+iA=sub2ind(size(AA),iA,jA);
+idx=bsxfun(@plus,iA(:).',is(:)-1);
+
+B = sort(AA(idx),1);
+j=any(isnan(B),1);
+last = zeros(1,size(B,2))+size(B,1);
+[~, last(j)]=max(isnan(B(:,j)),[],1);
+last(j)=last(j)-1;
+
+M = nan(1,size(B,2));
+valid = find(last>0);
+mid = (1 + last)/2;
+i1 = floor(mid(valid));
+i2 = ceil(mid(valid));
+i1 = sub2ind(size(B),i1,valid);
+i2 = sub2ind(size(B),i2,valid);
+M(valid) = 0.5*(B(i1) + B(i2));
+M = reshape(M,size(A));
+
+end % medianna
+
diff --git a/filter/ndnanfilter.m b/filter/ndnanfilter.m
new file mode 100644
index 0000000..e7d3e96
--- /dev/null
+++ b/filter/ndnanfilter.m
@@ -0,0 +1,484 @@
+function [Y,W] = ndnanfilter(X,HWIN,F,DIM,WINOPT,PADOPT,WNAN)
+% NDNANFILTER   N-dimensional zero-phase digital filter, ignoring NaNs.
+%
+%   Syntax:
+%         Y = ndnanfilter(X,HWIN,F);
+%         Y = ndnanfilter(X,HWIN,F,DIM);
+%         Y = ndnanfilter(X,HWIN,F,DIM,WINOPT);
+%         Y = ndnanfilter(X,HWIN,F,DIM,WINOPT,PADOPT);
+%         Y = ndnanfilter(X,HWIN,F,DIM,WINOPT,PADOPT,WNAN);
+%     [Y,W] = ndnanfilter(...);
+%
+%   Input:
+%     X      - Data to be filtered with/without NaNs.
+%     HWIN   - Window function handle (or name) or numeric multidimensional
+%              window to be used (without NaNs). See WINDOW for details.   
+%              Default:   @rectwin  or 'rectwin' (moving average).
+%     F      - A vector specifying the semi-width of the window for each
+%              dimension. The final window's width will be 2*F+1. 
+%              Default: 3 (i.e. a 1-dimensional window of width 6).
+%     DIM    - If F is a single scalar, the window will be applied through
+%              this dimension; otherwise, this will be ignored. 
+%              Default: columns (or the first non-singleton dimension).
+%     WINOPT - Cell array specifying optional arguments for the window
+%              function HWIN (in addition to the width). 
+%              Default: {} (window's defaults).
+%     PADOPT - Cell array specifying the optional arguments for the
+%              PADARRAY MATLAB's function (in addition to the array X and
+%              the padsize: 2*F+1). If the function is not found, data is
+%              padded with zeros or the specified value: try {mean(X(:))}
+%              for example.
+%              Default: {'replicate'} (repeats border elements of X).
+%              Default: {0} (pads with zeros if PADARRAY not found).
+%     WNAN   - Integer indicating NaNs treatment and program behaviour!:
+%              0: Filters data and interpolates NaNs         (default). 
+%              1: Filters data but do not interpolates NaNs 
+%              2: "Do not filters data" but interpolates NaNs!
+%              See the NOTEs below
+%
+%   Output:
+%     Y      - Filtered X data (same size as X!).
+%     W      - N-dimensional window with central symmetry generated by a
+%              special subfunction called NDWIND. See the description below
+%              for details.
+%
+%   Description:
+%     This function applies a N-dimensional convolution of X with W, using
+%     the MATLAB's IMFILTER or CONVN function. One important aspect of the
+%     function is the generation of the N-dimensional window (W) from the
+%     specified function and width, which cannot be done with MATLAB's
+%     functions. Besides, unlike MATLAB's FILTER, FILTER2 and IMFILTER,
+%     NaNs elements are taken into account (ignored).
+%
+%     The N-dimensional window is generated from rotating the 1-dimensional
+%     output of the HWIN function, through each of the N-dimensions, and
+%     then shrinking it through each of its axes in order to fit the
+%     specified semi-widths (F). This is done in the included subfunction
+%     named NDWIND. In this way, the window has central symmetry and do not
+%     produce a phase shift on X data.
+%
+%     By default, the edges are padded with the values of X at the borders
+%     with the PADARRAY MATLAB's function. In this way, the edges are
+%     treated smoothly. When PADARRAY is not found, the program performs
+%     zero-padding.
+%
+%   Notes: 
+%     * The use of semi-widths F's is to force the generated window to be
+%       even and, therefore, the change of phase is null.  
+%     * The window function HWIN should output an even function, otherwise,
+%       it won't generate an error but the user should be aware that this
+%       program will consider only the last half of it.
+%     * The function window should return a monotonically decreasing
+%       result, this restriction is because I try to avoid the use of FZERO
+%       function, for example, to find the expanding/shrinking factors.
+%     * If the user has an already generated window, it can be used in HWIN
+%       instead of a function handle or name. 
+%     * Accepts empty value for any input. When X is empty, the program can
+%       be used as a N-dimensional window generator.
+%     * NaNs elements surrounded by no-NaNs elements (which will depend on
+%       window width) are the ones that will be interpolated. The others
+%       are leaved untouched.
+%     * When WNAN=2, the programs acts like an NAN-interpolat/GAP-filling,
+%       leaving untouched the no-NaNs elements but the filtering is
+%       perfomed anyway. I recomend the default behaviour (WNAN=0) in order
+%       to keep the filtered data in the workspace, and then use the code
+%       at the end of this function to get/remove the interpolated NaNs
+%       (see the example).
+%     * The program looks for the IMFILTER and PADARRAY functions from the
+%       Image Processing Toolbox. If not found, then CONVN is used instead
+%       (slower) and pads with zeros or the given value. In this latter
+%       case, if border elements are NaNs, the window won't work properly.
+%
+%   Example:
+%     FWIN = 'hamming';
+%     F = [13 8];
+%     N = 100;
+%     Pnoise = 0.30;
+%     PNaNs  = 0.20;
+%     X = peaks(N);                                     % original
+%     Y = X + ((rand(size(X))-0.5)*2)*max(X(:))*Pnoise; % add noise
+%     Y(round(1 + (N^2-1).*rand(N^2*PNaNs,1))) = NaN;   % add NaNs
+%     [Z0,W] = ndnanfilter(Y,FWIN,F);                   % filters
+%     Z1 = Z0; Z2 = Y; inan = isnan(Y);
+%     Z1(inan) = NaN;
+%     Z2(inan) = Z0(inan);  
+%     subplot(231), imagesc(X), clim = caxis; axis equal tight
+%                   title('Original data')
+%     subplot(232), imagesc(Y),  caxis(clim), axis equal tight 
+%                   title('Data + NOISE + NaNs')
+%     subplot(234), imagesc(Z0), caxis(clim), axis equal tight 
+%                   title('FILTERS + NaNs interpolation')
+%     subplot(235), imagesc(Z1), caxis(clim), axis equal tight 
+%                   title('FILTERS ignoring NaNs')
+%     subplot(236), imagesc(Z2), caxis(clim), axis equal tight 
+%                   title('GAP-filling with interpolated NaNs')
+%     subplot(233), imagesc(-F(1):F(1),-F(2):F(2),W), axis equal tight, 
+%                    title([upper(FWIN) ' 2D window']), view(2) 
+%
+%   See also: FILTER, FILTER2 and CONVN; WINDOW from the Signal Processing
+%   Toolbox; and FWIND1, FWIND2, FSPECIAL, IMFILTER and PADARRAY from the
+%   Image Processing Toolbox. 
+
+%   Copyright 2008 Carlos Adrian Vargas Aguilera
+%   $Revision: 1.2 $  $Date: 2008/06/30 18:00:00 $
+
+%   Written by
+%   M.S. Carlos Adrian Vargas Aguilera
+%   Physical Oceanography PhD candidate
+%   CICESE 
+%   Mexico, 2008
+%   nubeobscura@hotmail.com
+%
+%   Download from:
+%   http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objec
+%   tType=author&objectId=1093874
+
+%   1.0    Release (2008/06/23 10:30:00)
+%   1.1    Fixed Bug adding an extra dimension of unitary width. 
+%   1.2    Fixed Bug with ynan.
+
+% Use the IMFILTER function? (faster than CONVN):
+yimfilter = (exist('imfilter','file')==2);
+
+% Use the PADARRAY function (or zero padding): 
+ypadarray = (exist('padarray','file')==2);
+
+% Check inputs and sets defaults of principal arguments:
+if nargin<3 || nargin>7
+ error('Filtern:IncorrectNumberOfInputs',...
+  'At least three inputs are needed and less than 7.')
+end
+if isempty(HWIN)
+ HWIN = 'rectwin';
+end
+if isempty(F)
+ F = 3;
+end
+N = length(F);
+S = size(X);
+% Secondary arguments:
+if N && (nargin<4 || isempty(DIM))
+ DIM = find(S~=1,1);   %  DIM = min(find(S~=1));
+ if isempty(DIM), DIM = 1; end
+end
+if nargin<5 || isempty(WINOPT)
+ WINOPT = {};
+end
+if nargin<6 || isempty(PADOPT)
+ if ypadarray
+  PADOPT = {'replicate'};
+ else
+  PADOPT = {0};
+ end
+elseif ~ypadarray && ~isnumeric(PADOPT{1})
+ PADOPT = {0};
+end
+if nargin<7 || isempty(WNAN)
+ WNAN = 0;
+end
+
+% Selects the 1-dimensional filter or set a row vector: 
+if N==1
+ a = zeros(1,DIM);
+ a(DIM) = F;
+ F = a;
+ clear a
+end
+
+% Checks if the window input is a function or an array:
+if ~isa(HWIN,'function_handle') && ~ischar(HWIN)
+ W = HWIN;
+else
+ W = [];
+end
+
+% If no input data but two outputs then generates the window only:
+if isempty(X)
+ Y = [];
+ if nargout==2 && ~isempty(W)
+  W = ndwind(HWIN,F,WINOPT{:});
+ end
+ return
+end
+
+% Generates the window:
+if isempty(W)
+ W = ndwind(HWIN,F,WINOPT{:});
+end
+
+% Check for NaN's:
+inan = isnan(X);
+ynan = any(inan(:));                       % Bug fixed 30/jun/2008
+if ynan
+ X(inan) = 0;
+else
+ factor = sum(W(:));
+end
+
+% Filtering:
+if yimfilter                                % Use IMFILTER (faster)
+ if ~isfloat(X)
+  X = double(X);
+ end
+ if ~isfloat(W)
+  W = double(W);
+ end
+ if ynan
+  Y = imfilter(X,W       ,PADOPT{:},'conv');
+ else
+  Y = imfilter(X,W/factor,PADOPT{:},'conv');
+ end
+else                                        % Use CONVN
+ % Sets F and S of equal sizes.
+ F = reshape(F,1,N);
+ Nx = numel(S);
+ if N<Nx
+  F(N+1:Nx) = 0;
+ elseif N>Nx
+  S(Nx+1:N) = 1;
+ end
+ F2 = 2*F;
+ % Pads the borders:
+ if ypadarray
+  ind    = padarray(false(S),F2,true     );    % Index of the padding.
+  Y      = padarray(X       ,F2,PADOPT{:});
+ elseif length(PADOPT{1})==1
+  ind2 = cell(N,1);
+  for n = 1:N
+   ind2{n} = F2(n) + (1:S(n)).';
+  end
+  ind          = repmat(true     ,2*F2+S);
+  Y            = repmat(PADOPT{1},2*F2+S);
+  ind(ind2{:}) = false;
+    Y(ind2{:}) = X;
+ else % No padding at all
+  Y    = X;
+  ind  = repmat(false,S); 
+  warning('Ndnanfilter:PaddingOption','Do not perfom any padding.')
+ end
+ % Convolutes both arrays:
+ if ynan
+  Y = convn(Y,W       ,'same');
+ else
+  Y = convn(Y,W/factor,'same');
+ end
+ %  Eliminates the padding:
+ Y(ind) = [];
+ Y      = reshape(Y,S);   
+end
+
+% Estimates the averages when NaNs are present:
+if ynan
+ if yimfilter
+  factor       = imfilter(double(~inan),W,PADOPT{:},'conv');
+ else
+  if ypadarray
+   factor      = padarray(~inan,F2,PADOPT{:});
+  elseif length(PADOPT{1})==1 % (won't work properly with NaNs at borders)
+   factor          = ind;
+   factor(ind2{:}) = ~inan;
+  else
+   factor = ~inan;
+  end
+  factor      = convn(factor,W,'same');
+  factor(ind) = [];
+  factor      = reshape(factor,S);
+ end
+ Y = Y./factor;
+end
+
+% What about NaNs?:
+if     WNAN == 1       % Leave NaNs elements untouched!
+ Y(inan) = NaN;
+elseif WNAN == 2       % Leave no-NaNs elements untouched!!!
+ X(inan) = Y(inan);
+ Y = X;
+end  
+
+
+function W = ndwind(HWIN,F,varargin)
+% NDWIND Generate a N-Dimensional zero-phase window.
+%
+%   Syntax:
+%     W = ndwind(HWIN,F);
+%     W = ndwind(HWIN,F,OPT);
+%
+%   Input:
+%     HWIN - Window function handle. See WINDOW for details. By default
+%            uses: @rectwin (a rectangular window).
+%     F    - A vector specifying the semiwidth of the window for each
+%            dimension. The window's width will be 2*F+1. By default uses:
+%            3 (i.e. a window of width 6). 
+%     OPT  - Cell array specifying optional arguments for the window
+%            function. By default uses: {[]} (window's defaults).
+%
+%   Output:
+%     W    - N-Dimensional window with central symmetry.
+%
+%   Description:
+%     In the axes of each dimension, W has a 1-D window defined as
+%              feval(HWIN,2*F(n)+1), n = 1,...,N.
+%     That is, they are defined by the same window function but have
+%     different widths. So, this program creates another widther window (at
+%     least 201 points), with the same definition, and finds how much the
+%     former windows should be expanded in order to fit the latter one. 
+%
+%     Afterwards, the coordinates of every point are expanded accordingly
+%     and the value of the window in those points are found by linear
+%     interpolation with the bigger window. 
+%
+%     In resume, it is like rotating this big window through every
+%     dimension and then shrinking it through each of its axes to fix the
+%     specified widths.
+%
+%   Notes: 
+%     * Because of the use of the semi-widths F's, all the generated
+%       windows are even. Therefore the change of phase is null. 
+%     * The window function HWIN should output an even function, otherwise,
+%       it won't generate an error but this program will consider only the
+%       last half of it.
+%     * The window should be monotonically decreasing.
+%     * Instead of the handle window, it can be given as a string:
+%       'hamming' instead of @hamming, for example.
+%     * Uses the MATLAB's function FUNC2STR.
+%
+%   Example:
+%     W = ndwind(@hamming,[3 2])
+%     % Results:
+%     W =
+%     
+%              0         0    0.0800         0         0
+%              0    0.1417    0.3100    0.1417         0
+%              0    0.3966    0.7700    0.3966         0
+%         0.0800    0.5400    1.0000    0.5400    0.0800
+%              0    0.3966    0.7700    0.3966         0
+%              0    0.1417    0.3100    0.1417         0
+%              0         0    0.0800         0         0
+%
+%
+%   See also: WINDOW from the Signal Processing Toolbox; and FWIND1,
+%   FWIND2, and FSPECIAL from the Image Processing Toolbox.
+
+%   Copyright 2008 Carlos Adrian Vargas Aguilera
+%   $Revision: 1.1 $  $Date: 2008/06/26 19:30:00 $
+
+%   Written by
+%   M.S. Carlos Adrian Vargas Aguilera
+%   Physical Oceanography PhD candidate
+%   CICESE 
+%   Mexico, 2008
+%   nubeobscura@hotmail.com
+%
+%   Download from:
+%   http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objec
+%   tType=author&objectId=1093874
+
+%   1.0    Release (2008/06/23 10:30:00)
+%   1.1    Fixed Bug adding an extra dimension of unitary width.  
+
+% Check inputs:
+if nargin<1 || isempty(HWIN)
+ HWIN = 'rectwin';
+end
+if nargin<2 || isempty(F)
+ F = 3;
+end
+
+% Rectangular wind?:
+if isa(HWIN,'function_handle')
+ HWIN = func2str(HWIN);
+end
+if strcmpi(HWIN,'rectwin')
+ W = ones([2*F(:).'+1 1]);
+ return
+end
+
+% Generate the BIG window (only the last half):
+FBIG         = max([100; F(:)]);
+BIGw         = feval(HWIN,2*FBIG+1,varargin{:});
+BIGw(1:FBIG) = [];       % Deletes the first half.
+rBIGw        = 0:FBIG;   % Window argument (distance).
+
+% Axial windows widths:
+N  = numel(F);
+F  = reshape(F,1,N); 
+F  = [F 0];             % BUG fixed by adding an extra dimension.
+N  = N+1;
+F2 = 2*F+1;
+
+
+% Pre-allocates the final window and the expanded axis:
+W  = zeros(F2);
+An = cell(N,1);
+Ae = An;
+
+% Generates the index and expanded axes:
+for n = 1:N
+ 
+ % Generate temporally the window in the n-axis:
+ wn = feval(HWIN,F2(n),varargin{:});
+ 
+ % Finds the expansion factors (Note: the window should tends to zero):
+ if F(n)
+  piv = wn(end);
+  ind = (BIGw == piv);
+  if ~any(ind)
+   ind1 = (BIGw >= piv); ind1 = length(ind1(ind1));
+   ind2 = (BIGw <= piv); ind2 = length(ind2(~ind2))+1;
+   if ind2>FBIG+1
+    r = rBIGw(ind1);
+   else
+    r = interp1(BIGw([ind1 ind2]), rBIGw([ind1 ind2]),piv);
+   end
+  else
+   r = rBIGw(ind);
+  end
+  Ef = r/F(n);
+ else
+  Ef = 1;
+ end
+ 
+ % Reversed index and expanded n-axis (for the following grid):
+ An{n} = (F(n):-1:0);
+ Ae{n} = An{n}*Ef;
+ 
+end
+
+% Estimates the expanded distances outside the axes (only at the 1st
+% quarter):
+% Note: In a 2-Dimensional matrix, by the 1st quarter of a matrix I mean
+% the first 1/4 piece of the matrix after you divided it throuh the middle
+% row and column. In N-dimensions it would be the 1st 1/2^N part.
+gride4      = cell(N,1);
+[gride4{:}] = ndgrid(Ae{:});
+R4          = sqrt(sum(reshape([gride4{:}],prod(F+1),N).^2,2));
+
+% Generates the window and linear index in the 1st quarter:
+grid4     = cell(N,1);
+[grid4{:}]= ndgrid(An{:});
+in        = (R4<=rBIGw(end));           % Looks for elements inside window.
+W4        = zeros(F+1);                 % 1st quarter of the window.
+W4(in)    = interp1(rBIGw,BIGw,R4(in)); % Interpolates the window values.
+for n=1:N                               % Linear index on the 1st quarter.
+ grid4{n} = flipdim(grid4{n}+1,n);
+end
+ind4      = sub2ind(F2,grid4{:});
+
+% Index of permutations to fill the N-D window:
+np = 2^N-1;
+ip = zeros(1,np);
+for n = 1:N
+ ini  = 2^(n-1);
+ step = ini*2;
+ ip(ini:step:np) = n;
+end
+
+% Fills the N-D window by flipping W4 and the index: 
+ones4       = repmat(false,F2);    % Avoids using new FALSE function
+ones4(ind4) = true;
+W(ones4)    = W4;
+for kp = ip
+ W4         = flipdim(W4,kp);
+ ones4      = flipdim(ones4,kp);
+ W(ones4)   = W4;
+end
\ No newline at end of file
diff --git a/filter/spikeRemoval.m b/filter/spikeRemoval.m
new file mode 100644
index 0000000..bd1d561
--- /dev/null
+++ b/filter/spikeRemoval.m
@@ -0,0 +1,601 @@
+function [x,spikes] = spikeRemoval(w,varargin)
+%   DESPIKING DISCRETE-TIME SIGNAL USING HISTOGRAM METHOD (a.k.a. DELETE OUTLIERS)
+%
+%   Time series may contain undesired transients and spikes. This function
+%   replaces sudden spikes (outliers) exceeding the threshold value by
+%   either:
+%
+%   (1) interpolating among previous and subsequent data points, or
+%   (2) replacing them with NaN (delete outliers)
+%
+%   User needs to choose one of the above options; default is interpolation.
+%
+%   "w" is the input array and "x" is the spike removed output array. In
+%   addition, output argument "spikes" is a structure array, which returns
+%   the indices of spikes, corresponding values, and replacement values. If
+%   npass parameter is selected as 2 (it means the code will run twice),
+%   then "spikes" will be the structure array of results corresponding to
+%   each run.
+%
+%   The threshold is defined as mean +/- a number of standard deviations of
+%   windowed data centered at spike locations.
+%
+%   USAGE:
+%   [x,spikes] = spikeRemoval(w) 
+%
+%   or
+%   
+%   [x,spikes] = spikeRemoval(w,prop_name,prop_val)
+%
+%   STATIC INPUT:
+%            w = vector of time series data (1xn or nx1)
+%
+%   VALID PROP_NAME / PROP_VAL PAIRS:
+%   -----------------------------------------
+%   'wnsz'   --> (1x1)-[numeric]-[default:25]
+%   'nstd'   --> (1x2)-[numeric]-[default: 3]
+%   'nbins'  --> (1x2)-[numeric]-[default: 4]
+%   'ndata'  --> (1x1)-[numeric]-[default: 5]
+%   'npass'  --> (1x1)-[numeric]-[default: 2]
+%   'method' --> [text]-[default: 'Linear']
+%   'debug'  --> [text]-[default: 'True']
+%
+%  NOTES:
+%         x = spike removed time series
+%
+%    spikes = 3xn structure array, first column contains indexes of spikes,
+%             second column contains values of spikes and third column
+%             contains replacement values
+%
+%      wnsz = window size to compute mean and standard deviation for
+%             thresholding
+%
+%      nstd = number of standard deviations above and below mean
+%
+%     nbins = number of bins used in histogram method
+%
+%     ndata = number of neighbouring data points for spike replacement
+%
+%     npass = number of passes for spike detection and replacement (e.g.,
+%             1); Suggest npass = 2 or 3 for noisy data with back-to-back
+%             spikes
+%
+%    method = interpolation method (linear, spline, cubic, nearest). To
+%             delete outliers chose method as 'delete' 
+%             (e.g., 'method','delete').
+%
+%     debug = 'True' for print debug messages and plot results or '' for no
+%             debug messages or plots
+%
+%   WARNING: This function is for detecting sudden spikes, it may not be
+%            effective detecting long duration spike-like pulses.
+%
+%   EXAMPLE 1 (using defaults values):
+%
+%    w = rand(500,1)*10; w(ceil(rand(1,20)*500))=rand(1,20)*100;
+%    [x,spikes] = spikeRemoval(w);
+%
+%   EXAMPLE 2 (using user defined values for interpolation):
+%
+%    w = rand(500,1)*10; w(ceil(rand(1,20)*500))=rand(1,20)*100;
+%    [x,spikes] = spikeRemoval(w,'wnsz',20,'nstd',2,'npass',1);
+%
+%   EXAMPLE 3 (using user defined values for deleting outliers):
+%
+%    w = rand(500,1)*10; w(ceil(rand(1,20)*500))=rand(1,20)*100;
+%    [x,spikes] = spikeRemoval(w,'nstd',2,'npass',1,'method','delete');
+%
+%   REQUIREMENTS:
+%   spikeRemoval function doesn't require any MatLAB toolbox.
+%
+%   ACKNOWLEDGEMENTS:
+%   I would like thank Dr. Ayal Anis of Texas AM and Jeanne Jones of USGS
+%   for their valuable comments and suggestions, which helped improving
+%   this code significantly. Also, thanks to Jan Glscher of University of
+%   Hamburg for making his nanmean and nanstd codes available. Finally,
+%   thanks to James for making his cent_diff_n function available at MatLAB
+%   FEX.
+%
+%   REFERENCES:
+%
+%   Solomon, O. M., D. R. Larson, and N. G. Paulter (2001). Comparison of
+%   some algorithms to estimate the low and high state level of pulses,
+%   IEEE Instrumentation and Measurement Technology Conference, Vol. 1,
+%   Budapest, Hungary, 21?23 May 2001, 96?101.
+%
+%   cent_diff_n function is available at
+%   https://www.mathworks.com/matlabcentral/fileexchange/36123-n-point-central-differencing
+%   
+%   nanmean and nanstd functions are available at
+%   https://www.mathworks.com/matlabcentral/fileexchange/6837-nan-suite
+%
+%   THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED
+%   WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+%   MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
+%   NO EVENT SHALL THE COPYRIGHT OWNER BE LIABLE FOR ANY DIRECT, INDIRECT,
+%   INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+%   BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+%   OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+%   ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
+%   TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
+%   USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
+%   DAMAGE.
+%
+%   Written by Dr. Erol Kalkan, P.E. (kalkan76@gmail.com)
+%   Revision: 2.1.0
+%   Date: 2019/03/20 08:58:00 $
+
+%% DEFAULT PROPERTIES
+if (nargin >= 1)
+    wnsz = 25;
+    nstd = 3;
+    ndata = 5;
+    nbins = 4;
+    method = 'Linear';
+    debug = 'True';
+    npass = 2;
+else
+    error('spikeRemoval: First argument must be a waveform')
+end
+
+%% USER-DEFINED PROPERTIES
+if (nargin > 1)
+    v = varargin;
+    nv = nargin-1;
+    if ~rem(nv,2) == 0
+        error(['spikeRemoval: Arguments after wave must appear in ',...
+            'property name/val pairs'])
+    end
+    for n = 1:2:nv-1
+        name = lower(v{n});
+        val = v{n+1};
+        switch name
+            case 'nstd'
+                nstd = val;
+            case 'wnsz'
+                wnsz = val;
+            case 'ndata'
+                ndata = val;
+            case 'method'
+                method = val;
+            case 'debug'
+                debug = val;
+            case 'nbins'
+                nbins = val;
+            case 'npass'
+                npass = val;
+            otherwise
+                error('spikeRemoval: Property name not recognized')
+        end
+    end
+end
+
+% Initialize
+nspk  = 0; % number of spikes (outliers)
+x = w;
+
+for i = 1:npass
+    [x,spikes(i).results,nspk] = run(x,nspk,nstd,wnsz,ndata,method,debug,nbins);
+end
+if debug
+    fprintf('Number of spikes found and replaced = %4d\n',nspk);
+    figure(i+1);
+    subplot(311); plot(w,'r'); hold on; title('Original'); lim = ylim;
+    subplot(312); plot(x); title('SpikeRemoval'); ylim(lim);
+    subplot(313); plot(medfilt1_perso(w,nstd),'g');
+    title('MatLAB 1-D Median Filter'); ylim(lim);
+end
+end
+
+function [x,spikes,nspk] = run(x,nspk,nstd,wnsz,ndata,method,debug,nbins)
+% Enforce input as row vector
+if ~isrow(x)
+    x = x';
+end
+
+% Compute first derivative using central differentiation with order 3
+xdot = cent_diff_n(x,1,3);
+
+% Define spike locations via histogram method
+R = statelevel(abs(xdot),nbins);
+locs = find(abs(xdot) > R(1));
+
+if debug
+    figure(1); plot(xdot); hold on; plot(locs,xdot(locs),'ro');
+    title('Spikes to be checked on first derivative of time series');
+end
+
+for i = 1:length(locs)
+    [x,spike,flag] = spikeFix(x,locs(i),nstd,wnsz,ndata,method,debug);
+    nspk = nspk + flag;
+    if nspk ~= 0 && flag == 1
+        spikes(nspk,:) = spike;
+    end
+end
+if ~exist('spikes','var')
+    spikes = -1;
+end
+end
+
+function [x,spike,flag] = spikeFix(x,idx,nstd,wnsz,ndata,method,debug)
+%   Input:
+%         x = vector of time series (1xn)
+%
+%       idx = spike location
+%
+%   Output:
+%         x = vector of spiked removed time series
+%
+%      flag = 1 or 0; 1 means spike is found and fixed, 0 means no spike
+%             detected
+
+% If spike is at the beginning, replace it with the mean value of next
+% ndata points
+if idx <= 2*ndata
+    w = x(idx);
+    if method == 'delete'
+        x(idx) = NaN;
+    else
+        x(idx) = nanmean(x(idx+1:idx+ndata));
+    end
+    flag = 1;
+    spike = [idx,w,x(idx)];
+    if debug
+        fprintf('Peak of %5.2f @%5.0f is replaced by %5.2f\n',w,idx,x(idx));
+    end
+    % If spike is at the end, replace it with the mean value of previous
+    % ndata points
+elseif idx >= length(x)-2*ndata
+    w = x(idx);
+    if method == 'delete'
+        x(idx) = NaN;
+    else
+        x(idx) = nanmean(x(idx-ndata:idx-1));
+    end
+    flag = 1;
+    spike = [idx,w,x(idx)];
+    if debug
+        fprintf('Peak of %5.2f @%5.0f is replaced by %5.2f\n',w,idx,x(idx));
+    end
+else
+    % Set spike search region. Note: arr needs to be updated if n-point
+    % central difference method is used for differentiation, for instance
+    % if 5-point central difference method is used then arr =
+    % (idx-2:idx+1); the current implementation uses a 3-point central
+    % difference
+    arr = (idx-1:idx+1);
+    [~, idx] = max(abs(x(arr)));
+    i = arr(idx);
+    
+    % Overwrite window size if spike occurs within first or last data
+    % window
+    if i <= floor(wnsz/2)
+        wnsz = (i-1)*2;
+    elseif i >= length(x) - floor(wnsz/2)
+        wnsz = length(x) - i;
+    end
+    
+    % Overwrite ndata if it is bigger than half of window size
+    if ndata > floor(wnsz/2)
+        ndata = floor(wnsz/2)-1;
+    end
+    
+    % Compute mean and standard deviation of data window centered at spike
+    % by excluding spike
+    arr1 = [x((i - floor(wnsz/2)):(i-1)), x((i+1):(i + floor(wnsz/2)))];
+    mu_x = nanmean(arr1); std_x = nanstd(arr1);
+    
+    if debug
+        fprintf('wnsz =%5.0f, ndata =%5.0f, i = %5.0f, x(i) = %5.2f, mu_x = %5.2f, std_x = %5.2f\n', ...
+            wnsz,ndata,i,x(i),mu_x,std_x);
+    end
+    
+    % Check threshold exceedance
+    if x(i) >= mu_x + std_x*nstd || x(i) <= mu_x - std_x*nstd
+        w = x(i);
+        
+        if method == 'delete'
+            x(i) = NaN;
+        else
+            % Replace spike region using interpolation of neighbouring data
+            % points
+            x(i-1:i+1) = interp1([i-ndata-1:i-2,i+2:i+ndata+1],...
+                [x(i-ndata-1:i-2),x(i+2:i+ndata+1)],i-1:i+1,method);
+        end
+        flag = 1;
+        spike = [i,w,x(i)];
+        if debug
+            fprintf('Peak of %5.2f @%5.0f is replaced by %5.2f\n',w,i,x(i));
+        end
+    else
+        flag = 0;
+        spike = -1;
+    end
+end
+end
+
+function [levels, histogram, bins] = statelevel(y,n)
+ymax = max(y);
+ymin = min(y)-eps;
+
+% Compute histogram
+idx = ceil(n * (y-ymin)/(ymax-ymin));
+idx = idx(idx>=1 & idx<=n);
+histogram = zeros(n, 1);
+for i=1:numel(idx)
+    histogram(idx(i)) = histogram(idx(i)) + 1;
+end
+
+% Compute center of each bin
+ymin = min(y);
+Ry = ymax-ymin;
+dy = Ry/n;
+bins = ymin + ((1:n)-0.5)*dy;
+
+% Compute state levels
+iLowerRegion = find(histogram > 0, 1, 'first');
+iUpperRegion = find(histogram > 0, 1, 'last');
+
+iLow  = iLowerRegion(1);
+iHigh = iUpperRegion(1);
+
+% Define the lower and upper histogram regions halfway between the lowest
+% and highest nonzero bins.
+lLow  = iLow;
+lHigh = iLow + floor((iHigh - iLow)/2);
+uLow  = iLow + floor((iHigh - iLow)/2);
+uHigh = iHigh;
+
+% Upper and lower histograms
+lHist = histogram(lLow:lHigh, 1);
+uHist = histogram(uLow:uHigh, 1);
+
+levels = zeros(1,2);
+[~, iMax] = max(lHist(2:end));
+[~, iMin] = max(uHist);
+levels(1) = ymin + dy * (lLow + iMax(1) - 1.5);
+levels(2) = ymin + dy * (uLow + iMin(1) - 1.5);
+end
+
+function y = nanmean(x,dim)
+% FORMAT: Y = NANMEAN(X,DIM)
+% 
+%    Average or mean value ignoring NaNs
+%
+%    This function enhances the functionality of NANMEAN as distributed in
+%    the MATLAB Statistics Toolbox and is meant as a replacement (hence the
+%    identical name).  
+%
+%    NANMEAN(X,DIM) calculates the mean along any dimension of the N-D
+%    array X ignoring NaNs.  If DIM is omitted NANMEAN averages along the
+%    first non-singleton dimension of X.
+%
+%    Similar replacements exist for NANSTD, NANMEDIAN, NANMIN, NANMAX, and
+%    NANSUM which are all part of the NaN-suite.
+%
+%    See also MEAN
+
+% -------------------------------------------------------------------------
+%    author:      Jan Glscher
+%    affiliation: Neuroimage Nord, University of Hamburg, Germany
+%    email:       glaescher@uke.uni-hamburg.de
+%    
+%    $Revision: 1.1 $ $Date: 2004/07/15 22:42:13 $
+
+if isempty(x)
+	y = NaN;
+	return
+end
+
+if nargin < 2
+	dim = min(find(size(x)~=1));
+	if isempty(dim)
+		dim = 1;
+	end
+end
+
+% Replace NaNs with zeros.
+nans = isnan(x);
+x(isnan(x)) = 0; 
+
+% denominator
+count = size(x,dim) - sum(nans,dim);
+
+% Protect against a  all NaNs in one dimension
+i = find(count==0);
+count(i) = ones(size(i));
+
+y = sum(x,dim)./count;
+y(i) = i + NaN;
+% $Id: nanmean.m,v 1.1 2004/07/15 22:42:13 glaescher Exp glaescher $
+end
+
+function y = nanstd(x,dim,flag)
+% FORMAT: Y = NANSTD(X,DIM,FLAG)
+% 
+%    Standard deviation ignoring NaNs
+%
+%    This function enhances the functionality of NANSTD as distributed in
+%    the MATLAB Statistics Toolbox and is meant as a replacement (hence the
+%    identical name).  
+%
+%    NANSTD(X,DIM) calculates the standard deviation along any dimension of
+%    the N-D array X ignoring NaNs.  
+%
+%    NANSTD(X,DIM,0) normalizes by (N-1) where N is SIZE(X,DIM).  This make
+%    NANSTD(X,DIM).^2 the best unbiased estimate of the variance if X is
+%    a sample of a normal distribution. If omitted FLAG is set to zero.
+%    
+%    NANSTD(X,DIM,1) normalizes by N and produces the square root of the
+%    second moment of the sample about the mean.
+%
+%    If DIM is omitted NANSTD calculates the standard deviation along first
+%    non-singleton dimension of X.
+%
+%    Similar replacements exist for NANMEAN, NANMEDIAN, NANMIN, NANMAX, and
+%    NANSUM which are all part of the NaN-suite.
+%
+%    See also STD
+
+% -------------------------------------------------------------------------
+%    author:      Jan Glscher
+%    affiliation: Neuroimage Nord, University of Hamburg, Germany
+%    email:       glaescher@uke.uni-hamburg.de
+%    
+%    $Revision: 1.1 $ $Date: 2004/07/15 22:42:15 $
+
+if isempty(x)
+	y = NaN;
+	return
+end
+
+if nargin < 3
+	flag = 0;
+end
+
+if nargin < 2
+	dim = min(find(size(x)~=1));
+	if isempty(dim)
+		dim = 1; 
+	end	  
+end
+
+
+% Find NaNs in x and nanmean(x)
+nans = isnan(x);
+avg = nanmean(x,dim);
+
+% create array indicating number of element 
+% of x in dimension DIM (needed for subtraction of mean)
+tile = ones(1,max(ndims(x),dim));
+tile(dim) = size(x,dim);
+
+% remove mean
+x = x - repmat(avg,tile);
+
+count = size(x,dim) - sum(nans,dim);
+
+% Replace NaNs with zeros.
+x(isnan(x)) = 0; 
+
+
+% Protect against a  all NaNs in one dimension
+i = find(count==0);
+
+if flag == 0
+	y = sqrt(sum(x.*x,dim)./max(count-1,1));
+else
+	y = sqrt(sum(x.*x,dim)./max(count,1));
+end
+y(i) = i + NaN;
+
+% $Id: nanstd.m,v 1.1 2004/07/15 22:42:15 glaescher Exp glaescher $
+end
+
+function df = cent_diff_n(f,h,n)
+% df = cent_diff_n(f,h,n)
+% Computes an n-point central difference of function f with spacing h.
+% Returns a vector df of same size as f.
+% Input f must be a vector with evenly spaced points.
+% Input n must be 3,5,7, or 9.
+% All three inputs are required.
+%
+% Differences for points near the edges are calculated with lower order.
+% For example, if n=5 and length(f)=10, then 3-point central differencing is used
+% to calculate values at points 2 and 9, 2-point forward differencing is used for
+% point 1, 2-point backward differencing is used for point 10, and 5-point central
+% differencing is used for points 3-7.
+%
+% If f contains less than n points, the order will be downgraded to the
+% maximum possible.  Ex: if length(f) = 6, n will be downgraded to 5.
+%
+% Differencing formulae from: http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/central-differences/
+% Accessed 4/10/12.
+
+if nargin < 3
+    error('Not enough inputs.  See help documentation.')
+end
+
+if ~isscalar(h)
+    error('Input h must be a scalar value.')
+end
+
+possible_ns = [3,5,7,9];
+if ~ismember(n,possible_ns)
+    error('Input n must be 3,5,7, or 9')
+end
+
+numPts = length(f);
+if numPts < n
+    newN = max(possible_ns(possible_ns<=numPts));
+    warnstr = [num2str(n) '-point differencing was requested,\n'...
+        'but input function only has ' num2str(numPts) ' points.\n'...
+        'Switching to ' num2str(newN) '-point differencing.'];
+    warning(warnstr,'%s')
+    n = newN;
+end
+
+df_1 = b_diff(f,h);
+df_End = f_diff(f,h);
+
+% Calculate 3-point for all
+df_3pt = c_diff(f,h,3);
+
+if n >=5
+    df_5pt = c_diff(f,h,n);
+    % For the 2nd and next-to-last grid point, use 3-point differencing.
+    df_2 = df_3pt(1);
+    df_Endm1 = df_3pt(end);
+end
+if n >=7
+    df_7pt = c_diff(f,h,7);
+    % For the 3nd and 2nd from last grid point, use 5-point differencing.
+    df_3 = df_5pt(1);
+    df_Endm2 = df_5pt(end);
+end
+if n>= 9
+    df_9pt = c_diff(f,h,9);
+    % For the 4nd and 3rd from last grid point, use 7-point differencing.
+    df_4 = df_7pt(1);
+    df_Endm3 = df_7pt(end);
+end
+
+switch n
+    case 3
+        df = [df_1 df_3pt df_End];
+    case 5
+        df = [df_1 df_2 df_5pt df_Endm1 df_End];
+    case 7
+        df = [df_1 df_2 df_3 df_7pt df_Endm2 df_Endm1 df_End];
+    case 9
+        df = [df_1 df_2 df_3 df_4 df_9pt df_Endm3 df_Endm2 df_Endm1 df_End];
+end
+end
+
+function df = c_diff(f,h,n)
+midStartPoint = ceil(n/2); % First point at which full n points can be used
+midEndPoint = length(f)-midStartPoint+1; % Last point at which full n points can be used
+
+df = [];
+for k = midStartPoint:midEndPoint
+    switch n
+        case 3
+            df_k = (f(k+1) - f(k-1))/(2*h);
+        case 5
+            df_k = (f(k-2) - 8*f(k-1) + 8*f(k+1) - f(k+2))/(12*h);
+        case 7
+            df_k = (-f(k-3) + 9*f(k-2) - 45*f(k-1) + 45*f(k+1) - 9*f(k+2) + f(k+3))/(60*h);
+        case 9
+            df_k = (3*f(k-4) - 32*f(k-3) + 168*f(k-2) - 672*f(k-1) + 672*f(k+1) - 168*f(k+2) + 32*f(k+3) - 3*f(k+4))/(840*h);
+    end
+    df = [df df_k];
+end
+end
+
+function df1 = b_diff(f,h)
+df1 = (f(2)-f(1))/h;
+end
+
+function dfEnd = f_diff(f,h)
+dfEnd = (f(end)-f(end-1))/h;
+end
\ No newline at end of file
diff --git a/moored_adcp_proc/adcp_check_surface.m b/moored_adcp_proc/adcp_check_surface.m
index 3f19322..2dec06f 100644
--- a/moored_adcp_proc/adcp_check_surface.m
+++ b/moored_adcp_proc/adcp_check_surface.m
@@ -10,7 +10,7 @@ set(gca,'YTick',bins);
 ylabel('Bins');
 hold on;
 gregtick;
-title('Meridional velocity before correction');
+title('Zonal velocity before correction');
 
 subplot 323;
 imagesc(time,bin(bins),u1(bins,:),[-1 1]);
@@ -18,7 +18,7 @@ set(gca,'YTick',bins);
 hold on;
 ylabel('Bins');
 gregtick;
-title('Meridional velocity after correction');
+title('Zonal velocity after correction');
 
 subplot 322;
 imagesc(time,bin(bins),v(bins,:),[-1 1]);
@@ -26,7 +26,7 @@ set(gca,'YTick',bins);
 hold on;
 ylabel('Bins');
 gregtick;
-title('Zonal velocity before correction');
+title('Meridional velocity before correction');
 
 subplot 324;
 imagesc(time,bin(bins),v1(bins,:),[-1 1]);
@@ -34,10 +34,10 @@ set(gca,'YTick',bins);
 hold on;
 ylabel('Bins');
 gregtick;
-title('Zonal velocity after correction');
+title('Meridional velocity after correction');
 
 subplot 325;
-for i=1:length(bins);
+for i=1:length(bins)
     plot(time,hammfilter_nodec(z(bins(i),:),73),'color','k');
     hold on;
 end
@@ -47,7 +47,7 @@ ylabel('Bins');
 title('Bin depth after correction');
 
 subplot 326;
-for i=1:length(bins);
+for i=1:length(bins)
     plot(time,hammfilter_nodec(z(bins(i),:),73),'color','k');
     hold on;
 end
diff --git a/moored_adcp_proc/adcp_filt_sub.m b/moored_adcp_proc/adcp_filt_sub.m
index 310e7c7..ae84377 100644
--- a/moored_adcp_proc/adcp_filt_sub.m
+++ b/moored_adcp_proc/adcp_filt_sub.m
@@ -1,77 +1,112 @@
-function [uintfilt,vintfilt,inttim]=adcp_filt_sub(data,ui,vi,intdepvec,nhours)
+function [uintfilt,vintfilt,inttim,utid_baro,vtid_baro]=adcp_filt_sub(data,ui,vi,intdepvec,nhours)
 
-sf=1/((data.time(2)-data.time(1))*24);
+%sf = 1/((data.time(2)-data.time(1))*24);
 
 for iidep = 1:length(intdepvec)
     
-  adcpui = ui(iidep,:);
-  adcpvi = vi(iidep,:);
-  valid = find(~isnan(adcpui));
+  adcpui  = ui(iidep,:);
+  adcpvi  = vi(iidep,:);
+  valid   = find(~isnan(adcpui));
   invalid = find(isnan(adcpui));
  
-  if length(valid)>1 && length(valid(1):valid(end))>240; % length(invalid)>900 & 
+  if length(valid)>1 && length(valid(1):valid(end))>240 % length(invalid)>900 & 
 
       % Horizontal interpolation to fill gaps
-      timval = data.time(valid);
-      adcpui = interp1(timval,adcpui(valid),data.time);
-      adcpvi = interp1(timval,adcpvi(valid),data.time);
+      timval                           = data.time(valid);
+      adcpui                           = interp1(timval,adcpui(valid),data.time);
+      adcpvi                           = interp1(timval,adcpvi(valid),data.time);
       sum(isnan(adcpui));
       sum(isnan(adcpvi));
-      % filtering
-      adcpui(valid(1):valid(end)) = mfilter(adcpui(valid(1):valid(end)),sf,1/nhours,0,2*nhours,1);
-      adcpvi(valid(1):valid(end)) = mfilter(adcpvi(valid(1):valid(end)),sf,1/nhours,0,2*nhours,1);
-      sum(isnan(adcpui));
-      sum(isnan(adcpvi));
-%       adcpui(invalid) = NaN;
-%       adcpvi(invalid) = NaN;
-      sum(isnan(adcpui));
-      uifilt(iidep,1:length(adcpui)) = adcpui;
-      vifilt(iidep,1:length(adcpvi)) = adcpvi;
-      % subsampling
-      inttim = [ceil(data.time(1)):0.25:floor(data.time(end))];
+      
+      
+      %% Hamming filter
+      adcpui(valid(1):valid(end))      = mfilter(adcpui(valid(1):valid(end)),sf,1/nhours,0,2*nhours,1);
+      adcpvi(valid(1):valid(end))      = mfilter(adcpvi(valid(1):valid(end)),sf,1/nhours,0,2*nhours,1);
+      
+      %% Fourier filter
+%       [adcpui(valid(1):valid(end))]  = fft_filter(adcpui(valid(1):valid(end)),2,[24 Inf]);
+%       [adcpvi(valid(1):valid(end))]  = fft_filter(adcpvi(valid(1):valid(end)),2,[24 Inf]);
+          
+      
+      %% subsampling
+      uifilt(iidep,1:length(adcpui))   = adcpui;
+      vifilt(iidep,1:length(adcpvi))   = adcpvi;
+      inttim                           = [ceil(data.time(1)):0.25:floor(data.time(end))];
       uintfilt(iidep,1:length(inttim)) = interp1(data.time,transpose(uifilt(iidep,:)),inttim);
       vintfilt(iidep,1:length(inttim)) = interp1(data.time,transpose(vifilt(iidep,:)),inttim);
-      inttim = inttim;  
    
   elseif length(valid)<=1 || length(valid(1):valid(end))<=240
-      % subsampling
-      uifilt(iidep,1:length(adcpui)) = nan;
-      vifilt(iidep,1:length(adcpvi)) = nan;
-      inttim = [ceil(data.time(1)):0.25:floor(data.time(end))];
-      uintfilt(iidep,1:length(inttim)) =nan;
-      vintfilt(iidep,1:length(inttim)) =nan;
+      
+      %% subsampling
+      uifilt(iidep,1:length(adcpui))   = NaN;
+      vifilt(iidep,1:length(adcpvi))   = NaN;
+      inttim                           = [ceil(data.time(1)):0.25:floor(data.time(end))];
+      uintfilt(iidep,1:length(inttim)) = NaN;
+      vintfilt(iidep,1:length(inttim)) = NaN;
+      
       
   else
-      % Horizontal interpolation to fill gaps
-      timval = data.time(valid);
-      adcpui = interp1(timval,adcpui(valid),data.time);
-      adcpvi = interp1(timval,adcpvi(valid),data.time); 
-      % filtering
-      adcpui = mfilter(adcpui,sf,1/nhours,0,2*nhours,1);
-      adcpvi = mfilter(adcpvi,sf,1/nhours,0,2*nhours,1);
-%       adcpui(invalid) = NaN;
-%       adcpvi(invalid) = NaN;
-      uifilt(iidep,1:length(adcpui)) = adcpui;
-      vifilt(iidep,1:length(adcpvi)) = adcpvi;
-      %subsampling
-      inttim = [ceil(data.time(1)):0.25:floor(data.time(end))];
+      
+      %% Horizontal interpolation to fill gaps
+      timval                           = data.time(valid);
+      adcpui                           = interp1(timval,adcpui(valid),data.time);
+      adcpvi                           = interp1(timval,adcpvi(valid),data.time); 
+      
+      %% Hamming filter
+      adcpui                           = mfilter(adcpui,sf,1/nhours,0,2*nhours,1);
+      adcpvi                           = mfilter(adcpvi,sf,1/nhours,0,2*nhours,1);
+
+      %% Subsampling
+      uifilt(iidep,1:length(adcpui))   = adcpui;
+      vifilt(iidep,1:length(adcpvi))   = adcpvi;
+      inttim                           = [ceil(data.time(1)):0.25:floor(data.time(end))];
       uintfilt(iidep,1:length(inttim)) = interp1(timval,transpose(uifilt(iidep,valid)),inttim);
       vintfilt(iidep,1:length(inttim)) = interp1(timval,transpose(vifilt(iidep,valid)),inttim);
-      inttim = inttim;
+      inttim                           = inttim;
       
   end
       
 end
 
+% [uintfilt] = fft_filter(uintfilt,2,[40 Inf]);
+% [vintfilt] = fft_filter(vintfilt,2,[40 Inf]);
+% uintfilt   = real(uintfilt);
+% vintfilt   = real(vintfilt);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Correct tide using model
+% %% Calc tide
+% addpath(genpath('/home/proussel/Documents/OUTILS/TOOLS/TMD'))
+% addpath(genpath('/home/proussel/Documents/OUTILS/TOOLS/TIDES/tpxo9_atlas'))
+% addpath(genpath('/home/proussel/Documents/OUTILS/TOOLS/TIDES/tpxo8_atlas'))
+% 
+% % u component
+% %utid      = tmd_tide_pred('/home/proussel/Documents/OUTILS/TOOLS/TIDES/tpxo8_atlas_compact/Model_tpxo8.v1',datenum(julian(inttim')),data.lat,data.lon,'u')/100;
+% %utid      = tmd_tide_pred('C:\Users\proussel\Documents\outils\TOOLS\TIDES\tpxo9.1\Model_tpxo9.v1',datenum(julian(inttim')),data.lat,data.lon,'u')/100;
+% utid      = tmd_tide_pred('/home/proussel/Documents/OUTILS/TOOLS/TMD/DATA/Model_tpxo7.2',datenum(julian(inttim')),data.lat,data.lon,'u')/100;
+% utid_baro = utid'.*ones(length(uintfilt(:,1)),1);
+% % v component
+% %vtid      = tmd_tide_pred('/home/proussel/Documents/OUTILS/TOOLS/TIDES/tpxo8_atlas_compact/Model_tpxo8.v1',datenum(julian(inttim')),data.lat,data.lon,'v')/100;
+% %vtid      = tmd_tide_pred('C:\Users\proussel\Documents\outils\TOOLS\TIDES\tpxo9.1\Model_tpxo9.v1',datenum(julian(inttim')),data.lat,data.lon,'v')/100;
+% vtid      = tmd_tide_pred('/home/proussel/Documents/OUTILS/TOOLS/TMD/DATA/Model_tpxo7.2',datenum(julian(inttim')),data.lat,data.lon,'v')/100;
+% vtid_baro = vtid'.*ones(length(vintfilt(:,1)),1);
+
+% %% Correct tide
+% % uintfilt = uintfilt-utid_baro;
+% % vintfilt = vintfilt-vtid_baro;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+utid_baro = 0;
+vtid_baro = 0;
+
 hf=figure('position', [0, 0, 1400, 1000]);
 subplot 121;
-imagesc(data.time,intdepvec,ui,[-1 1]);
+imagesc(data.time,intdepvec,ui(intdepvec,:),[-1 1]);
 %set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10]);
 gregtick;
 ylabel('Bins');
 title('U field - data raw');
 subplot 122;
-imagesc(data.time,intdepvec,uintfilt,[-1 1]);
+imagesc(inttim,intdepvec,uintfilt,[-1 1]);
 %set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10]);
 gregtick;
 ylabel('Bins');
@@ -79,14 +114,16 @@ title('U field - data interpolated, filtered and subsampled');
 
 hf=figure('position', [0, 0, 1400, 1000]);
 subplot 121;
-imagesc(data.time,intdepvec,vi,[-1 1]);
+imagesc(data.time,intdepvec,vi(intdepvec,:),[-1 1]);
 %set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10])
 gregtick;
 ylabel('Bins');
 title('V field - data raw');
 subplot 122;
-imagesc(data.time,intdepvec,vintfilt,[-1 1]);
+imagesc(inttim,intdepvec,vintfilt,[-1 1]);
 %set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10]);
 gregtick;
 ylabel('Bins');
 title('V field - data interpolated, filtered and subsampled');
+
+end
diff --git a/moored_adcp_proc/adcp_shadowzone.m b/moored_adcp_proc/adcp_shadowzone.m
index dc5ba49..3fa8432 100644
--- a/moored_adcp_proc/adcp_shadowzone.m
+++ b/moored_adcp_proc/adcp_shadowzone.m
@@ -1,4 +1,4 @@
-function z_shadowzone=adcp_shadowzone(adcp_dpt,beamangle);
+function z_shadowzone=adcp_shadowzone(adcp_dpt,beamangle)
 
 % z_shadowzone=adcp_shadowzone(adcp_dpt,beamangle);
 %
@@ -8,4 +8,6 @@ function z_shadowzone=adcp_shadowzone(adcp_dpt,beamangle);
 %
 % R. Kopte, 2014/11/20
 
-z_shadowzone=adcp_dpt*(1-cosd(beamangle));
+z_shadowzone = adcp_dpt*(1-cosd(beamangle));
+
+end
diff --git a/moored_adcp_proc/adcp_surface_fit.m b/moored_adcp_proc/adcp_surface_fit.m
index 560fe2c..f0dc31a 100644
--- a/moored_adcp_proc/adcp_surface_fit.m
+++ b/moored_adcp_proc/adcp_surface_fit.m
@@ -1,49 +1,60 @@
-function [zbins,zadcp1,offset,x_null]=adcp_surface_fit(zadcp,ea,surface_bins,blen,blnk,nbin);
+function [zbins,zadcp1,offset,x_null]=adcp_surface_fit(zadcp,ea,surface_bins,blen,blnk,nbin)
     
     % Bin depth matrix
-    dpt1 = repmat(zadcp,nbin,1);
+    dpt1   = repmat(zadcp,nbin,1);
     binmat = repmat((1:nbin)',1,length(dpt1));   
-    z = dpt1-(binmat-0.5)*blen-blnk; 
+    z      = dpt1-(binmat-0.5)*blen-blnk; 
     
     % Loop over time, determine bin of maximum ea in surface bin range and 
     % do quadratic fit over 2 neighbouring and center bins   
-    easurf=ea(surface_bins,:);
+    easurf = ea(surface_bins,:);
     for ii=1:length(ea)
-        [eamax,maxind]=max(easurf(:,ii));
+        [eamax,maxind] = max(easurf(:,ii));
         
-        if length(eamax)>1; maxind=maxind(1); end
+        if length(eamax)>1 
+            maxind = maxind(1); 
+        end
         
-        if maxind>1 & eamax>80
-            if surface_bins(maxind)==nbin; 
+        if maxind>1 && eamax>80
+            if surface_bins(maxind)==nbin 
                 xtofit(1:2) = easurf(maxind-1:maxind,ii);
-                xtofit(3) = 0;
+                xtofit(3)   = 0;
             else           
-                xtofit = easurf(maxind-1:maxind+1,ii);
+                xtofit      = easurf(maxind-1:maxind+1,ii);
             end
 
             for jj=1:3
-                A(jj,:)=[(surface_bins(maxind)+jj-2)^2, surface_bins(maxind)+jj-2, 1];
+                A(jj,:) = [(surface_bins(maxind)+jj-2)^2, surface_bins(maxind)+jj-2, 1];
             end 
-            coef(:,ii)= A\xtofit;
+            coef(:,ii) = A\xtofit;
         else
-            coef(:,ii)=NaN;
+            coef(:,ii) = NaN;
         end
     end
     
     % Find maximum of quadratic fit ax^2+bx+c: 2ax+b=0
     x_null = -coef(2,:)./2./coef(1,:);
-
-    %offset=-round(nmedian((x_null-0.5)*blen+blnk)-nmedian(zadcp));
-    offset=-round(((x_null)*blen+blnk)-(zadcp)); %% P. Rousselot - offset over time
+    
+    
+    %% Calculate offset
+    offset = round(((x_null-0.5)*blen+blnk)-(zadcp));     
     disp('-------------------------------');
-    disp(['Depth offset is ' num2str(nanmean(offset)) ' m']);
+    disp(['Depth offset is ' num2str(round(nanmean(offset))) ' m']);
     disp('-------------------------------');
     
+    % offset over time cleaned (median filter)
+    [offset_clean,~] = clean_median(offset,20,2.8,[0.5 5],2,NaN);
+    lin_offset       = linspace(1,length(offset),length(offset));
+    offset           = interp1(lin_offset(~isnan(offset_clean)),...
+        offset_clean(~isnan(offset_clean)),lin_offset,'linear','extrap');
+    % offset median
+    %offset = nanmedian(offset);
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
-    % Plot histogram of differences
-    dz=((x_null)*blen+blnk)-zadcp;
-    count=[-100:1:100];
-    ncount=hist(-dz,count);
+    %% Plot histogram of differences
+    dz     = ((x_null-0.5)*blen+blnk)-zadcp;
+    count  = [-100:1:100];
+    ncount = hist(-dz,count);
     figure(1);
     bar(count,ncount);
     grid on
@@ -53,38 +64,33 @@ function [zbins,zadcp1,offset,x_null]=adcp_surface_fit(zadcp,ea,surface_bins,ble
     plot(-zadcp,'b');
     grid on
     hold on;
-    plot(-(x_null)*blen+blnk,'r');
+    plot(-((x_null-0.5)*blen+blnk),'r');
     legend('Original','Reconstructed from surface reflection');
     
-%     if abs(offset)>15
-%        reply = input('Do you want to overwrite offset? 1/0 [0]:');
-%        if isempty(reply)
-%           reply = 0;
-%        end
-%        
-%        if reply==1
-%            offset=input('Enter new offset:');
-%        end
-%     end
+    if abs(offset)>15
+       reply = input('Do you want to overwrite offset? 1/0 [0]:');
+       if isempty(reply)
+          reply = 0;
+       end
+       
+       if reply==1
+           offset=input('Enter new offset:');
+       end
+    end
     
-    disp(['Offset of ' num2str(offset) ' m is applied']);
-    offset = 19;
-    % Apply offset to get correct bin depth and instrument depth:
-    zbins=z-offset;
-    zadcp1=zadcp-offset;
+    %% Apply offset to get correct bin depth and instrument depth:
+    zbins  = z+offset;
+    zadcp1 = zadcp+offset;
     
     figure(2);
     plot(-zadcp1,'y');
     text(300, max(zadcp),['Offset applied: ' num2str(offset) ' m']);
-    %,'fonts',12,'fontw','bold','backgroundc','w');
-    legend('Original','Reconstructed from surface reflection','Offset applied');
-    
+    legend('Original','Reconstructed from surface reflection','Offset applied');   
     %print -dpng surface_fit;
     
     figure(3);
     pcolor([1:length(x_null)],-zbins,ea); shading flat;
     title('Amplitude'); colorbar; ylabel('Depth [m]');xlabel('Time index');
-    
     %print -dpng surface_ea;
 
 end
\ No newline at end of file
diff --git a/moored_adcp_proc/filtfilt.m b/moored_adcp_proc/filtfilt.m
index d4e2acd..6b09ada 100644
--- a/moored_adcp_proc/filtfilt.m
+++ b/moored_adcp_proc/filtfilt.m
@@ -56,7 +56,7 @@
 %   Copyright 1988-2011 The MathWorks, Inc.
 %   $Revision: 1.7.4.10 $  $Date: 2012/12/25 21:34:47 $
 
-error(nargchk(3,3,nargin,'struct'))
+narginchk(3,3)
 
 % Only double precision is supported
 if ~isa(b,'double') || ~isa(a,'double') || ~isa(x,'double')
diff --git a/template_get_adcp_data.m b/template_get_adcp_data.m
index d68e946..205e94a 100644
--- a/template_get_adcp_data.m
+++ b/template_get_adcp_data.m
@@ -1,147 +1,158 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % template_get_adcp_data.m
 % -------------------------------
-% Author : Jérémie HABASQUE - IRD
+% Author : Jeremie HABASQUE / Pierre Rousselot - IRD
 % -------------------------------
 % INPUTS:
 % - binary raw file with .000 extension
 % OUTPUTS:
 % - U and V fields interpolated on a regulard grid, filtered and subsampled
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+close all; clear
+addpath(genpath('../ADCP_mooring_data_processing'));
+addpath('/home/proussel/Documents/OUTILS/TOOLS/nansuite'); % NaNSuitePath
 
 % First part --------------------------------------------------------------------------------------------------------------------
-close all
-clear all
- 
 %% META information:
-
-% Path
-addpath('.\moored_adcp_proc'); % ou par exemple ('C:\Users\IRD_US_IMAGO\TRAITEMENTS\ADCP_MOUILLAGE\01_DATA_PROCESSING\moored_adcp_proc');
-addpath('.\backscatter'); % (Optionnel) / ou par exemple ('C:\Users\IRD_US_IMAGO\TRAITEMENTS\ADCP_MOUILLAGE\01_DATA_PROCESSING\backscatter');
-
 % Location rawfile
-fpath = '';
-rawfile='.\FR26_000.000'; % binary file with .000 extension
+rawfile          = './FR26_000.000';        % binary file with .000 extension
+fpath_output     = './FR28/';               % Output directory
  
-% Directory for outputs
-fpath_output = '.\FR28\';
-
 % Cruise/mooring info
-cruise.name = 'PIRATA-FR28';
-mooring.name='0N0E'; % 0N10W par exemple
-mooring.lat=00+00/60; %latitude en degrés décimaux
-mooring.lon=00+00/60; %longitude en degrés décimaux
-clock_drift = 253/3600; % convert into hrs
+cruise.name      = 'PIRATA-FR28';           % cruise name
+mooring.name     = '0N0W';                  % '0N10W'
+mooring.lat      = '00°01.060';             % latitude [°']
+mooring.lon      = '-000°00.330';           % longitude [°']
+clock_drift      = 0/3600;                  % [seconds]
 
 % ADCP info
-adcp.sn=22545;
-adcp.type='150 khz Quartermaster'; % Type : ‘Quartermaster’, ‘longranger’
-adcp.direction='up';        % upward-looking 'up', downward-looking 'dn'
-adcp.instr_depth=300;       % nominal instrument depth
-instr=1;                    % this is just for name convention and sorting of all mooring instruments
-
-% If ADCP was not set up to correct for magnetic deviation internally
-% ("EA0" code in configuration file), use http://www.ngdc.noaa.gov/geomag-web/#declination
-% Magnetic deviation: Mean of deviations at time of deployment and time of recovery 
-
-% Magnetic deviation values
-magnetic_deviation_ini = -5.27;
-magnetic_deviation_end = -4.99;
-rot=(magnetic_deviation_ini+magnetic_deviation_end)/2;  
-
-% Read rawfile
+adcp.sn          = 15258;                   % ADCP serial number
+adcp.type        = '150 khz Quartermaster'; % Type : Quartermaster, longranger
+adcp.direction   = 'up';                    % upward-looking 'up', downward-looking 'dn'
+adcp.instr_depth = 300;                     % nominal instrument depth
+instr            = 1;                       % this is just for name convention and sorting of all mooring instruments
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Convert variables
+latDegInd               = strfind(mooring.lat,'°');  
+lonDegInd               = strfind(mooring.lon,'°');  
+mooring.lat             = str2double(mooring.lat(1:latDegInd-1))+str2double(mooring.lat(latDegInd+1:end-1))/60;
+mooring.lon             = str2double(mooring.lon(1:lonDegInd-1))+str2double(mooring.lon(lonDegInd+1:end-1))/60;
+clock_drift             = clock_drift/3600;  % convert into hrs
+
+%% Read rawfile
 fprintf('Read %s\n', rawfile);
-raw=read_os3(rawfile,'all');
-
-% Correct clock drift
-time0          = julian(raw.juliandate);
-clockd         = linspace(0, clock_drift, length(time0));
-raw.juliandate = raw.juliandate - clockd / 24;                 % P. Rousselot - change clock drift /24
-
-% Magnetic devitation over time %P . Rousselot
-mag_dev = linspace(magnetic_deviation_ini, magnetic_deviation_end, length(time0)); 
-
-
-figure;plot(raw.pressure);
+raw                     = read_os3(rawfile,'all');
+
+%% Correct clock drift
+time0                   = julian(raw.juliandate);
+clockd                  = linspace(0, clock_drift, length(time0));
+raw.juliandate          = raw.juliandate - clockd / 24;              
+disp('Correct clock drift')
+
+%% Calculate Magnetic deviation values
+[a,~]                   = gregorian(raw.juliandate(1));
+magnetic_deviation_ini  = magdev(mooring.lat,mooring.lon,0,a+(raw.juliandate(1)-julian(a,1,1,0,0,0))/365.25);
+[a,~]                   = gregorian(raw.juliandate(end));
+magnetic_deviation_end  = magdev(mooring.lat,mooring.lon,0,a+(raw.juliandate(end)-julian(a,1,1,0,0,0))/365.25);
+rot                     = (magnetic_deviation_ini+magnetic_deviation_end)/2;  
+mag_dev                 = linspace(magnetic_deviation_ini, magnetic_deviation_end, length(time0)); 
+disp('Correct magnetic deviation')
+
+%% Plot pressure and temperature sensor
+figure;
+subplot(2,1,1)
+plot(raw.pressure);
 detrend_sdata = detrend(raw.pressure);
-trend = raw.pressure - detrend_sdata;
+trend         = raw.pressure - detrend_sdata;
 hold on
 plot(trend, 'r--')
 hold off
-title('Pressure sensor');ylabel('Depth [m]');xlabel('Time index');grid on; 
+title('Pressure sensor');
+ylabel('Depth [m]');
+xlabel('Time index');
+grid on; 
 saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Pressure_sensor'],'fig')
 
-figure;plot(raw.temperature);
-title('Temperature sensor');ylabel('Temperature [°C]');xlabel('Time index');grid on; 
+subplot(2,1,2)
+plot(raw.temperature);
+title('Temperature sensor');
+ylabel('Temperature [°C]');
+xlabel('Time index');
+grid on; 
 % Second part --------------------------------------------------------------------------------------------------------------------
 
-% Determine first and last indiced when instrument was at depth (you can do this by plotting 'raw.pressure' for example           
-first = 11; 
-last = 17161; 
-
-% amplitude of the bins / Correction ADCP's depth
-ea = squeeze(mean(raw.amp(:,:,first:last),2));   
-figure; colormap jet; pcolor(ea); shading flat; title('Amplitude of the bins'); colorbar;
-ylabel('Bins');xlabel('Time index'); 
+%% Determine first and last indiced when instrument was at depth (you can do this by plotting 'raw.pressure' for example            
+first               = input('Determine first indice when instrument was at depth (with pres/temp plot): ');
+last                = input('Determine last indice when instrument was at depth (with pres/temp plot): ');
+
+%% amplitude of the bins / Correction ADCP's depth
+ea                  = squeeze(mean(raw.amp(:,:,first:last),2));   
+figure; 
+colormap jet; 
+pcolor(ea); 
+shading flat; 
+title('Amplitude of the bins'); colorbar;
+ylabel('Bins');
+xlabel('Time index'); 
 saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Amplitude_bins'],'fig')
 
 % Third part --------------------------------------------------------------------------------------------------------------------
 
-% If upward looking: range of surface bins used for instrument depth correction below!
-sbins= 29:34;%30:35; % here a range of bins is given which cover the surface reflection
-
-% Exclude data with percent good below prct_good
-prct_good = 20;
-
-%% Read data
-freq = raw.config.sysconfig.frequency;
-
-u2 = squeeze(raw.vel(:,1,first:last));
-v2 = squeeze(raw.vel(:,2,first:last));
-w = squeeze(raw.vel(:,3,first:last));
-vel_err = squeeze(raw.vel(:,4,first:last));  % the difference in vertical velocity between the two pairs of transducers 
-time = raw.juliandate(first:last);
-ang = [raw.pitch(first:last) raw.roll(first:last) raw.heading(first:last)]; 
-soundspeed = raw.soundspeed(first:last);
-T = raw.temperature(first:last);
-press = raw.pressure(first:last);
-mag_dev = mag_dev(first:last);
-
-nbin = raw.config.ncells;  % number of bins
-bin  = 1:nbin;
-blen = raw.config.cell;    % bin length
-blnk = raw.config.blank;   % blank distance after transmit
-
-dt=(time(2)-time(1))*24;   % Sampling interval in hours
-
+%% If upward looking: determine range of surface bins used for instrument depth correction below!
+sbins               = input('Determine range of surface bins used for instrument depth correction (with aplitude plot, ie. 30:35): ');
+
+%% Exclude data with percent good below prct_good
+prct_good           = input('Determine prct_good threshold (generally 20): ');
+
+%% Extract data
+freq                = raw.config.sysconfig.frequency;
+u2                  = squeeze(raw.vel(:,1,first:last));
+v2                  = squeeze(raw.vel(:,2,first:last));
+w                   = squeeze(raw.vel(:,3,first:last));
+vel_err             = squeeze(raw.vel(:,4,first:last));  % the difference in vertical velocity between the two pairs of transducers 
+time                = raw.juliandate(first:last);
+ang                 = [raw.pitch(first:last) raw.roll(first:last) raw.heading(first:last)]; 
+soundspeed          = raw.soundspeed(first:last);
+temp                = raw.temperature(first:last);
+press               = raw.pressure(first:last);
+if press < 0
+    press = -press;
+end
+mag_dev             = mag_dev(first:last);
+nbin                = raw.config.ncells;  % number of bins
+bin                 = 1:nbin;
+blen                = raw.config.cell;    % bin length
+blnk                = raw.config.blank;   % blank distance after transmit
+dt                  = (time(2)-time(1))*24;   % Sampling interval in hours
+
+%% Correction of magnetic deviation
 for ii = 1 : length(mag_dev)
-    [u(:,ii),v(:,ii)]=uvrot(u2(:,ii), v2(:,ii), -mag_dev(ii));   % Correction of magnetic deviation (over time-P. Rousselot)
+    [u(:,ii),v(:,ii)] = uvrot(u2(:,ii), v2(:,ii), -mag_dev(ii));   
 end
 
-pg = squeeze(raw.pg(:,4,first:last));  % percent good
-
-igap=find(pg<prct_good);           % Exclude data with percent good below prct_good
-u(igap)=nan;
-v(igap)=nan;
-w(igap)=nan;
-vel_err(igap)=nan;
+%% Correct percent good
+pg                  = squeeze(raw.pg(:,4,first:last));  
+igap                = find(pg<prct_good);           % Exclude data with percent good below prct_good
+u(igap)             = NaN;
+v(igap)             = NaN;
+w(igap)             = NaN;
+vel_err(igap)       = NaN;
 
 %% Calculate depth of each bin:
-dpt = sw_dpth(press,mooring.lat)';  % convert pressure to depth, press needs to have dimension (n x 1)
-dpt1 = repmat(dpt,nbin,1);
-binmat = repmat((1:nbin)',1,length(dpt1));
+dpt                 = sw_dpth(press,mooring.lat)';  % convert pressure to depth, press needs to have dimension (n x 1)
+dpt1                = repmat(dpt,nbin,1);
+binmat              = repmat((1:nbin)',1,length(dpt1));
 
-% If ADCP is upward-looking a depth correction can be inferred from the
-% surface reflection, which is done in adcp_surface_fit
+% If ADCP is upward-looking a depth correction can be inferred from the surface reflection, which is done in adcp_surface_fit
 if strcmp(adcp.direction,'up')  
-    [z,dpt1,offset,xnull]=adcp_surface_fit(-dpt,ea,sbins,blen,blnk,nbin);
+    [z,dpt1,offset,xnull] = adcp_surface_fit(dpt,ea,sbins,blen,blnk,nbin);
 elseif strcmp(adcp.direction,'dn')
-    z = dpt1+(binmat-0.5)*blen+blnk;
+    z                     = dpt1+(binmat-0.5)*blen+blnk;
 else
     error('Bin depth calculation: unknown direction!');
 end
 
-saveas(figure(1),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Hist_diff_orig-depth_recon-depth'],'fig')
+saveas(figure(1),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Histdiff_depth'],'fig')
 saveas(figure(2),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Offset_depth'],'fig')
 saveas(figure(3),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Amplitude_bins_2'],'fig')
 
@@ -149,30 +160,29 @@ saveas(figure(3),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2s
 u1=u; v1=v; w1=w; vel_err1=vel_err; ea1=ea;
 
 if strcmp(adcp.direction,'up')
-    for i=1:length(time)
-        sz_dpt(i)=adcp_shadowzone(dpt(i),raw.config.sysconfig.angle); % depending on the instrument depth and the beam angle the shadow zone, i.e. the depth below the surface which is contaminated by the surface reflection is determined
-
-        iz(i)=find(z(:,i)<-sz_dpt(i),1,'first')-1;
-        sbin(i)=bin(iz(i));
+    for i = 1:length(time)
+        sz_dpt(i) = adcp_shadowzone(dpt1(i),raw.config.sysconfig.angle); % depending on the instrument depth and the beam angle the shadow zone, i.e. the depth below the surface which is contaminated by the surface reflection is determined
+        iz(i)     = find(z(:,i)>sz_dpt(i),1,'last');
+        sbin(i)   = bin(iz(i));
         
         %sbin(i)=30; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         % here a manual criterion should be hard-coded if
         % adcp_check_surface (below) shows bad velocities close to the
         % surface
 
-        fz(i)=z(iz(i),i);    
+        fz(i)      = z(iz(i),i);    
     end
 
-    for i=1:length(time) 
-        u1(sbin(i)+1:end,i)=nan;
-        v1(sbin(i)+1:end,i)=nan;
-        w1(sbin(i)+1:end,i)=nan;
-        vel_err1(sbin(i)+1:end,i)=nan;
-        ea1(sbin(i)+1:end,i)=nan;
+    for i = 1:length(time) 
+        u1(sbin(i)+1:end,i)       = NaN;
+        v1(sbin(i)+1:end,i)       = NaN;
+        w1(sbin(i)+1:end,i)       = NaN;
+        vel_err1(sbin(i)+1:end,i) = NaN;
+        ea1(sbin(i)+1:end,i)      = NaN;
     end
 
     if 1
-        bins=nmedian(sbin)-4:nmedian(sbin)+2;
+        bins = nmedian(sbin)-4:nmedian(sbin)+2;
         adcp_check_surface(bins,u,u1,v,v1,time,bin,z); 
         % here the closest bins below the surface are plotted that are supposed to have good velocities, if there are still bad velocities a manual criterion needs to be found
     end
@@ -182,59 +192,60 @@ saveas(figure(4),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2s
 %% SAVE DATA
 % More meta information
 %adcp.comment='';
-adcp.config=raw.config;
-adcp.z_offset=offset;
-adcp.ang=ang;
-adcp.mag_dev=rot;
+adcp.config         = raw.config;
+adcp.z_offset       = offset;
+adcp.ang            = ang;
+adcp.mag_dev        = rot;
 
 % Data structure
-data.u=u1;
-data.v=v1;
-data.w=w1;
-data.e=vel_err1;
-data.ea=ea1;
-data.pg=pg;
-data.time=time;
-data.z_bins=z;
-data.depth=dpt;
-data.temp=T;
-data.sspd = soundspeed;
+data.u              = u1;
+data.v              = v1;
+data.w              = w1;
+data.e              = vel_err1;
+data.ea             = ea1;
+data.pg             = pg;
+data.time           = time;
+data.z_bins         = z;
+data.depth          = dpt;
+data.temp           = temp;
+data.sspd           = soundspeed;
+data.lat            = mooring.lat;
+data.lon            = mooring.lon;
 
 % Save
-save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '.mat'],'adcp','mooring','data','raw');
+%save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '.mat'],'adcp','mooring','data','raw','-v7.3');
 
 %% Interpolate data on a regular vertical grid
-Z = fliplr(blen/2:blen:max(z(:))+blen);
-Zmax = max(Z);
-
-u_interp = NaN(length(time),length(Z));
-v_interp = NaN(length(time),length(Z));
+Z                   = fliplr(blen/2:blen:max(z(:))+blen);
+Zmax                = max(Z);
+u_interp            = NaN(length(time),length(Z));
+v_interp            = NaN(length(time),length(Z));
 for i=1:length(time)
     % indice correspondant sur la grille finale Z
     ind = round((Zmax-z(1,i))/blen)+1;
     % filling the grid
-    npts = min([length(Z)-ind+1 length(bin)]);
+    npts = min([length(Z)+ind+1 length(bin)]);
     u_interp(i,ind:ind+npts-1) = u1(1:npts,i);
     v_interp(i,ind:ind+npts-1) = v1(1:npts,i);
 end
   
 %% Horizontal interpolation, filtering and subsampling
-[uintfilt,vintfilt,inttim] = adcp_filt_sub(data,u_interp',v_interp',1:length(Z),40);
+[uintfilt,vintfilt,inttim,utid_baro,vtid_baro] = adcp_filt_sub(data,u_interp',v_interp',1:length(Z),40);
 saveas(figure(5),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','data_raw_filt_subsampled_1'],'fig')
 saveas(figure(6),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','data_raw_filt_subsampled_2'],'fig')
 
-% Save interpolated data
-bin_start = 1; % bin indice where good interpolated data for the whole dataset start
-bin_end   = length(Z);
-data.uintfilt=uintfilt(bin_start:bin_end,:);
-data.vintfilt=vintfilt(bin_start:bin_end,:);
-data.Z = Z(bin_start:bin_end);
-data.inttim = inttim;
+%% Save interpolated data
+bin_start           = input('Determine first bin indice with good interpolated data: '); % bin indice where good interpolated data for the whole dataset start
+bin_end             = input('Determine last bin indice with good interpolated data: '); % bin indice where good interpolated data for the whole dataset start
+data.uintfilt       = uintfilt(bin_start:bin_end,:);
+data.vintfilt       = vintfilt(bin_start:bin_end,:);
+data.Z              = Z(bin_start:bin_end);
+data.inttim         = inttim;
 save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '_int_filt_sub.mat'],'adcp','mooring','data','raw');
 
 %% Figure
-niv_u = (-1:0.05:1);
-niv_v = (-1:0.05:1);
+niv_u               = (-1:0.05:1);
+niv_v               = (-1:0.05:1);
 
 hf=figure('position', [0, 0, 1400, 1000]);
 %u
@@ -257,7 +268,7 @@ subplot(2,1,2);
 [C,h] = contourf(inttim,Z(bin_start:bin_end),vintfilt(bin_start:bin_end,:),niv_v); 
 set(h,'LineColor','none');
 caxis(niv_v([1 end]));
-h=colorbar;
+h     = colorbar;
 ylabel(h,'V [m s^-^1]');
 set(gca,'ydir', 'reverse');
 ylabel('Depth (m)');
@@ -268,25 +279,59 @@ title({[mooring.name, ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz']});
 
 graph_name = [fpath_output, mooring.name '_U_V_int_filt_sub'];
 set(hf,'Units','Inches');
-pos = get(hf,'Position');
+pos        = get(hf,'Position');
 set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
 print(hf,graph_name,'-dpdf','-r300');
- 
+
+% %% Plot tide
+% hf=figure('position', [0, 0, 1400, 1000]);
+% niv_tide = (-5:0.2:5);
+% utid_baro = utid_baro * 100;
+% vtid_baro = vtid_baro * 100;
+% %u
+% subplot(2,1,1);
+% colormap jet
+% [C,h] = contourf(inttim,Z(bin_start:bin_end),utid_baro(bin_start:bin_end,:),niv_tide); 
+% set(h,'LineColor','none');
+% caxis(niv_tide([1 end]));
+% h=colorbar;
+% ylabel(h,'U [cm s^-^1]');
+% set(gca,'ydir', 'reverse');
+% ylabel('Depth (m)');
+% ylim([0,adcp.instr_depth]);
+% %change figure label in HH:MM
+% gregtick;
+% title({[mooring.name, ' - ZONAL TIDE VELOCITY - RDI ',num2str(freq),' kHz']});
+% 
+% %v
+% subplot(2,1,2);
+% [C,h] = contourf(inttim,Z(bin_start:bin_end),vtid_baro(bin_start:bin_end,:),niv_tide); 
+% set(h,'LineColor','none');
+% caxis(niv_tide([1 end]));
+% h     = colorbar;
+% ylabel(h,'V [cm s^-^1]');
+% set(gca,'ydir', 'reverse');
+% ylabel('Depth (m)');
+% ylim([0,adcp.instr_depth]);
+% %change figure label in HH:MM
+% gregtick;
+% title({[mooring.name, ' - MERIDIONAL TIDE VELOCITY - RDI ',num2str(freq),' kHz']});
+
 %%  Write netcdf file           
 [yr_start , ~, ~] = gregorian(inttim(1));
-[yr_end,  ~, ~] = gregorian(inttim(length(inttim)));
+[yr_end,  ~, ~]   = gregorian(inttim(length(inttim)));
 
-ncid=netcdf.create([fpath_output,'ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'_1d.nc'],'NC_WRITE');
+ncid     = netcdf.create([fpath_output,'ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'_1d.nc'],'NC_WRITE');
  
 %create dimension
-dimidt = netcdf.defDim(ncid,'time',length(inttim));
-dimidz = netcdf.defDim(ncid,'depth',length(Z));
+dimidt   = netcdf.defDim(ncid,'time',length(inttim));
+dimidz   = netcdf.defDim(ncid,'depth',length(Z));
 %Define IDs for the dimension variables (pressure,time,latitude,...)
-time_ID=netcdf.defVar(ncid,'time','double',dimidt);
-depth_ID=netcdf.defVar(ncid,'depth','double',dimidz);
+time_ID  = netcdf.defVar(ncid,'time','double',dimidt);
+depth_ID = netcdf.defVar(ncid,'depth','double',dimidz);
 %Define the main variable ()
-u_ID = netcdf.defVar(ncid,'u','double',[dimidt dimidz]);
-v_ID = netcdf.defVar(ncid,'v','double',[dimidt dimidz]);
+u_ID     = netcdf.defVar(ncid,'u','double',[dimidt dimidz]);
+v_ID     = netcdf.defVar(ncid,'v','double',[dimidt dimidz]);
 %We are done defining the NetCdf
 netcdf.endDef(ncid);
 %Then store the dimension variables in
@@ -297,9 +342,6 @@ netcdf.putVar(ncid,u_ID,uintfilt);
 netcdf.putVar(ncid,v_ID,vintfilt);
 %We're done, close the netcdf
 netcdf.close(ncid);
+% -------------------------------------------------------------------------------------------
 
-% rmpath
-% rmpath('..\moored_adcp_proc');
-clear all; close all;
-
-% -------------------------------------------------------------------------------------------
\ No newline at end of file
+clear all; close all;
\ No newline at end of file
-- 
GitLab