stats_v

 STATS_V   Some basic stats of one or two vector quantities.
    Usage: [mean1, std1, mean2, std2, corr, transfn]=stats(u1,[u2],[printout]);
       or  [mean, std, major, minor]=stats(u,printout);

  Inputs:
         u1 = vector 1 (can be complex)
         u2 = vector 2 (can be complex)
        printout = optional printout control.  If present and positive (e.g. ==1)
                   then stats are printed out to the screen. 

  Outputs: 
        mean1=mean of vector 1
        std1=standard deviation of vector 1
        mean2=mean of vector 2
        std2=standard deviation of vector 2
 
        corr= correlation (complex for complex vectors)
        transfn = transfer function (complex for complex vectors)

        The printout (if printout==1 on input) contains "formatted"
        output, where the angles are printed out relative to True North
        (eg. Eastward=90), and the correlation is expressed as the 
        complex correlation and rotation angle. A positive angle indicates
        that series 1 is rotated counterclockwise from series 2.

0001 function [mean1, std1, mean2, std2, corr, transfn]=stats(u1,u2,printout);
0002 % STATS_V   Some basic stats of one or two vector quantities.
0003 %    Usage: [mean1, std1, mean2, std2, corr, transfn]=stats(u1,[u2],[printout]);
0004 %       or  [mean, std, major, minor]=stats(u,printout);
0005 %
0006 %  Inputs:
0007 %         u1 = vector 1 (can be complex)
0008 %         u2 = vector 2 (can be complex)
0009 %        printout = optional printout control.  If present and positive (e.g. ==1)
0010 %                   then stats are printed out to the screen.
0011 %
0012 %  Outputs:
0013 %        mean1=mean of vector 1
0014 %        std1=standard deviation of vector 1
0015 %        mean2=mean of vector 2
0016 %        std2=standard deviation of vector 2
0017 %
0018 %        corr= correlation (complex for complex vectors)
0019 %        transfn = transfer function (complex for complex vectors)
0020 %
0021 %        The printout (if printout==1 on input) contains "formatted"
0022 %        output, where the angles are printed out relative to True North
0023 %        (eg. Eastward=90), and the correlation is expressed as the
0024 %        complex correlation and rotation angle. A positive angle indicates
0025 %        that series 1 is rotated counterclockwise from series 2.
0026 
0027 if nargin==1,
0028     nvar=1;
0029     printout=0;
0030     u1=denan(u1);
0031 elseif nargin==2,
0032     if length(u2)==1; 
0033       printout=u2; 
0034       nvar=1; 
0035     else 
0036       printout=0;
0037       nvar=2;
0038     end
0039 else 
0040     printout=1; 
0041     nvar=2;
0042 end
0043 if nvar==1,
0044   u1=denan(u1);
0045 else
0046   ind=find(~isnan(u1(:)+u2(:)));
0047   disp('denanning both')
0048   if(isempty(ind)),disp('All NANs'),return,end
0049   u1=u1(ind);
0050   u2=u2(ind);
0051 end
0052 [n,m]=size(u1);
0053 if m>n; u1=conj(u1'); [n,m]=size(u1);
0054 end
0055 
0056 mean1=gmean(u1);
0057 
0058 
0059 if nvar==2
0060 [n2,m2]=size(u2);
0061 if n2~=n; u2=conj(u2');
0062 end
0063 
0064 mean2=gmean(u2);
0065 cv=cov(u1,u2);
0066 std1=sqrt(abs(cv(1,1)));
0067 std2=sqrt(abs(cv(2,2)));
0068 transfn=cv(1,2)/cv(1,1);
0069 corr=cv(1,2)/std1/std2;
0070 else
0071 
0072 cc=cov(real(u1),imag(u1));
0073 std1=sqrt(cc(1,1)+cc(2,2));
0074 theta=0.5*atan2(2.*cc(2,1),(cc(1,1)-cc(2,2)) );
0075 major=sqrt(cc(1,1)*cos(theta)^2+2*cc(1,2)*cos(theta)*sin(theta)...
0076    +cc(2,2)*sin(theta)^2);
0077 minor=sqrt(cc(1,1)*sin(theta)^2-2*cc(1,2)*cos(theta)*sin(theta)...
0078    +cc(2,2)*cos(theta)^2);
0079 theta=theta*180./pi;
0080 
0081 end
0082 
0083 if printout~=0
0084 if nvar==2
0085 if printout>0
0086 disp(' mean1 theta  std1 mean2 theta  std2   corr theta  transfn theta')
0087 end
0088 fprintf('%6.2f%6.1f%6.2f',abs(mean1),90-angle(mean1)*180/pi,std1);
0089 fprintf('%6.2f%6.1f%6.2f',abs(mean2),90-angle(mean2)*180/pi,std2);
0090 fprintf('%8.2f%6.1f',abs(corr),-angle(corr)*180/pi);
0091 fprintf('%7.2f%8.1f',abs(transfn),-angle(transfn)*180/pi);
0092 fprintf('\n')
0093 else
0094 
0095 if printout>0
0096     disp('   mean   theta      std   major    minor   theta')
0097 end
0098 fprintf('%8.2f%8.1f%8.2f',abs(mean1),90-angle(mean1)*180/pi,std1);
0099 fprintf('%8.2f%8.2f%8.1f',major,minor,90-theta);
0100 fprintf('\n')
0101 theta=theta*pi/180;
0102 mean2=major*cos(theta)+i*major*sin(theta);
0103 std2=-minor*sin(theta)+i*minor*cos(theta);
0104 end
0105 
0106 end
0107