function[THX]=PSD_octave(THF,fc,imax,oex)
%
%    PSD_octave.m    ver 1.0    July 3, 2007
%
clear sum;
clear ssum;
clear f;
clear ff;
clear amp;
clear n;
clear count;
clear s;
clear grms;
%
f=THF(:,1);
amp=THF(:,2);
%
n = length(amp);
%
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);
%            
                THX(ijk,1)=  ff(ijk);
                THX(ijk,2)=ssum(ijk);
%             
                ijk=ijk+1;
% 
		end
    end