diff --git a/calculate_v5.m b/calculate_v5.m index 560772f..c62ea5d 100644 --- a/calculate_v5.m +++ b/calculate_v5.m @@ -1,367 +1,371 @@ function [dh_max, dh_min, dh_f, dh_pp, r_f,fdom]= calculate_v5(folder_save,freq,time,nb,nbsensor,dat_cp,hm,length_sig,lowp,highp,removefrequ,percentage,Em,rho,nm,folder,del_pump, F_pump, F_pump2,harmo1, harmo2, harmo3,log,norm, MP, projectname,dele,channel,save_post,folder_res,rel) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%Function defintition %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %This function is called by the function main after calling the read %function. It receives a matrix dat_cp with all the pressure fluctuation %data already centered % This is the normal version % I am modifying this version + +%Hello tanguy + + %% Here is what the function do : % - Add low pass filter to the temporal signal % - Add a high pass filter to the temporal signal % - Remove unwanted frequencies to the temporal signal % - Apply fft with proper windowing % - Delete the extremums in function of the confidence interval chosen. % - Plot and store the plot in the time domain. % - Plot and store the plot in the frequency domain. % - Calculate the data such as DH/H_max, DH/H_min DH/H_pp etc... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% filter the data %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 2.1 : add low pass filer lowp_exist = lowp>1; % Check if a lowpass filter if lowp_exist ==1 % if a lowpass is requested dlp = designfilt('lowpassiir','FilterOrder',7, ... % Definition of the lowpass filter 'PassbandFrequency',lowp,'PassbandRipple',0.2, ... 'SampleRate',freq); dat_cp_lp = filtfilt(dlp,dat_cp); % Application of the lowpass filter on the data. else dat_cp_lp=dat_cp; % data centered, permuted and "lowpassed" end %% 2.2 add high pass filer highp_exist = highp>0; % Check if a highpass filter was demanded in GUI if highp_exist ==1 % if a highpass is requested dhp = designfilt('highpassiir','FilterOrder',7, ... % Definition of the highpass filter 'PassbandFrequency',highp,'PassbandRipple',0.2, ... 'SampleRate',freq); dat_cp_lp_hp = filtfilt(dhp,dat_cp_lp); % Application of the highpass filter on the data. else dat_cp_lp_hp=dat_cp_lp; % data centered, permuted "lowpassed" and "highpassed" end %% 2.3 : Remove specific frequencies (50 Hz, 100 hz etc...) d1=split(removefrequ,' ') delfrequ=str2double(d1); % removefrequ = the frequencies that the user ask to remove in the GUI dat_cp_lp_hp_f=dat_cp_lp_hp; % Predifine the variable if delfrequ>1 % if we want to delete specified frequencies for i=1:length(delfrequ) % The Filterorder is very important ! If set too small too much data will be removed d(i) = designfilt('bandstopiir','FilterOrder',4, ... % definition of the filter, order 4. From delfrequ-0.5 Hz to delfrequ + 0.5 Hz 'HalfPowerFrequency1',delfrequ(i)-0.5,'HalfPowerFrequency2',delfrequ(i)+0.5, ... 'DesignMethod','butter','SampleRate',freq); dat_cp_lp_hp_f = filtfilt(d(i),dat_cp_lp_hp_f); % Application of the filter if harmo1>0 % If we want to also delete the harmonic of the specified frequency. d1(i) = designfilt('bandstopiir','FilterOrder',4, ... % definition of the filter, order 4. From delfrequ-1 Hz to delfrequ + 1 Hz 'HalfPowerFrequency1',2*(delfrequ(i)-0.5),'HalfPowerFrequency2',2*(delfrequ(i)+0.5), ... 'DesignMethod','butter','SampleRate',freq); dat_cp_lp_hp_f = filtfilt(d1(i),dat_cp_lp_hp_f); % Application of the filter end if harmo2>0 d2(i) = designfilt('bandstopiir','FilterOrder',4, ... % definition of the filter, order 4. From delfrequ-1.5 Hz to delfrequ + 1.5 Hz 'HalfPowerFrequency1',3*(delfrequ(i)-0.5),'HalfPowerFrequency2',3*(delfrequ(i)+0.5), ... 'DesignMethod','butter','SampleRate',freq); dat_cp_lp_hp_f = filtfilt(d2(i),dat_cp_lp_hp_f); % Application of the filter end if harmo3>0 d3(i) = designfilt('bandstopiir','FilterOrder',4, ... % definition of the filter, order 4. From delfrequ-2 Hz to delfrequ + 2 Hz 'HalfPowerFrequency1',4*(delfrequ(i)-0.5),'HalfPowerFrequency2',4*(delfrequ(i)+0.5), ... 'DesignMethod','butter','SampleRate',freq); dat_cp_lp_hp_f = filtfilt(d3(i),dat_cp_lp_hp_f); % Application of the filter end end end save([folder_save 'calc_p']) %% 2.4 :Remove circuit pump frequency dat_cp_lp_hp_f_pump=dat_cp_lp_hp_f; if del_pump>0 % if we want to delete specified frequencies. Defined in the GUI. dat_cp_lp_hp_f_pump1=dat_cp_lp_hp_f; for i=1:nb d_pump(i) = designfilt('bandstopiir','FilterOrder',4, ... % define a filter to delete the pump frequency. 'HalfPowerFrequency1',abs(F_pump(i)-0.4),'HalfPowerFrequency2',abs(F_pump(i)+0.4), ... 'DesignMethod','butter','SampleRate',freq); d_pump6(i) = designfilt('bandstopiir','FilterOrder',4, ... % define a filter to delete the 6th pump frequency. 'HalfPowerFrequency1',abs(6*F_pump(i)-0.4),'HalfPowerFrequency2',abs(6*F_pump(i)+0.4), ... 'DesignMethod','butter','SampleRate',freq); % delete the sixth harmonic (6 blades) if F_pump2 ~=0 %There are two pumps on the the third platform. If we are on the third platform we have to delete the frequency of the second pump d_pump2(i) = designfilt('bandstopiir','FilterOrder',4, ... % define the filter of the second pump 'HalfPowerFrequency1',abs(F_pump2(i)-0.4),'HalfPowerFrequency2',abs(F_pump2(i)+0.4), ... 'DesignMethod','butter','SampleRate',freq); d_pump26(i) = designfilt('bandstopiir','FilterOrder',4, ... % define a filter to delete the 6th pump frequency of the second pump. 'HalfPowerFrequency1',abs(6*F_pump2(i)-0.4),'HalfPowerFrequency2',abs(6*F_pump2(i)+0.4), ... 'DesignMethod','butter','SampleRate',freq); end dat_cp_lp_hp_f_cut=squeeze(dat_cp_lp_hp_f_pump1(:,i,:)); dat_cp_lp_hp_f_pump_cut = filtfilt(d_pump(i),dat_cp_lp_hp_f_cut); % remove the phase shift dat_cp_lp_hp_f_pump6_cut = filtfilt(d_pump6(i),dat_cp_lp_hp_f_pump_cut);% remove the phase shift of the 6th harmonic if F_pump2 ~=0 dat_cp_lp_hp_f_pump2_cut = filtfilt(d_pump2(i),dat_cp_lp_hp_f_pump6_cut); dat_cp_lp_hp_f_pump26_cut = filtfilt(d_pump26(i),dat_cp_lp_hp_f_pump2_cut); dat_cp_lp_hp_f_pump(:,i,:)=dat_cp_lp_hp_f_pump26_cut; else dat_cp_lp_hp_f_pump(:,i,:)=dat_cp_lp_hp_f_pump6_cut; % data centered, permuted, lowpassed, highpassed, filtered, pump_filtered. end end end save([folder_save 'calc_p2']) %% 2.5 :Apply windowing for proper spectrum (for dh_f only) %pwelch provides the PSD but not the amplitude of different frequences therefore we will not use windowing for df/f calc Method = 'Hamming'; fs=2000; Overlap=0.5; Window_Size=2^13; Time = [0:1/fs:20]; for i=1:nb for j=1:nbsensor Signal = dat_cp_lp_hp_f_pump(:,i,j); [fhf(:,i,j), pxxf(:,i,j)] = Sxx_Function_Local(Time, Signal, Method, Overlap, Window_Size); end end % pxxf=sqrt(pxxf); fhf_r=fhf(:,1,1); % fhf is always the same [dh_f1,ind_f]= max(pxxf); fdom= squeeze(round(fhf_r(ind_f),3)); r_f= abs((round(60*fdom./nm',3))); %% 2.6 :Take out xx % of the data (interval confidence) k=round((1-(percentage/100))*length_sig/2); % nb of data to remove. In fucntion of the chose percentage of confidence interval chosen in GUI for i=1:nb % nb of measurement for j=1:nbsensor % nb of sensor in each measurement temp=dat_cp_lp_hp_f_pump(:,i,j); [~, I(:,i,j)] = maxk(temp,k); % fing the k largest elements temp(I(:,i,j))=[]; % delete the largest elements dat_f1(:,i,j)=temp; end end for i=1:nb % nb of measurement for j=1:nbsensor % nb of sensor in each measuremen temp=dat_f1(:,i,j); [Bmin(:,i,j), Imin(:,i,j)] = mink(temp,k); % fing the k samllest elements temp(Imin(:,i,j))=[]; % delete the smallest elements dat_cp_f97(:,i,j)=temp; end end %% 2.7 : Calculate the key values : dh_max1=max(dat_cp_f97); %dat_cp_f97 = data centered permuted filtered with confidence interval 97% dh_min1=min(dat_cp_f97); %[dh_f1,ind_f]= max(ft); %% Convert the signal from bars to DH/H dh_max = squeeze(round(dh_max1.*100.*100000./(Em.*rho),10)); dh_min = squeeze(round(dh_min1.*100.*100000./(Em.*rho),10)); dh_f = squeeze(round(dh_f1.*100.*100000./(Em.*rho),3)); dh_pp=dh_max-dh_min; % if there is only one MP the matrix have one less dimension so we need to % do this trick to avoid an error if nb ==1 dh_max=dh_max'; dh_min =dh_min'; dh_f=dh_f'; dh_pp=dh_pp'; fdom=fdom'; r_f=r_f'; end save([folder_save 'calc']) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% Save postprocessed signals %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % We know have done all the signal analysis. We will save the post-treated % signals with the following lines. if save_post ==1 chan= channel{1,:}; for i = 1 :nb xlswrite([folder_res,'\Excel\signal_post\', projectname,'_MP',num2str(MP(i)),'_post_pro.xlsx'],chan,['MP', num2str(MP(i))],'A1') xlswrite([folder_res,'\Excel\signal_post\', projectname,'_MP',num2str(MP(i)),'_post_pro.xlsx'],squeeze(dat_cp_f97(:,i,:)),['MP', num2str(MP(i))],'A2') end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% PLOT TIME DOMAIN %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k=0; t= (1:freq*time)/freq; % time vector %% define the min and the max value of the ordinate axis. if nb==1 % if we have only one sensor the matrix have one less dimension it requires special attention ma=max(max(squeeze(dat_cp))*100000*100./(rho*Em')); mi=min(min(squeeze(dat_cp))*1000./hm'); else ma=max(squeeze(max(dat_cp)),[],2)*1000./hm'; % ma is the maximal DH/H value of each MP mi=min(squeeze(min(dat_cp)),[],2)*1000./hm'; % mi is the minimal DH/H value of each MP end timegraph =figure('Name','tempplot','visible','off'); %timegraph =figure('Name','tempplot') set(0, 'CurrentFigure', timegraph); for i=1:nb % Nb of MP files to treat for j= 1:nbsensor % nb of sensors per MP file k=k+1; % for the name of the plot set(0, 'CurrentFigure', timegraph); plot(downsample(t,10),downsample(squeeze(dat_cp(:,i,j)*1000./hm(i)),10),'Color',[1,0,0]); grid on grid minor set(gca,'Fontsize',18); ylab=ylabel('$\Delta H/H_{\rm{M}}\ (\%)$','Interpreter','latex'); ylim([1.1*mi(i) 1.1*ma(i)]); yticks([round(mi(i),1),0,round(ma(i),1)]); if mod(j,4)==0 || mod(j,nbsensor)==0 % If it is the last graph on the page we want to include xlabel set(gcf,'units','points','position',[10,10,550,245]) % as the xlabel is added we need to slightly change the shape of the graph xlab=xlabel('$t \ $(s)','Interpreter','latex'); set(xlab,'Fontsize',20); set(ylab,'Fontsize',20); set(ylab,'Fontweight','bold'); else set(gca, 'XTickLabel', []) set(gcf,'units','points','position',[10,10,550,213]) % set size of the plot end %namefig=[folder '\Results\images\' projectname '\' projectname '-timeplot-MP-' num2str(MP(i)) '-' num2str(j)]; namefig=[folder '\Results\images\' projectname '\' projectname '-timeplot-MP-' num2str(MP(i)) '-' num2str(j)]; set(0, 'CurrentFigure', timegraph); save2pdf(namefig); % save the plot end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% PLOT FREQ DOMAIN %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% TO save time, we discard the data after the lowpass % if lowp_exist ==1 % indice_lowp=find(fhf_r==lowp); % end % fhf_r(indice_lowp+100:end-1)=[]; % pxxf(indice_lowp+100:end-1,:,:)=[]; freqgraph=figure('Name','frequ','visible','off'); set(0, 'CurrentFigure', freqgraph); k=0; rel_decided=rel; for i=1:nb rel=rel_decided; %% if the rotational speed is less than 50 rpm, the x-plot is set to absolute. if abs(nm(i))<50 rel=0; end for j= 1:nbsensor k=k+1; maxf=max(squeeze(pxxf(:,i,j))); % max value for each sensor of each MP %plot(fr*60./abs(nm(i)),squeeze(ft(:,i,j)/maxf),'Color',[1,0,0]) % frequency domain is normalizedby maximal value if norm ==0 % In case we don't want to have normalized freqency plots. maxf= Em(i)*rho(i)*maxf./(maxf *100*100000); % If we don't want to normallize we want to have DH/H values. 100000 is to convert bars to Pa 100 is to convert to percent Em * rho (=rho *g*H) is to get the value of the Head in Pa ylab=ylabel('$\Delta H/{H} \ $(%)','Interpreter','latex'); end set(0, 'CurrentFigure', freqgraph); if log>0 % If we want a semilog plot if rel>0 % If we want f at the x-axis instead of f/f. ab=bar(downsample(fhf_r*60./abs(nm(i)),2),downsample(squeeze(pxxf(:,i,j)/maxf),2),'r'); set(gca,'XScale','log') else ab=bar(downsample(fhf_r*60,2),downsample(squeeze(pxxf(:,i,j)/maxf),2),'r'); set(gca,'XScale','log') end %ab=semilogx(fhf_r*60,squeeze(pxxf(:,j,i)/maxf),'Color',[1,0,0]) else if rel>0 % If we want f at the x-axis instead of f/f. plot(fhf_r*60./abs(nm(i)),squeeze(pxxf(:,i,j)/maxf),'r'); else plot(fhf_r,squeeze(pxxf(:,i,j)/maxf),'r'); end end grid on grid minor set(gca,'Fontsize',16); if norm ==0 ylab=ylabel('$\Delta H/{H} \ $(-)','Interpreter','latex'); else ylab=ylabel('$\Delta H/{\Delta H}_{\rm{f,dom}} \ $(-)','Interpreter','latex'); end if log>0 xlim([0.1,min(1.1*freq*60./abs(nm(i)),1.2*lowp*60./abs(nm(i)))]); %xlim([0.1,min((1.1*freq*60),1.2*lowp*60)]); elseif rel>0 xlim([0.1,min(1.1*freq*60./abs(nm(i)),1.2*lowp*60./abs(nm(i)))]); %xlim([0.1,min((1.1*freq*60),1.2*lowp*60)]); else xlim([0.1,min(1000,1.2*lowp)]); end % ylim([0 5]); if mod(j,4)==0 || mod(j,nbsensor)==0 % If it is the last graph of the page... set(gcf,'units','points','position',[10,10,550,250]); %set(gcf,'units','points','position',[10,10,110,49]); if rel>0 xlab=xlabel('${f}/{f_{\rm{n}}}\ $ (-) ','Interpreter','latex'); else xlab=xlabel('${f}$ (Hz) ','Interpreter','latex'); xlim([0.1,min(1000,1.2*lowp)]); end set(xlab,'Fontsize',25); set(ylab,'Fontsize',18); else %xlab=xlabel(' ') %set(gcf,'units','points','position',[10,10,110,43]); set(gcf,'units','points','position',[10,10,550,213]); set(gca, 'XTickLabel', []); end %namefig=[folder '\Results\images\' projectname '\' projectname '-freqplot-MP-' num2str(MP(i)) '-' num2str(j)]; namefig=[folder,'\Results\images\' projectname '\' projectname '-freqplot-MP-' num2str(MP(i)) '-' num2str(j)]; % index_dot=strfind(namefig,'.') % namefig(index_dot(end))='-'; %saveas(gcf,namefig,'epsc') set(0, 'CurrentFigure', freqgraph); save2pdf(namefig,freqgraph,20); %exportgraphics(freqgraph,[namefig,'.pdf']) % exportgraphics %appears to be slower. save2pdf is preferred. if dele ==1 savefig(namefig) end end end end \ No newline at end of file