%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Hybrid despiking method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Amir Hossein Birjandi June 2010 % University of Manitoba % Department of mechanical and Manufacturing Engineering %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The present code removes high amplitude spikes form the Acoustic Dopple Velocimetry (ADV) % signal or any other signals and replace them by a linear interpolation between two contiguous % valid data points. The advantage of this method is its accuracy in dealing with highly contaminated signals. % This method removes high amplitude spikes in two steps. In the first step spikes with high a deviation from the mean % value are removed by the thresholding and acceleration mothods. In the % second step spikes and abnormal data points with a low deviation from the % mean vale are explored. Removed data points are replaced by a linear % interpolation between two contiguous valid data points. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc clear %Read the input data file %Insert your file name here [data]=textread('Raw ADV dataset.dat'); fs=200; % Sampling frequency t=data(:,1); %Time series tb=1; % Begin time index te=round(length(t)); % End time index %Velocity vectors U_x=data(tb:te,3); U_y=data(tb:te,4); U_z1=data(tb:te,5); U_z2=data(tb:te,6); U_z=(U_z1+U_z2)/2; %Fluctuating velocities u_x=U_x-mean(U_x); u_y=U_y-mean(U_y); u_z=U_z-mean(U_z); %Plot the original signal U1=U_x; figure(1) plot(t,U1) %%%%%%%%%%%%%%%%% %Signal segmentation %%%%%%%%%%%%%%%%% % If you need a segment of the signal for the process you can select your % segment in this section. You should plug in the start and end point of % the segment. %%%%%%%%%%%%%%%%% s_segment= input('Do you want select a segment of your signal? y/n [y]: ', 's'); close figure 1; if isempty(s_segment) s_segment= 'y'; end % Segmentation if s_segment=='y' satisfaction='n'; while satisfaction=='n' t=data(:,1); tb=input('Please enter the beginning time: '); %Start time te=input('Please enter the end time: '); % End time tb=round((tb-t(1))*fs+1);%Index of the start time te=round((te-t(1))*fs+1);%Index of the end time t=t(tb:te); U1=data(tb:te,3); figure(2) plot(t,U1) satisfaction=input('Are you satisfied with this selection? y/n [y]: ','s'); if isempty(satisfaction) satisfaction = 'y'; end end close figure 2; end % Despiking on all 3 velocity components for i=1:1:3 % Velocity vectors (u, v, w) U1=data(tb:te,2+i); % Fluctuations u1=U1-mean(U1); u_prime=u1; tshift=t-t(1); % shifted time series (starts from zero) figure(1) plot(tshift,u1,'b*') hold on % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Hybrid despiking method % If your signal contains high amplitude spikes it may affect your % statistics. So before any further calculation you require to remove % spikes from your signal. For that perpose we offer the Hybrid despiiking method here. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% despiking = input('Do you want to despike your signal? y/n [y]: ', 's'); if isempty(despiking) despiking = 'y'; end if despiking=='y' uxnew1=u_prime; % Stores the valid data points u_spike=zeros(length(u_prime),1);% Stores the spikes % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Threshold filtering % This method filtes high amplitude spikes % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T_filter = input('Do you want to use the thresholding filtering method first? y/n [y]: ', 's'); if isempty(T_filter) T_filter = 'y'; end if T_filter=='y' V_threshold = input('Please enter the the velocity threshold [3*std]: '); if isempty(V_threshold) V_threshold = 3*std(u1); % defult for velocity threshold is 3*std end u_spike=((abs(u_prime)>=V_threshold)) .*u1+u_spike;% Filters points beyond the threshold u_prime=(abs(u_prime)(0.0001*length(u_prime))) niter=niter+1; n_spikes=length(uxnew1); u_x2=circshift(u_prime,[-1,0]); du_x=(u_x2-u_prime).*fs; %Acceleration (first derivative of the signal) threshold=acth*std(uxnew1); %Filtering threshold %%%%%%%%%%%%%%%%%%%% %Despiking criteria %%%%%%%%%%%%%%%%%%%% u_spike=((abs(u_prime)>=threshold & abs(du_x)>=(acfa*9.81))) .*u1+u_spike; u_prime=(abs(u_prime)(0.0001*length(u_prime))) n_spikes=length(uxnew1); %First derivative u_x2=circshift(u_prime,[-1,0]); u_x3=circshift(u_prime,[1,0]); du_x=(u_x2-u_x3)/2; %First derivative %Second derivative du_x2=circshift(du_x,[-1,0]); du_x3=circshift(du_x,[1,0]); ddu_x=(du_x2-du_x3)/2; %Second derivative lambda=sqrt(2*log(length(u_prime))); % Universal threshold %Principal axes and rotation angle teta=atan(sum((u_prime-mean(u_prime)).*ddu_x)/sum((u_prime-mean(u_prime)).^2)); % Rotation angle e1=lambda*std(u_prime); % Primcipal axis in u-du plane e2=lambda*std(du_x); % Primcipal axis in u-du and du-ddu planes e3=lambda*std(ddu_x); % Primcipal axis in du-ddu plane e4=sqrt(((lambda*std(u_prime)*cos(teta))^2-(lambda*std(ddu_x)*sin(teta))^2)/((cos(teta))^4-(sin(teta))^4)); % Primcipal axis in u-ddu plane e5=sqrt(((lambda*std(ddu_x)*cos(teta))^2-(lambda*std(u_prime)*sin(teta))^2)/((cos(teta))^4-(sin(teta))^4)); % Primcipal axis in u-ddu plane u_x22=u_prime.*cos(teta)+ddu_x.*sin(teta); % Rotates the data by theta ddu_x22=ddu_x.*cos(teta)-u_prime.*sin(teta); % Rotates the data by theta % %If you want to see the filtering ellipsoids you can activate this % %section % %Plot the ellipses and filter the outside points % %%%%%%%%%%%% % %u-du plane % %%%%%%%%%%%% % [r1,c1,filter1]=find((((u_prime-mean(u_prime)).^2)/e1^2+((du_x).^2)/e2^2-1)>0); % detect spikes % %Filtering ellipse % xc = mean(u_prime); yc = mean(du_x); % tt = 0:pi/200:2*pi; % xu = e1*cos(tt); % yv = e2*sin(tt); % figure(5) % plot(xu+xc,yv+yc) % plot the filtering ellipse % hold on % plot(u_prime,du_x,'b*') % plot the whole dataset % hold on % plot(u_prime(r1),du_x(r1),'rO') % plot the spikes % %%%%%%%%%%%% % %du-ddu plane % %%%%%%%%%%%% % [r2,c2,filter2]=find(((du_x).^2/e2^2+(ddu_x).^2/e3^2-1)>0); % detect spikes % %Filtering ellipse % xc = 0; yc = 0; % tt = 0:pi/200:2*pi; % xu = e2*cos(tt); % yv = e3*sin(tt); % figure(6) % plot(xu,yv)% plot the filtering ellipse % hold on % plot(du_x,ddu_x,'b*') % plot the whole dataset % hold on % plot(du_x(r2),ddu_x(r2),'rO')% plot the spikes % % %%%%%%%%%%%% % %u-ddu plane % %%%%%%%%%%%% % [r3,c3,filter3]=find(((u_x22).^2/e4^2+(ddu_x22).^2/e5^2-1)>0); % detect spikes % %Filtering ellipse % xc = 0; yc = 0; % tt = 0:pi/200:2*pi; % xu = e4*cos(tt); % yv = e5*sin(tt); % xx = xu*cos(teta) - yv*sin(teta); % yy = xu*sin(teta) + yv*cos(teta); % x = xc + xx; % y = yc + yy; % figure(7) % plot(x,y)% plot the filtering ellipse % hold on % plot(u_prime,ddu_x,'b*') % plot the whole dataset % hold on % plot(u_prime(r3),ddu_x(r3),'rO')% plot the spikes % %%%%%%%%%%%%%%%%%%%%%% % Filtering criteria % %%%%%%%%%%%%%%%%%%%%%% u_spike=(((u_prime).^2/e1^2+(du_x).^2/e2^2-1)>0 | ((du_x).^2/e2^2+(ddu_x).^2/e3^2-1)>0 | ((u_x22).^2/e4^2+(ddu_x22).^2/e5^2-1)>0).*u1+u_spike; u_prime=(((u_prime).^2/e1^2+(du_x).^2/e2^2-1)<0 & ((du_x).^2/e2^2+(ddu_x).^2/e3^2-1)<0 & ((u_x22).^2/e4^2+(ddu_x22).^2/e5^2-1)<0).*u_prime; [r,c,uxnew1]=find(u_prime); u_prime=u_prime-(u_prime ~=0).*mean(uxnew1); % Removes the mean value [r,c,uxnew1]=find(u_prime); % uxnew1 is a vector of clean data points n_spikes=n_spikes-length(uxnew1); end [r2,c2,uxnew2]=find((u_prime ~=0).*(u1-u_prime)); s_dif=mean(uxnew2); if u_prime(length(u1))==0 u_prime(length(u1))=u1(length(u1)); end if u_prime(1)==0 u_prime(1)=u1(1); end [r,c,uxnew1]=find(u_prime); %Valid data points vector uxnewi=r(1):1:r(length(r)); %Evenly distributed vector from 1 to the length of the signal uxnewj=interp1(r,uxnew1,uxnewi);% Linear interpolation to replace spikes u_prime=transpose(uxnewj)+s_dif; %Filtered data points [r3,c3,u_spike]=find(u_spike); u_spike=transpose(u_spike); %Spikes vector figure(4) plot(tshift,u1,'b*') hold on plot(tshift(r3),u1(r3),'rO') hold on plot(tshift(r3),u_prime(r3),'r*') figure(5) plot(tshift,u_prime,'r*') hold on plot(tshift,u1,'b*') end %%%%%%%%%%%%%%%%%%% %Write the new value %%%%%%%%%%%%%%%%%%% if i==1 u_x=u1; U_x=U1; u_x_new=u_prime; U_x_new=u_x_new+mean(U1); figure(6) plot(tshift,U_x,'b*') hold on plot(tshift,U_x_new,'r*') end if i==2 u_y=u1; U_y=U1; u_y_new=u_prime; U_y_new=u_y_new+mean(U1); figure(6) plot(tshift,U_y,'b*') hold on plot(tshift,U_y_new,'r*') end if i==3 u_z=u1; U_z=U1; u_z_new=u_prime; U_z_new=u_z_new+mean(U1); figure(6) plot(tshift,U_z,'b*') hold on plot(tshift,U_z_new,'r*') end ckey=input('Press any key to continue ','s'); close all end fidout = fopen('Despiked ADV dataset.dat','w'); fprintf(fidout, 'variables= "time","U_x","U_y","U_z","U_x_new","U_y_new","U_z_new" \r\n'); % fprintf(fidout, '%f %f %f %f %f %f %f % \r\n',transpose(tshift),transpose(U_x),transpose(U_y),transpose(U_z),transpose(U_x_new),transpose(U_y_new),transpose(U_z_new)); new_vector=[transpose(tshift);transpose(U_x);transpose(U_y);transpose(U_z);transpose(U_x_new);transpose(U_y_new);transpose(U_z_new)]; fprintf(fidout, '%f %f %f %f %f %f %f \r\n',new_vector); status = fclose('all');