disp(' ');
disp('   PSDoct.m    version 1.2,   July 3, 2007 ');
disp('   by Tom Irvine ');
disp('   Email:  tomirvine@aol.com ');
disp(' ');
disp('  This program converts a PSD with constant bandwidth');
disp('  to an octave format.');
disp(' ');
disp('  The input file must be frequency(Hz) and PSD(G^2/Hz) ');
disp('  The format is free. ');
%
clear fc;
clear imax;
clear ioct;
clear sum;
clear ssum;
clear f;
clear ff;
clear amp;
clear n;
clear count;
clear s;
clear grms;
%
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);
%      
    fid = fopen(filename,'r');
    THM = fscanf(fid,'%g %g',[2 inf]);
    THM=THM';
else
    THM = input(' Enter the matrix name:  ');
end
%
f=THM(:,1);
amp=THM(:,2);
%
n = length(amp);
%
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);    
%
disp(' ');
disp(' set limits... ');
%
	for(i=1:imax)
			sum(i)=0.;
			fl(i)=fc(i)/(2.^oex);
	end
	for(i=1:(imax-1))
	
			fu(i)=fl(i+1);	
	end
	fu(imax)=fc(i)*(2.^oex);
%
    disp(' ');
	disp('  counts...');
%
	for(i=1:imax)	
		sum(i)=0.;
		count(i)=0;
	end
%
    for(k=1:n)
%        
		for(i=1:imax)
%		
			if( f(k)>= fl(i) & f(k) < fu(i))
%
				sum(i)=sum(i)+ amp(k);
				count(i)=count(i)+1;
			end
		end
	end
%
   disp(' ');
   disp('  calculate output data...');
%
    iflag=0;
    ijk=1;
	for(i=1:imax)
%	
		if( fl(i) > f(n))
		   break;
		end
         
		if(iflag==0)
			if(count(i)>=1 && count(i+1)>=2)
				iflag=1;
            end
        end

		if( iflag==1 & sum(i) > 1.0e-20) 
%		
			ssum(ijk)=sum(i)/count(i);
            ff(ijk)=fc(i);
            ijk=ijk+1;
% 
		end
    end
%
if( ijk>=2)
%    
    [s,grms]=calculate_PSD_slopes(ff,ssum);
%
    disp(' ');
    disp(' Plot the PSD? ')
    choice = input(' 1=yes 2=no ');
%
    if(choice == 1)
        figure(1);     
        plot(ff,ssum)
        grid;
        set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
        xlabel(' Frequency (Hz)');
        ylabel(' Accel (G^2/Hz)');
        at = sprintf(' Power Spectral Density  Overall Level = %6.3g GRMS ',grms);   
        title(at);      
    end
else
    disp(' ')
    disp(' Error ')
end
%
%  Rename matlab matrix with descriptive name
%
    clear psd_oct;
    psd_oct=[ff;ssum]'; 
    out5 = sprintf('\n\n The PSD matrix is renamed as "psd_oct" ');
    disp(out5);   
%