disp(' ')
disp(' signal_stats.m ')
disp(' ver 3.5   August 28, 2011 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program gives statistics and plots for ')  
disp(' PSD, SRS and time history functions.')
disp(' ')
%
clear THM;
clear amp;
clear pa;
clear tim;
clear table;
%
clear a2;
clear b2;
clear f2;
clear n2;
clear sz2;
clear THF2;
clear size;
clear length;
clear max;
clear min;
clear diff;
clear disp;
clear abs;
clear temp;
clear y;
clear mat;
%
mina= 1.0e+90;
maxa=-1.0e+90;
%
iflag=0;
%
disp(' Select function: '); 
disp('      1=PSD  ');
disp('      2=SRS  ');
disp('      3=VRS  ');
disp('      4=SPL  ');
disp('      5=time history');
m=input(' ');
disp(' ');
%
if(m<=3)
%
    disp(' Select data input method ')
    disp('  1=data preloaded in matlab ')  
    disp('  2=keyboard entry ');
    disp('  3=external ASCII file ');
    disp('  4=Excel file ');
%
    method=input('');                 
    disp(' ')
%
    clear a;
    clear b;
%
    if(method==1)
        THF = input(' Enter the input filename:  ');
        sz=size(THF);
        f=double(THF(:,1));        
        a=double(THF(:,2));
        n=length(f);       
%
        if(m==2)
            if(sz(2)==2)
                table=[f ; a]';               
            else
                b=double(THF(:,2)); 
                table=[f ; a ; b]';    
            end
        else
            table=[f ; a]';
        end
%        
    end   
    if(method==2)
        n=input(' Input the number of coordinates  ');	
%
        for(i=1:n)
            if(m==1)
                out = sprintf('\n Enter Frequency %ld ',i);
            else
                out = sprintf('\n Enter Natural Frequency %ld ',i);      
            end
            disp(out);
            f(i)=input('');
%
            out = sprintf('\n Enter amplitude %ld ',i);
            disp(out);
            a(i)=input('');  
        end
        n=length(f);     
        table=[f ; a]';  
        THF=table;
    end
    if(method==3)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
%      
        fid = fopen(filename,'r');
%        THF = fscanf(fid,'%g %g',[2 inf]);
%        THF=THF';
         THF=load(filename);
        sz=size(THF);
        f=double(THF(:,1));
        a=double(THF(:,2));
        if(sz(2)==3)
            b=double(THF(:,3));
            table=[f , a , b];   
        else
            table=[f , a];
        end
        n=length(f);            
    end 
%
    if(method==4)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THF = xlsread(xfile);
%        
        f=double(THF(:,1));
        a=double(THF(:,2));  
%       
        n=length(f);
        table=[f ; a]; 
    end   
%
end
if(m==1)  % PSD
    [n,mina,maxa]=coordinates_ss(m,THF,mina,maxa,iflag);
    io=0;
    while(io < 3 )
        disp(' ');
        disp(' Additional Options         ');
        disp(' 1=add dB margin            ');
        disp(' 2=convert to octave format ');
        disp(' 3=scale to GRMS            ');
        disp(' 4=exit                     ');   
        io=input(' ');
        disp(' ')
        if(io==1)
            dB=input(' Input Margin (dB) ');
            scale=10^(dB/10);
            THF(:,2)=THF(:,2)*scale;
            [n,mina,maxa]=coordinates_ss(m,THF,mina,maxa,iflag);       
        end   
        if(io==2)
            disp(' ');
            disp(' Enter output format ')
            disp(' 1= octave  2= 1/3 octave  3= 1/6 octave  4= 1/12 octave ');
            ioct=input(' ');
%
            disp(' ');
            disp(' set octave... ');
% 
	        if(ioct==1)     % octave
		        oex=1./2.;
	        end
	        if(ioct==2)     % one-third octave
		        oex=1./6.;
	        end
	        if(ioct==3)     % one-sixth octave
		        oex=1./12.;
	        end
	        if(ioct==4)     % one-twelfth octave
		        oex=1./24.;
    	    end
%                
            [fc,imax]=octaves(ioct);            
            clear THX;
            [THX]=PSD_octave(THF,fc,imax,oex);
            clear THF;
            THF=THX;
            [n,mina,maxa]=coordinates_ss(m,THF,mina,maxa,iflag);    
        end
        if(io==3)
            disp(' ');
            GRMS=input(' Input GRMS ');
            [n,mina,maxa,THF]=coordinates_ss_scale(m,THF,mina,maxa,iflag,GRMS);       
        end
%
        if(io~=4)
            psd_new=THF;
            disp(' ');
            disp(' Matlab output array = psd_new ');
        end
%
    end
end
%
if(m==2 | m==3)  % SRS or VRS
    [n,mina,maxa]=coordinates_ss(m,THF,mina,maxa,iflag);
    grid on;
    io=0;
    while(io < 2 )
        disp(' ');    
        disp(' Additional Options ');
        disp(' 1=add dB margin    ');
        if(m==2)
            disp(' 2=Superimpose another SRS curve');
        else
            disp(' 2=Superimpose another VRS curve');           
        end
        disp(' 3=exit             ');   
        io=input(' ');
        disp(' ')
        if(io==1)
            dB=input(' Input Margin (dB) ');
            scale=10^(dB/20);
            THF(:,2)=THF(:,2)*scale;
            if(sz(2)==3)
                THF(:,3)=THF(:,3)*scale;      
            end
            [n,mina,maxa]=coordinates_ss(m,THF,mina,maxa,iflag);       
        end   
        if(io==2)
            hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            disp(' Select data input method ')
            disp('  1=data preloaded in matlab ')  
            disp('  2=keyboard entry ');
            disp('  3=external ASCII file ');
            disp('  4=Excel file ');
%
            method2=input('');                 
            disp(' ')
%
            if(method2==1)
                THF2 = input(' Enter the input filename:  ');
                sz2=size(THF2);
                f2=double(THF2(:,1));        
                a2=double(THF2(:,2));
                n2=length(f2);       
%
                if(m==2)
                    if(sz2(2)==2)
                        table=[f2 ; a2]';               
                    else
                        b2=double(THF2(:,2)); 
                        table=[f2 ; a2 ; b2]';    
                    end
                else
                    table=[f2 ; a2]';
                end
%        
            end   
%
            if(method2==2)
                n2=input(' Input the number of coordinates  ');	
%
                for(i=1:n2)
                    if(m==1)
                        out = sprintf('\n Enter Frequency %ld ',i);
                    else
                        out = sprintf('\n Enter Natural Frequency %ld ',i);      
                    end
                    disp(out);
                    f2(i)=input('');
%
                    out = sprintf('\n Enter amplitude %ld ',i);
                    disp(out);
                    a2(i)=input('');  
                end
                n2=length(f2);     
                table=[f2 ; a2]';
                THF2=table; 
            end
            if(method2==3)
                [filename, pathname] = uigetfile('*.*');
                filename = fullfile(pathname, filename);
%      
                fid = fopen(filename,'r');
%        THF = fscanf(fid,'%g %g',[2 inf]);
%        THF=THF';
                THF2=load(filename);
                sz2=size(THF2);
                f2=double(THF2(:,1));
                a2=double(THF2(:,2));
                if(sz2(2)==3)
                    b2=double(THF2(:,3));
                    table=[f2 , a2 , b2];   
                else
                    table=[f2 , a2];
                end
                n2=length(f2);            
            end 
%
            if(method2==4)
                [filename, pathname] = uigetfile('*.*');
                xfile = fullfile(pathname, filename);
%        
                THF2 = xlsread(xfile);
%        
                f2=double(THF2(:,1));
                a2=double(THF2(:,2));  
%       
                n2=length(f2);
                table=[f2 ; a2]; 
            end               
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
            iflag=1;
            [n2,mina,maxa]=coordinates_ss(m,THF2,mina,maxa,iflag);  
        end
    end
    hold off;
end
if(m==4)  % SPL
     [n,mina,maxa]=coordinates_ss(m,THF,mina,maxa,iflag); 
     io=0;
     while(io < 2 )
        disp(' ');    
        disp(' Additional Options ');
        disp(' 1=add dB margin    ');
        disp(' 2=exit             ');   
        io=input(' ')
        disp(' ')
        if(io==1)
            dB=input(' Input Margin (dB) ');
            THF(:,2)=THF(:,2)+dB;
            [n,mina,maxa]=coordinates_ss(m,THF,mina,maxa,iflag);            
        end
    end   
end    
%
if( m==5) % Time History    
    disp(' The time history must be in a two-column matrix format: ')
    disp(' Time(sec)  & amplitude ')
    disp(' ')    
%
    disp(' Select file input method ');
    disp('   1=external ASCII file ');
    disp('   2=file preloaded into Matlab ');
    disp('   3=Excel file ');
    file_choice = input('');
%
    if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g',[2 inf]);
        THM=THM';
    end
    if(file_choice==2)
        THM = input(' Enter the matrix name:  ');
    end
    if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
    end
%
    amp=double(THM(:,2));
    tim=double(THM(:,1));
    n = length(amp);
%disp(' mean values ')
    mu=mean(amp);
    sd=std(amp);
    mx=max(amp);
    mi=min(amp);
    rms=sqrt(sd^2+mu^2);
    kt=0.;
    sk=0.;
%1
    for i=1:n
        if( amp(i)==mx)
            tt_max=tim(i);
        end
        if( amp(i)==mi)
            tt_min=tim(i);
        end
        kt=kt+(amp(i)-mu)^4;
        sk=sk+(amp(i)-mu)^3;
    end      
    kt=kt/(n*sd^4);
    sk=sk/(n*sd^3);
%1
    disp(' ')
    disp(' Time Stats ')
    tmx=max(tim);
    tmi=min(tim);
    duration=(tmx-tmi);
% 
    out3 = sprintf('   start  = %g sec    end = %g sec   duration = %g sec  ',tmi,tmx,duration);
    disp(out3)
%
    dt=duration/(n-1);
    out5 = sprintf('\n number of samples = %d  ',n);
    disp(out5)
%
    clear t;
    t=tim;
    disp(' ')
    disp(' Time Step ');
    clear difft;
    difft=diff(t);
    dtmin=min(difft);
    dtmax=max(difft);
%
    out4 = sprintf(' dtmin  = %8.4g sec  ',dtmin);
    out5 = sprintf(' dt     = %8.4g sec  ',dt);
    out6 = sprintf(' dtmax  = %8.4g sec  ',dtmax);
    disp(out4)
    disp(out5)
    disp(out6)
%
    disp(' ')
    disp(' Sample Rate ');
    out4 = sprintf(' srmin  = %8.4g samples/sec  ',1/dtmax);
    out5 = sprintf(' sr     = %8.4g samples/sec  ',1/dt);
    out6 = sprintf(' srmax  = %8.4g samples/sec  \n',1/dtmin);
    disp(out4)
    disp(out5)
    disp(out6)
    clear t;
%
    disp(' ')
    disp(' Amplitude Stats ')
     out0 = sprintf(' number of points = %d ',n);
    out1 = sprintf(' mean = %8.4g    std = %8.4g    rms = %8.4g ',mu,sd,rms);
    out2a = sprintf(' max  = %9.4g  at  = %8.4g sec             ',mx,tt_max);
    out2b = sprintf(' min  = %9.4g  at  = %8.4g sec            \n',mi,tt_min);
    disp(out0);
    disp(out1);
    disp(out2a);
    disp(out2b);
%
    amax=abs(mx);
    amin=abs(mi);
    if(amax<amin)
        amax=amin;
    end
    crest=amax/rms;
    out4 = sprintf(' crest factor = %8.4g ',crest);
    out5 = sprintf(' skewness = %8.4g  kurtosis  = %8.4g ',sk,kt);
    disp(out4);
    disp(out5);
%
    nzc=0;
    pzc=0;
    for(i=1:n-1)
%
       if(amp(i) <=0 & amp(i+1) > 0)
           pzc=pzc+1;
       end
%       
       if(amp(i) >=0 & amp(i+1) < 0)
           nzc=nzc+1;      
       end 
    end
%
    out6=sprintf('\n positive zero crossings = %d \n negative zero crossings = %d \n',pzc,nzc);
    disp(out6);    
    %
    j=1;
%
    clear np;
    clear pp;
    clear abs;
    clear pa;
%%%    pa=zeros(n,1);
    np=0;
    pp=0;
%
    for(i=2:(n-1))
        if( amp(i)>=amp(i-1) & amp(i)>=amp(i+1) & amp(i)>0 )
            pa(j)=abs(amp(i));
            j=j+1;
            pp=pp+1;
        end
        if( amp(i)<=amp(i-1) & amp(i)<=amp(i+1)& amp(i)<0 )
            pa(j)=abs(amp(i));
            j=j+1;
            np=np+1;
        end
    end
%
%
    out9 =sprintf(' number of positive peaks = %d \n number of negative peaks = %d ',pp,np);
    disp(out9);
    out10 =sprintf(' total peaks = %d ',pp+np);
    disp(out10);  
    disp(' ')
    fignum=0;
    disp(' Plot time history file? ')
    choice = input(' 1=yes 2=no ');
    if(choice == 1)
        fignum=fignum+1;
        figure(fignum);
        plot(tim,amp)
        xlabel(' Time (sec)');
        zoom on;
    end
%
    disp(' ')
    disp(' Plot histogram? ')
    choice = input(' 1=yes 2=no ');
%
    if(choice == 1)
        while(1)
            disp(' ');
            nbar=input(' Enter the number of bars ');
            fignum=fignum+1;
            figure(fignum);    
            hist(amp,nbar)
            title(' Histogram ');
            ylabel(' Counts ');
            xlabel('Amplitude');
            disp(' ');
            disp(' Re-plot with different number of bars? ');
            disp(' 1=yes  2=no ');
            irp=input(' ');
            if(irp==2)
                break;
            end
        end
%
        disp(' ')
        disp(' Plot cumulative histogram? ')
            cum_choice = input(' 1=yes 2=no '); 
% 
            if(cum_choice==1)
                [nnn,xout]=hist(amp,nbar);
                c_elements = cumsum(nnn);
                fignum=fignum+1;
                figure(fignum);
                bar(xout,c_elements);
                title(' Cumulative Histogram ');
                ylabel(' Counts ');
                xlabel('Amplitude');
                grid on;
            end
%
        end
%
    disp(' ')
    disp(' Plot absolute value peak distribution? ')
    choice = input(' 1=yes 2=no ');
%
    if(choice==1)
            while(1)
                disp(' ');
                nbar=input(' Enter the number of bars ');
                fignum=fignum+1;
                figure(fignum);  
                hist(pa,nbar)
                title(' Histogram of Absolute Values of Peaks');
                ylabel(' Counts ');
                disp(' ');
                disp(' Re-plot with different number of bars? ');
                disp(' 1=yes  2=no ');
                irp=input(' ');
                if(irp==2)
                    break;
                end
            end
    end
end
%
if(iflag==1)
    hold off;
end