Skip to content
Snippets Groups Projects
Commit c6bd6abf authored by pierre.rousselot_ird.fr's avatar pierre.rousselot_ird.fr
Browse files

Cleaned version

parent 81776aa2
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 3482 additions and 54 deletions
.gitignore 0 → 100644
+ 4
− 0
View file @ c6bd6abf
*.000
*.pdf
*.csv
*.mat
filter/clean_median.m 0 → 100644
+ 201
− 0
View file @ c6bd6abf
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);
filter/diurne.m 0 → 100644
+ 176
− 0
View file @ c6bd6abf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
filter/fft_filter.m 0 → 100644
+ 80
− 0
View file @ c6bd6abf
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
filter/filterfft.m 0 → 100644
+ 71
− 0
View file @ c6bd6abf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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)')
filter/filtfilt.m 0 → 100644
+ 310
− 0
View file @ c6bd6abf
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.
filter/hamm.m 0 → 100644
+ 14
− 0
View file @ c6bd6abf
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
filter/hammfilter_nodec.m 0 → 100644
+ 29
− 0
View file @ c6bd6abf
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;
filter/hamming.m 0 → 100644
+ 36
− 0
View file @ c6bd6abf
% 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;
filter/lisse.m 0 → 100644
+ 131
− 0
View file @ c6bd6abf
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
filter/medfilt1_perso.m 0 → 100644
+ 22
− 0
View file @ c6bd6abf
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
filter/median_filter.m 0 → 100644
+ 106
− 0
View file @ c6bd6abf
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
filter/mfilter.m 0 → 100644
+ 1075
− 0
View file @ c6bd6abf
This diff is collapsed.
filter/naninterp.m 0 → 100644
+ 5
− 0
View file @ c6bd6abf
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
filter/nanmedfilt2.m 0 → 100644
+ 44
− 0
View file @ c6bd6abf
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
filter/ndnanfilter.m 0 → 100644
+ 484
− 0
View file @ c6bd6abf
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
filter/spikeRemoval.m 0 → 100644
+ 601
− 0
View file @ c6bd6abf
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
moored_adcp_proc/adcp_check_surface.m
+ 6
− 6
View file @ c6bd6abf
...@@ -10,7 +10,7 @@ set(gca,'YTick',bins); ...@@ -10,7 +10,7 @@ set(gca,'YTick',bins);
ylabel('Bins'); ylabel('Bins');
hold on; hold on;
gregtick; gregtick;
title('Meridional velocity before correction'); title('Zonal velocity before correction');
subplot 323; subplot 323;
imagesc(time,bin(bins),u1(bins,:),[-1 1]); imagesc(time,bin(bins),u1(bins,:),[-1 1]);
...@@ -18,7 +18,7 @@ set(gca,'YTick',bins); ...@@ -18,7 +18,7 @@ set(gca,'YTick',bins);
hold on; hold on;
ylabel('Bins'); ylabel('Bins');
gregtick; gregtick;
title('Meridional velocity after correction'); title('Zonal velocity after correction');
subplot 322; subplot 322;
imagesc(time,bin(bins),v(bins,:),[-1 1]); imagesc(time,bin(bins),v(bins,:),[-1 1]);
...@@ -26,7 +26,7 @@ set(gca,'YTick',bins); ...@@ -26,7 +26,7 @@ set(gca,'YTick',bins);
hold on; hold on;
ylabel('Bins'); ylabel('Bins');
gregtick; gregtick;
title('Zonal velocity before correction'); title('Meridional velocity before correction');
subplot 324; subplot 324;
imagesc(time,bin(bins),v1(bins,:),[-1 1]); imagesc(time,bin(bins),v1(bins,:),[-1 1]);
...@@ -34,10 +34,10 @@ set(gca,'YTick',bins); ...@@ -34,10 +34,10 @@ set(gca,'YTick',bins);
hold on; hold on;
ylabel('Bins'); ylabel('Bins');
gregtick; gregtick;
title('Zonal velocity after correction'); title('Meridional velocity after correction');
subplot 325; subplot 325;
for i=1:length(bins); for i=1:length(bins)
plot(time,hammfilter_nodec(z(bins(i),:),73),'color','k'); plot(time,hammfilter_nodec(z(bins(i),:),73),'color','k');
hold on; hold on;
end end
...@@ -47,7 +47,7 @@ ylabel('Bins'); ...@@ -47,7 +47,7 @@ ylabel('Bins');
title('Bin depth after correction'); title('Bin depth after correction');
subplot 326; subplot 326;
for i=1:length(bins); for i=1:length(bins)
plot(time,hammfilter_nodec(z(bins(i),:),73),'color','k'); plot(time,hammfilter_nodec(z(bins(i),:),73),'color','k');
hold on; hold on;
end end
......
moored_adcp_proc/adcp_filt_sub.m
+ 83
− 46
View file @ c6bd6abf
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) for iidep = 1:length(intdepvec)
adcpui = ui(iidep,:); adcpui = ui(iidep,:);
adcpvi = vi(iidep,:); adcpvi = vi(iidep,:);
valid = find(~isnan(adcpui)); valid = find(~isnan(adcpui));
invalid = 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 % Horizontal interpolation to fill gaps
timval = data.time(valid); timval = data.time(valid);
adcpui = interp1(timval,adcpui(valid),data.time); adcpui = interp1(timval,adcpui(valid),data.time);
adcpvi = interp1(timval,adcpvi(valid),data.time); adcpvi = interp1(timval,adcpvi(valid),data.time);
sum(isnan(adcpui)); sum(isnan(adcpui));
sum(isnan(adcpvi)); 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); %% Hamming filter
sum(isnan(adcpui)); adcpui(valid(1):valid(end)) = mfilter(adcpui(valid(1):valid(end)),sf,1/nhours,0,2*nhours,1);
sum(isnan(adcpvi)); adcpvi(valid(1):valid(end)) = mfilter(adcpvi(valid(1):valid(end)),sf,1/nhours,0,2*nhours,1);
% adcpui(invalid) = NaN;
% adcpvi(invalid) = NaN; %% Fourier filter
sum(isnan(adcpui)); % [adcpui(valid(1):valid(end))] = fft_filter(adcpui(valid(1):valid(end)),2,[24 Inf]);
uifilt(iidep,1:length(adcpui)) = adcpui; % [adcpvi(valid(1):valid(end))] = fft_filter(adcpvi(valid(1):valid(end)),2,[24 Inf]);
vifilt(iidep,1:length(adcpvi)) = adcpvi;
% subsampling
inttim = [ceil(data.time(1)):0.25:floor(data.time(end))]; %% 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); uintfilt(iidep,1:length(inttim)) = interp1(data.time,transpose(uifilt(iidep,:)),inttim);
vintfilt(iidep,1:length(inttim)) = interp1(data.time,transpose(vifilt(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 elseif length(valid)<=1 || length(valid(1):valid(end))<=240
% subsampling
uifilt(iidep,1:length(adcpui)) = nan; %% subsampling
vifilt(iidep,1:length(adcpvi)) = nan; uifilt(iidep,1:length(adcpui)) = NaN;
inttim = [ceil(data.time(1)):0.25:floor(data.time(end))]; vifilt(iidep,1:length(adcpvi)) = NaN;
uintfilt(iidep,1:length(inttim)) =nan; inttim = [ceil(data.time(1)):0.25:floor(data.time(end))];
vintfilt(iidep,1:length(inttim)) =nan; uintfilt(iidep,1:length(inttim)) = NaN;
vintfilt(iidep,1:length(inttim)) = NaN;
else else
% Horizontal interpolation to fill gaps
timval = data.time(valid); %% Horizontal interpolation to fill gaps
adcpui = interp1(timval,adcpui(valid),data.time); timval = data.time(valid);
adcpvi = interp1(timval,adcpvi(valid),data.time); adcpui = interp1(timval,adcpui(valid),data.time);
% filtering adcpvi = interp1(timval,adcpvi(valid),data.time);
adcpui = mfilter(adcpui,sf,1/nhours,0,2*nhours,1);
adcpvi = mfilter(adcpvi,sf,1/nhours,0,2*nhours,1); %% Hamming filter
% adcpui(invalid) = NaN; adcpui = mfilter(adcpui,sf,1/nhours,0,2*nhours,1);
% adcpvi(invalid) = NaN; adcpvi = mfilter(adcpvi,sf,1/nhours,0,2*nhours,1);
uifilt(iidep,1:length(adcpui)) = adcpui;
vifilt(iidep,1:length(adcpvi)) = adcpvi; %% Subsampling
%subsampling uifilt(iidep,1:length(adcpui)) = adcpui;
inttim = [ceil(data.time(1)):0.25:floor(data.time(end))]; 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); uintfilt(iidep,1:length(inttim)) = interp1(timval,transpose(uifilt(iidep,valid)),inttim);
vintfilt(iidep,1:length(inttim)) = interp1(timval,transpose(vifilt(iidep,valid)),inttim); vintfilt(iidep,1:length(inttim)) = interp1(timval,transpose(vifilt(iidep,valid)),inttim);
inttim = inttim; inttim = inttim;
end end
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]); hf=figure('position', [0, 0, 1400, 1000]);
subplot 121; 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]); %set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10]);
gregtick; gregtick;
ylabel('Bins'); ylabel('Bins');
title('U field - data raw'); title('U field - data raw');
subplot 122; subplot 122;
imagesc(data.time,intdepvec,uintfilt,[-1 1]); imagesc(inttim,intdepvec,uintfilt,[-1 1]);
%set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10]); %set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10]);
gregtick; gregtick;
ylabel('Bins'); ylabel('Bins');
...@@ -79,14 +114,16 @@ title('U field - data interpolated, filtered and subsampled'); ...@@ -79,14 +114,16 @@ title('U field - data interpolated, filtered and subsampled');
hf=figure('position', [0, 0, 1400, 1000]); hf=figure('position', [0, 0, 1400, 1000]);
subplot 121; 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]) %set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10])
gregtick; gregtick;
ylabel('Bins'); ylabel('Bins');
title('V field - data raw'); title('V field - data raw');
subplot 122; subplot 122;
imagesc(data.time,intdepvec,vintfilt,[-1 1]); imagesc(inttim,intdepvec,vintfilt,[-1 1]);
%set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10]); %set(gca,'YLim',[min(intdepvec)-10 max(intdepvec)+10]);
gregtick; gregtick;
ylabel('Bins'); ylabel('Bins');
title('V field - data interpolated, filtered and subsampled'); title('V field - data interpolated, filtered and subsampled');
end
moored_adcp_proc/adcp_shadowzone.m
+ 4
− 2
View file @ c6bd6abf
function z_shadowzone=adcp_shadowzone(adcp_dpt,beamangle); function z_shadowzone=adcp_shadowzone(adcp_dpt,beamangle)
% 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); ...@@ -8,4 +8,6 @@ function z_shadowzone=adcp_shadowzone(adcp_dpt,beamangle);
% %
% R. Kopte, 2014/11/20 % R. Kopte, 2014/11/20
z_shadowzone=adcp_dpt*(1-cosd(beamangle)); z_shadowzone = adcp_dpt*(1-cosd(beamangle));
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment