disp(' ')
disp(' PSD.m ')
disp(' ver 2.6  January 11, 2012 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the power spectral density ')
disp(' of a time history that is pre-loaded into Matlab. ')
disp(' ')
disp(' The time history must be in a two-column matrix format: ')
disp(' Time(sec)  & amplitude ')
disp(' ')
%
clear n;
clear Y;
clear amp;
clear tim;
clear N;
clear df;
clear dt;
clear Mag;
clear mag_seg;
clear amp_seg;
clear full;
clear freq;
clear ss;
clear seg;
clear i_seg;
clear rlab;
clear ylab;
%
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);
%
while(1)
    disp(' ')
    disp(' Select amplitude metric ');
    disp('   1=acceleration ');
    disp('   2=velocity ');
    disp('   3=displacement ');
    disp('   4=force ');
    disp('   5=pressure ');
    disp('   6=volts ');
    iam = input('');
    if(iam==1 || iam==2 || iam==3 || iam==4 || iam==5 || iam==6)
        break;
    end
end
%
if(iam==1)
    ylab='Accel (G^2/Hz)';
    rlab='GRMS';
    disp(' Select unit ');
    disp('   1=G ');
    disp('   2=m/sec^2 ');  
    iunit=input(' ');
    if(iunit==2)
      ylab='Accel ((m/sec^2)^2/Hz)';
      rlab='(m/sec^2) RMS';
    end
end
if(iam==2)
    ylab='Velocity ((in/sec)^2/Hz)';
    rlab='(in/sec) RMS';   
    disp(' ')
    disp(' Select unit ');
    disp('   1=in/sec ');
    disp('   2=m/sec '); 
    iunit=input(' ');    
    if(iunit==2)
      ylab='Velocity ((m/sec)^2/Hz)';
      rlab='(m/sec) RMS';     
    end    
end
if(iam==3)
    ylab='Disp (in^2/Hz)';
    rlab='inch RMS';
    disp(' ')
    disp(' Select unit ');
    disp('   1=in ');
    disp('   2=m '); 
    iunit=input(' ');
    if(iunit==2)
      ylab='Disp (m^2/Hz)';
      rlab='m RMS'
    end      
end
if(iam==4)
    ylab='Force (lbf^2/Hz)';
    rlab='lbf RMS';
    disp(' ')
    disp(' Select unit ');
    disp('   1=lbf ');
    disp('   2=N '); 
    iunit=input(' ');
    if(iunit==2)
      ylab='Force (N^2/Hz)';
      rlab='N RMS';
    end         
end
if(iam==5)
    ylab='Pressure (psi^2/Hz)';  
    rlab='psi RMS';
    disp(' ')
    disp(' Select unit ');
    disp('   1=psi ');
    disp('   2=Pa '); 
    iunit=input(' ');
    if(iunit==2)
      ylab='Pressure (Pa^2/Hz)';
      rlab='Pa RMS';
    end         
end
if(iam==6)
    ylab='(Volts^2/Hz)'; 
    rlab='VRMS';
end
%
mu=mean(amp);
sd=std(amp);
mx=max(amp);
mi=min(amp);
rms=sqrt(sd^2+mu^2);
%
disp(' ')
disp(' Input Time History Statistics ')
disp(' ')
tmx=max(tim);
tmi=min(tim);
% 
out3 = sprintf(' start  = %g sec    end = %g sec            ',tmi,tmx);
disp(out3)
%
out0 = sprintf(' number of points = %d ',n);
out1 = sprintf(' mean = %8.4g    std = %8.4g    rms = %8.4g ',mu,sd,rms);
disp(out0)
disp(out1)
%
dt=(tmx-tmi)/(n-1);
sr=1./dt;
df=1./(tmx-tmi);
disp(' ')
disp(' Time Step ');
dtmin=min(diff(tim));
dtmax=max(diff(tim));
%
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)
%
ncontinue=1;
%
if(((dtmax-dtmin)/dt)>0.01)
    disp(' ')
    disp(' Warning:  time step is not constant.  Continue calculation? 1=yes 2=no ')
    ncontinue=input(' ');
end
if(ncontinue==1)
%
disp(' ')
disp(' mean removal? ');
mr_choice=input(' 1=yes   2=no  ' );
disp(' ')
%
disp(' ')
disp(' Select Window ')
h_choice = input(' 1=Rectangular 2=Hanning  ');
%
%%%%%%%%%%%%%%%  advise
%
NC=0;
for(i=1:1000)
%    
    nmp = 2^(i-1);
%   
    if(nmp <= n )
        ss(i) = 2^(i-1);
        seg(i) =n/ss(i);
        i_seg(i) = fix(seg(i));
        NC=NC+1;
    else
        break;
    end
end
%
disp(' ');
out4 = sprintf(' Number of   Samples per   Time per        df               ');
out5 = sprintf(' Segments     Segment      Segment(sec)   (Hz)      dof     ');
%
disp(out4)
disp(out5)
%        
for(i=1:NC)
    j=NC+1-i;
    if j>0
        if( i_seg(j)>0 )
%           str = int2str(i_seg(j));
            tseg=dt*ss(j);
            ddf=1./tseg;
            out4 = sprintf(' %8d \t %8d \t %10.3g  %10.3g \t %d',i_seg(j),ss(j),tseg,ddf,2*i_seg(j));
            disp(out4)
        end
    end
    if(i==12)
        break;
    end
end
%
disp(' ')
NW = input(' Choose the Number of Segments:  ');
disp(' ')
%
mmm = 2^fix(log(n/NW)/log(2));
%
df=1./(mmm*dt);
%
%%%  begin overlap
%
mH=((mmm/2)-1);
%
full=zeros(mH,1);
mag_seg=zeros(mH,1);
%
nov=0;
%
clear amp_seg
%
for ijk=1:(2*NW-1)
%
    amp_seg=zeros(mmm,1);
%
    amp_seg(1:mmm)=amp((1+ nov):(mmm+ nov));  
%
    nov=nov+fix(mmm/2);
%
    [mag_seg]=FFT_core(amp_seg,mmm,mH,mr_choice,h_choice);
%
    sz=size(mag_seg);
    if(sz(2)>sz(1))
        mag_seg=mag_seg';
    end
%
    full=full+mag_seg(1:mH);
end
%
den=df*(2*NW-1);
ms=0.;
fmax=(mH-1)*df;
freq=linspace(0,fmax,mH);
full=full/den;
clear sum;
ms=sum(full);
%
rms=sqrt(ms*df);
%
disp(' ');
out4 = sprintf(' Overall RMS = %10.3g ',rms);
out5 = sprintf(' Three Sigma = %10.3g ',3*rms);
disp(out4)
disp(out5)
disp(' ');
%
fmax=0;
zmax=0;
for(i=1:length(full))
        if(full(i)>zmax)
            zmax=full(i);
            fmax=freq(i);
        end
end
%
out5 = sprintf(' Peak occurs at %10.5g Hz \n',fmax);
disp(out5)
%
disp(' Write the PSD data to a file? ')
choice = input(' 1=yes 2=no ');
%
if(choice == 1)
   disp(' ');   
   disp(' Enter the output filename ');
   filename1 = input(' ','s');
   fid1 = fopen(filename1,'w');  
   for(k=1:mH)
      fprintf(fid1,'%11.4e \t %11.4e \n',freq(k), full(k));    
   end
   fclose(fid1);
end    
%
disp(' ');
disp(' Plot the PSD? ')
choice = input(' 1=yes 2=no ');
%
if(choice == 1)
    figure(1);     
    plot(freq,full)
    xlabel(' Frequency (Hz)');
    ylabel(ylab); 
    at = sprintf(' Power Spectral Density  Overall Level = %6.3g %s ',rms,rlab);   
    title(at);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log'); 
    grid;
%
    psd_max=max(full);
    psd_min=min(full); 
%    
    if(psd_max < 10000)
        ymax=10000;
    end
    
    for(i=-30:30)
        if(psd_max < 10^i)
            ymax=10^i;
            break;
        end
    end
%
    ymin=ymax/100;    
    if(psd_min < psd_max/1000)
       ymin=ymax/1000;    
    end 
    if(psd_min < psd_max/1000)
       ymin=ymax/1000;    
    end    
    fmax=max(freq);
    set(gca, 'XTickMode', 'manual');
%
    pfmax=10000;
    set(gca,'xtick',[10 100 1000 10000]); 
%
    pfmax=fmax;
    if(fmax<10000)
        pfmax=10000;
        set(gca,'xtick',[10 100 1000 10000]); 
    end
    if(fmax<5000)
        pfmax=5000;
        set(gca,'xtick',[10 100 1000 5000]);   
    end
    if(fmax<2000)
        pfmax=2000;
        set(gca,'xtick',[10 100 1000 2000]);      
    end
    if(fmax<1000)
        pfmax=1000;
        set(gca,'xtick',[1 10 100 1000]);      
    end
    if(fmax<100)
        pfmax=100;
        set(gca,'xtick',[0.1 1 10 100]);      
    end    
    if(fmax<10)
        pfmax=10;
        set(gca,'xtick',[0.01 0.1 1 10]);      
    end     
    if(fmax<1)
        pfmax=1;
        set(gca,'xtick',[0.001 0.01 0.1 1]);      
    end       
    fmin=df;
    axis([fmin,pfmax,ymin,ymax]);  
%
    disp(' ');
    disp(' Plot the PSD with a specification superimposed? ')
    choice = input(' 1=yes 2=no ');
%
    if(choice == 1)
%
        disp('');
        disp(' Select file input method ');
        disp('   1=external ASCII file ');
        disp('   2=file preloaded into Matlab ');
        file_choice = input('');
%
        if(file_choice==1)
%%
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);  
%%
%    disp(' Enter the input filename ');
%    filename = input(' ','s');
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g',[2 inf]);
        THM=THM';
    else
        THM = input('  Enter the matrix name:  ');
    end
%
    THM(:,2);
    THM(:,1);
    disp(' ')
 %           
    figure(2);   
    plot(freq,full,THM(:,1),THM(:,2),'r');
    legend ('Synthesis','Specification',2);
    axis([fmin,pfmax,ymin,ymax]);    
    xlabel(' Frequency (Hz)');
    ylabel(ylab); 
    at = sprintf(' Power Spectral Density  Overall Level = %6.3g %s ',rms,rlab);   
    title(at);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log'); 
    grid;
end
end
%
%  Rename matlab matrix with descriptive name
%
    clear power_spectral_density;
    sz=size(freq);
    if(sz(2)>sz(1))
        freq=freq';
    end  
    sz=size(full);
    if(sz(2)>sz(1))
        full=full';
    end
    power_spectral_density=[freq full]; 
    out5 = sprintf('\n\n The PSD matrix is renamed as "power_spectral_density" ');
    disp(out5);   
%    
end