psd_integ.m

功率谱估计

disp(' ')
disp(' psd_integ.m   ver 1.1   May 5, 2005')
disp(' by Tom Irvine   Email: [e-mail]')
disp(' ')
disp(' This program calculates the overall level of ')
disp(' a power spectral density function. ')
disp(' ')
disp(' Select data input method ')
m = input(' 1=input file  2=keyboard entry ');
disp(' ')
clear f;
clear a;
clear s;
%
if(m==1)
    disp(' The input file must be pre-loaded into Matlab.')
    disp(' ')
    disp(' The input file must have two columns: frequency(Hz) & accel(G^2/Hz).')
    disp(' ')
    THF = input(' Enter the input PSD filename:  ');
    f=THF(:,1);
    a=THF(:,2);
    n=length(f);
else
    n= input(' Enter the number of input points: ');
    disp(' ')
    disp(' The frequency unit is (Hz).  The amplitude unit is (G^2/Hz) ')
    disp(' ')
    for(i=1:n)
        out1=sprintf(' Enter frequency %d:  ',i);
        f(i)=input(out1);
        out2=sprintf(' Enter amplitude %d:  ',i);
        a(i)=input(out2);  
        disp(' ')
    end
end
disp(' ')
%
[s,grms]=calculate_PSD_slopes(f,a);
%
out5 = sprintf('\n Overall Acceleration = %10.3g GRMS',grms);
disp(out5);
%
vrms = APSD_VPSD(f,a);
out6 = sprintf('\n Overall Velocity     = %10.3g inch/sec RMS',vrms);
disp(out6);
%    
drms = APSD_DPSD(f,a);
out7 = sprintf('\n Overall Displacement = %10.3g inch RMS',drms);
disp(out7);
out7 = sprintf('                      = %10.3g inch 3-sigma\n',3*drms);
disp(out7);
%  
disp(' ')
disp(' Plot input PSD ?');
choice = input(' 1=yes  2= no  ');
ymin=min(a);
ymax=max(a);
pymax=ymax;
pymin=ymax/100.;
if(ymax<10)
    pymax=10;
end
if(ymax<1)
    pymax=1;
end
if(ymax<0.1)
    pymax=0.1;
end
if(ymax<0.01)
    pymax=0.01;
end
if(ymax<0.001)
    pymax=0.001;
end
if(ymax<0.0001)
    pymax=0.0001;
end
if(ymin <= 1 )
    pymin = 0.1;
end
if(ymin <= 0.1 )
    pymin = 0.01;
end
if(ymin <= 0.01 )
    pymin = 0.001;
end
if(ymin <= 0.001 )
    pymin = 0.0001;
end
if(ymin <= 0.0001 )
    pymin = 0.00001;
end
%
if(choice==1)
   plot(f,a);
   ylabel('Accel (G^2/Hz)');
   xlabel('Frequency (Hz)');
   out=sprintf(' Power Spectral Density %7.3g GRMS ',grms);
   title(out);
   grid;
   set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
   axis([f(1),f(n),pymin,pymax]);
   if(f(n)==2000 & f(1)==10)
        set(gca, 'XTickMode', 'manual');
        set(gca,'xtick',[10 100 1000 2000]);
   end
   if(f(n)==2000 & f(1)==20)
        set(gca, 'XTickMode', 'manual');
        set(gca,'xtick',[20 100 1000 2000]);
   end
end