diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..8a8a5431ace1b0f9339b36922d335c6825d10afd --- /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 0000000000000000000000000000000000000000..49b5f7a0a76844d9422cfb76947517d74dbdbf25 --- /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 0000000000000000000000000000000000000000..2a4beda03225682c64218bf1c30baad3ad89e14e --- /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 0000000000000000000000000000000000000000..12006812cdfad644a9f8036345e67a38a4b5cd61 --- /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 0000000000000000000000000000000000000000..4eb546bd9d5493802ecad3c61cf00be23a0f5ca5 --- /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 0000000000000000000000000000000000000000..d4e2acd2d025342071f055a67e3ab10618115a4e --- /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 0000000000000000000000000000000000000000..3a8c5ab44646fb81178fb738c166d492ec327e8d --- /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 0000000000000000000000000000000000000000..d3d616dd96dbd1cd8257394780580736b2e0cc68 --- /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 0000000000000000000000000000000000000000..b4123387158f8650a25ff10161fb2876f2563b44 --- /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 0000000000000000000000000000000000000000..fadfd79c9a231461c17b370eed97cd2fec7a44c2 --- /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 0000000000000000000000000000000000000000..3dbe693e936b33806a1f530081dea08a870033dc --- /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 0000000000000000000000000000000000000000..619a3a64a1496370cc3fb9282eb33b0b3179d8d6 --- /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 0000000000000000000000000000000000000000..654e0508e5f48b0660f6b8962c708c551f601077 --- /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 0000000000000000000000000000000000000000..8f5eebf33125c5b607d92f46109085087096b343 --- /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 0000000000000000000000000000000000000000..512be4c741fed006d62b977f9932ced006d1d42f --- /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 0000000000000000000000000000000000000000..e7d3e96747e18b127a1e78cfaef8854de3137820 --- /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 0000000000000000000000000000000000000000..bd1d56126a9b84c0e7b0d8d10d82d71e4929b9c8 --- /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 3f193220b36e9fd034551798e83a3d1d994042ca..2dec06fce0cd664a1449f2bb4f0e9cbffff450f4 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 310e7c7c5f716f7b279e0213fa5773a0ac8c64bd..ae84377c54728a04b0fef8383796816b6fd9cdea 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 dc5ba497b465ce1e042066a476a9d4f47127b21a..3fa843267553edee52fac02981afac239824a4ee 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 560fe2c3bf130a38b99291d6303882b3748ea5c1..f0dc31a72b676f7b7e9cef4f160c38c3f7c7590f 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 d4e2acd2d025342071f055a67e3ab10618115a4e..6b09adacb59669dbb3c184c1076fd205d3ed48b2 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 d68e9466925c6e828540fd855d09297401dca0b4..205e94abb6aa959589a9711096c45816770f65db 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