nass-micro-station-measurem.../2018-10-15 - Marc/analysis-marc/id31_microstation_15october2018.m

267 lines
8.2 KiB
Matlab

% Title: id31 microstation in EXP hutch
% Date: 15 october 2018
% Description: measure on id31 microstation in exp hutch
% FS: =256Hz
%% 15 october 2018 --------------
% L4-c sensor at 276V/m/s
% ch1: Tilt frame Z upstream
% ch2: Tilt frame Z downstream
% ch3: Ty frame Y
% TY motor off --> on at ~300sec
% capt1
% Tilt OFF --> ON at ~ 326sec
% capt2
% ----------------------------
% ch1: Hexa Z
% ch2: Tilt frame Z downstream
% ch3: Ty frame Y
%
% Hexa ON --> OFF at ~ 406sec (tilt ON)
% capt3
%
% Hexa OFF - Slip ring ON at ~ 300sec then spindle ON at ~ 620sec (tilt ON)
% capt4
%% Marble measurements ----
% ch1 floor Z
% ch2 marble Z
% ch3 floor Y
% ch4 marble Y
% capt5
%% PARAMETERS
beamline='ID31 Nanostation ';
% --------------------------------
%%----------OROS -----------------
ch_max=16;
% --------------------------------
mult=1e6/276*173; % --> m/s to micron/s and sensitivity correction
nyqhp=2.56; % nyquist
f_cut=0.5; % cut frequency for high pass filter
t_win=4; % window length in sec
t_ovlp=3; % overlap window in sec
d=1; % distance between vertical sensors.
warning off MATLAB:divideByZero
% specify capt # for which to run this
capt=1:5;
% specify channels for which shut correction must be applied
% shunt_ch_a=1:3;
% shunt_ch_b=1:4;
% in case of hammer inpacts specify capt # where it doesnt occur
no_hammer=1:5;
%no_hammer=0;
% specify hammer channel (or ch to find peak due to impacts)
shock_ch=1;
%% main loop --------
% ------------------
for i=capt
eval(['load Measurement',num2str(i)])
freq_max=Track1_TrueBandWidth;
dts=1/(freq_max*nyqhp);
freq=linspace(0,freq_max,t_win*freq_max);
wo=2*pi*freq;
for k=1:ch_max
vname=['Track',num2str(k)];
array_exist(k)=ismember(vname,who);
end
non_zero=find(array_exist);
for z=non_zero(1):length(non_zero)
track_nb=['Track',num2str(z)]';
eval(['data(:,z)=Track',num2str(z),';']);
end
c=data*mult;
%-------------
nbch=size(c,2);
%-------------
r=length(c);
if r/2~=fix(r/2) % loop to test for odd or even nb of samples
c=c(1:r-1,:); % take only even
else
end
%------------------------------
time=linspace(0,length(c)*dts,length(c));
for j=nbch %shunt_ch
[c(:,j),c_shut]=shut_c(c(:,j),1/dts); % correct for shunt
end
% compute differential level when necessary and store it as 4th column
if i<3
c(:,4)=(c(:,2)-c(:,1))/d; % divide by d to obtain Theta Y angle
end
b=find(no_hammer==i); % if i==1 | i==2 | i==6
if b~=0
[psd_v,integ_v,psd_d,integ_d]=integrated_psd(c,t_win,t_ovlp,nyqhp,dts);
[frz_cut,crsp,pwsp,coherz,nsp]=fqresp(c,1,t_win,t_ovlp,nyqhp,dts);
[frh_cut,crsp,pwsp,coherz,nsp]=fqresp(c,3,t_win,t_ovlp,nyqhp,dts);
else
thresh=0.5; % threshold of max value
sep=2.5; % separation minimum of peaks in sec
pre_ev=2; % pre event delay in sec
pos_ev=2; % post event delay in sec
[ti,t_impact]=findpeaks(c(:,shock_ch),'minpeakheight',max(c(:,shock_ch))*thresh,'minpeakdistance',ceil(sep/dts));
% find times at which there are impacts (threshold of max and separated by sep sec)
psd_v=zeros((pre_ev+pos_ev)/dts/nyqhp,nbch);
psd_d=zeros((pre_ev+pos_ev)/dts/nyqhp,nbch);
frz_cut=zeros((pre_ev+pos_ev)/dts/nyqhp,nbch);
for k=1:length(t_impact)
ibeg=fix(t_impact(k)-(pre_ev/dts));
iend=fix(t_impact(k)+(pos_ev/dts));
freq_s=linspace(0,freq_max,t_win/2*freq_max);
if ibeg>1 && iend<length(c) % eliminate indexes outside data range
[psd,integ_v,psd_int,integ_d]=integrated_psd(c(ibeg:iend,:),t_win,t_ovlp,nyqhp,dts);
psd_v=psd+psd_v;
psd_d=psd_int+psd_d;
[frz,crsp,pwsp,coherz,nsp]=fqresp(c(ibeg:iend,:),shock_ch,t_win,t_ovlp,nyqhp,dts);
frz_cut=frz+frz_cut;
end
end
psd_v=psd_v/length(t_impact);
psd_d=psd_d/length(t_impact);
frz_cut=frz_cut/length(t_impact);
end
drms=max(integ_d); % compute rms level
dc=hpfint(c,f_cut,dts); % filter and integrate in time domain
dppc=hpdpp(dc,t_win,t_ovlp,1,dts); % compute peak to peak level
% tranfer function, cross spectrum, power spectr. and coherence w.r.t. ch1
%----------------------------------------------------------------------------------------------------
eval(['c',num2str(i),'=c;'])
eval(['dc',num2str(i),'=dc;'])
eval(['dppc',num2str(i),'=dppc;'])
eval(['drms',num2str(i),'=drms;'])
eval(['psd_v',num2str(i),'=psd_v;']) % already integrated in OROS
eval(['psd_d',num2str(i),'=psd_d;'])
eval(['integ_v',num2str(i),'=integ_v;'])
eval(['integ_d',num2str(i),'=integ_d;'])
eval(['frz',num2str(i),'=frz_cut;'])
eval(['frh',num2str(i),'=frh_cut;'])
% eval(['frx',num2str(i),'=frx;'])
% eval(['coherz',num2str(i),'=coherz;'])
eval(['time',num2str(i),'=time;'])
clear data c dc psd psd_v psd_d time c_shut % clean up the mess
end
%% Plot settings for colors and linewidth----
proname(1)={'LineWidth'};
proname(2)={'Color'};
proname(3)={'LineStyle'};
val(1,1) = {.5} ;val(1,2) = {[0.6 0.2 1]} ;val(1,3) = {'-'};
val(2,1) = {2} ;val(2,2) = {[0 0 1]} ;val(2,3) = {'-'};
val(3,1) = {2} ;val(3,2) = {[0.25 0.9 0.65]} ;val(3,3) = {'-'};
val(4,1) = {2} ;val(4,2) = {[0 1 0]} ;val(4,3) = {'-'};
val(5,1) = {0.5} ;val(5,2) = {[1 0.4 0.4]} ;val(5,3) = {'-'};
val(6,1) = {2} ;val(6,2) = {[1 0 0]} ;val(6,3) = {'-'};
val(7,1) = {1} ;val(7,2) = {[0.8 0.8 0.8]} ;val(7,3) = {'-'};
val(8,1) = {2} ;val(8,2) = {[0.1 0.1 0.2]} ;val(8,3) = {'-'};
val(9,1) = {1} ;val(9,2) = {[0.7 0.8 0.4]} ;val(9,3) = {'-'};
val(10,1) = {2} ;val(10,2) = {[0.7 0.8 0.2]} ;val(10,3) = {'-'};
val(11,1) = {1} ;val(11,2) = {[0.9 0.7 0.35]} ;val(11,3) = {'-'};
val(12,1) = {2} ;val(12,2) = {[1 0.8 0.3]} ;val(12,3) = {'-'};
val(13,1) = {1} ;val(13,2) = {[0.5 0.4 0.3]} ;val(13,3) = {'-'};
val(14,1) = {2} ;val(14,2) = {[0.5 0.3 0.2]} ;val(14,3) = {'-'};
%% PLOT legends, titles,...
xlab1='Frequency in Hz';
xlab2='Time in sec';
ylab1='Amplification';
ylab2='PSD in $\frac{\mu{}m^{2}}{Hz}$';
ylab3='PSD in $\frac{\mu{}m s)^{2}/Hz';
ylab4='Displacement in ${\mu{}m}$';
ylab5='Displacement in \mum';
ylab7='Coherence';
font_s=14;
% ---------------------------------
% tit_1=[beamline,' - Amplification wrt Floor (Z)'];
% tit_2=[beamline,' - Amplification wrt Floor (Y)'];
% tit_3=[beamline,' - Amplification wrt Floor (X)'];
tit_4=[beamline,' - Vertical (Z) PSD'];
tit_7=[beamline,' - Horizontal (X) PSD'];
tit_6=[beamline,' - Horizontal (Y) PSD'];
legend1=['Floor','Marble','Location','NorthEast'];
% legend2=['''Floor OFF'',''Frame EM OFF'',''Floor ON'',''Frame EM ON'',''Location'',''NorthEast'''];
% legend3=['''EM ON'',''EM OFF'',''Location'',''NorthWest'''];
%% Response of Marble - Y
h1 = newFigure(16,12);
h=semilogy(freq,abs([psd_d5(:,[3 4])]));
set(h,proname,val([1 6],1:3))
eval(['leg1 = legend(',legend1,'); set(leg1, ''Interpreter'', ''latex'')'])
titlabel_font(tit_6,xlab1,ylab2,font_s);
axis([0 100 1e-11 1e-1])
grid
saveas(gcf,'psd_marble_y','fig')
print -dpng psd_marble_y
exportFigure(h1,'psd_marble_y', 'pdf')
%% Response of Marble - Z
h1 = newFigure(16,12);
h=semilogy(freq,abs([psd_d5(:,[1 2])]));
set(h,proname,val([1 6],1:3))
eval(['leg1 = legend(',legend1,'); set(leg1, ''Interpreter'', ''latex'')'])
titlabel_font(tit_4,xlab1,ylab2,font_s);
axis([0 100 1e-11 1e-1])
grid
saveas(gcf,'psd_marble_z','fig')
print -dpng psd_marble_z
exportFigure(h1,'psd_marble_z', 'pdf')
%% spectrograms
spt = strogram_h(c1(:,3),4,3,1/256,2.56,'ID31 nanostation Ty Y - Ty ON @ 300s',1,30,4,12,'egend1','legend2');
exportFigure(h1,'spectrogram_Ty_y', 'pdf')
spt = strogram_h(c2(:,1),4,3,1/256,2.56,'ID31 nanostation Tilt Z - Tilt ON @ 320s',1,30,4,12,'egend1','legend2');
exportFigure(h1,'spectrogram_Tilt_z', 'pdf')
spt = strogram_h(c3(:,1),4,3,1/256,2.56,'ID31 nanostation Hexa Z - Hexa OFF @ 410s',1,30,4,12,'egend1','legend2');
exportFigure(h1,'spectrogram_hexa_y', 'pdf')
spt = strogram_h(c4(:,1),4,3,1/256,2.56,'ID31 nanostation Hexa Z - SlipRing ON @ 300s; Spindle ON @ 620s',1,30,4,12,'egend1','legend2');
exportFigure(h1,'spectrogram_slip_spindle_y', 'pdf')