disp(' ')
disp(' PSD.m ')
disp(' ver 2.8  October 17, 2012 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the power spectral density ')
disp(' of a time history. ')
disp(' ')
disp(' The time history must be in a two-column matrix format: ')
disp(' Time(sec)  & amplitude ')
disp(' ')
%
close all;
%
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;
%
fig_num=1;
%
[tim,amp,dt,sr,tmx,tmi,n,ncontinue]=enter_time_history();
%
if(ncontinue==1)
%
    [rlab,ylab,iam]=PSD_amplitude_metric();
%
    disp(' ')
    disp(' mean removal? ');
    mr_choice=input(' 1=yes   2=no  ' );
    disp(' ')
%
    disp(' ')
    disp(' Select Window ')
    h_choice = input(' 1=Rectangular 2=Hanning  ');
%
%%%%%%%%%%%%%%%  advise
%
    [df,mmm,NW]=PSD_advise(dt,n);
%
%%%  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);
%
        mag_seg=fix_size(mag_seg);
%
        full=full+mag_seg(1:mH);
    end
%
    den=df*(2*NW-1);
    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(' ');
%
%  Rename matlab matrix with descriptive name
%
    clear power_spectral_density;
%
    freq=fix_size(freq);
    full=fix_size(full);
%
    power_spectral_density=[freq full]; 
%
    [~,fmax]=find_max(power_spectral_density);
%
    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,'%13.6e \t %11.4e \n',freq(k), full(k));    
        end
        fclose(fid1);
    end    
%
    t_string=...
    sprintf(' Power Spectral Density  Overall Level = %6.3g %s ',rms,rlab); 
    x_label=sprintf(' Frequency (Hz)');
    [fig_num]=plot_PSD(fig_num,x_label,ylab,t_string,power_spectral_density);     
%
    disp(' ');
    disp(' Plot the PSD with a specification superimposed? ')
    choice = input(' 1=yes 2=no ');
%
    if(choice == 1)
%
        [THM]=enter_PSD();
%        
        t_string=...
        sprintf(' Power Spectral Density  Overall Level = %6.3g %s ',...
                                                                 rms,rlab); 
        x_label=sprintf(' Frequency (Hz)');
        leg_a=sprintf('Synthesis');
        leg_b=sprintf('Specification');
        [fig_num]=...
        plot_PSD_two(fig_num,x_label,ylab,t_string,...
                                   power_spectral_density,THM,leg_a,leg_b);    
%
    end
%
    out5=...
    sprintf('\n\n The PSD matrix is renamed as "power_spectral_density" ');
    disp(out5);   
%    
end