nass-micro-station-measurements/2018-10-12 - Marc/analysis/id31_microstation_raw_12october2018.m

%% Title: id31 microstation in EXP hutch for RAW data
% Date: 12 october 2018

%% NOTE
% With this file you can load the raw mat files
 
%% Description: measure on id31 microstation in exp hutch
 
% FS: =256Hz



% all L28 (31V/m/s) except CH1 L4-C (276V/m/s)
% ch1: marble Z
% ch2: outer frame Ty Z
% ch3: inner frame Tilt Z
% ch4: hexa Z
% ch5: marble H
% ch6: outer frame TY H
% ch7: inner frame Tilt H
% ch8: hexa H
% ch9: hammer

%% measurements 12 october 2018

% excitation Y marble corner
% Measurement1

% excitation Y outer frame corner
% Measurement2

% excitation Y hexa
% Measurement3

% ------------------------------------------------

% excitation Z marble corner
% Measurement4

% excitation Z outer frame corner
% Measurement5

% excitation Z hexa
% Measurement6

% ------------------------------------------------

% excitation X marble corner
% Measurement7

% excitation X outer frame corner
% Measurement8

% excitation X hexa
% Measurement9

%%
microstation=['Marble     '; 'TY         ';'Tilt       '; 'Hexapod    '];

% 
%% PARAMETERS 

beamline='ID31 Nanostation - Hammer testing';
% --------------------------------
%%----------OROS -----------------
ch_max=16;
% --------------------------------

%mult=1e6/276*173;  % --> m/s to micron/s and sensitivity correction for L4-C

nyqhp=2.56; % nyquist 
f_cut=0.5; % cut frequency for high pass filter
t_win=4; % window length in sec
t_ovlp=0; % overlap window in sec

%d=1.3; % distance between vertical sensors


warning off MATLAB:divideByZero

% specify capt # for which to run this
capt=[1:9];

% specify channels for which shut correction must be applied
shunt_ch=[1];

% in case of hammer inpacts specify capt # where it doesnt occur
no_hammer=[];
%no_hammer=[0];
% specify hammer channel (or ch to find peak due to impacts)
shock_ch=9;

%% main loop --------
% ------------------
for i=capt
    
        
    eval(['load Measurement_raw',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
    %c(:,8)=(c(:,5)-c(:,4))/d; % differential Theta Y angle
    
     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,shock_ch,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
    
    
   %% transfer 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