function P=hspec_ss(x,y,Fs)
  
% 7/31/2002 - modify so that an even number of samples is used,
% since odd numbers end up looking for non-integer operands  
%
% with smoothing
%
% HSPEC  Power, co-, and quad-spectrum estimate of one or two data
% 	sequences with user specified confidence interval (eg. 90, 95, ...
%	percent).  The power and co-spectrum is defined such that its
% 	integral is the variance and co-variance (without the mean),
%	respectively.  The confidence intervals are determined using
% 	an estimated chi-squared distribution.  The supplementary FUNCTION
% 	zalph.m MUST accompany this m-file for confidence limits to operate. 
% 	P=HSPEC(X,Y,Fs) performs FFT analysis of the two signals
%	X and Y where Fs is the sampling frequency.  If Y=0 then
%	the FFT analysis is done only on signal X.
%	Smoothing is accomplished using a linearly increasing
%	interval averaging technique in the frequency domain.
%	HSPEC returns the array
%		P=[ f psdx psdy co qu]
%	where
%	  f    = frequency
%	  psdx = X-signal power spectral density, units of (x)^2 per unit of Fs
%	  psdy = Y-signal power spectral density, units of (y)^2 / unit of Fs
%	  co   = Co-spectrum, units of xy / unit of Fs
%	  qu   = Quad-spectrum, units of xy / unit of Fs

% Reference - Bendat, Julius and Piersol, Allan, Random Data, John Wiley
% 	      & Sons, Inc., 1986.

% Version 1.0 - 09/02/92
% Written by Harvey Chen
%%
psdy=[];
qu=[];
co=[];

% ensure even number of samples

if rem(length(x),2)~=0
  x=x(1:length(x)-1);
end

if y~=0 & rem(length(y),2)~=0
  y=y(1:length(y)-1);
end


% Set number of samples and sample rate

N=max(size(x));
SR=(Fs);         % fix(Fs); CHANGED 3/96 sm
tm=N/SR;
df=1/tm;

% Set Frequency axis

f=0:df:Fs;

x=x-mean(x);
y=y-mean(y);
x=x(:);
y=y(:);


% Calculation of Fourier transform and Power Spectral Density

xx=fft(x,N);
pxx=xx.*conj(xx)/N^2;
psdx=pxx/df;

% Size of psdx

di=max(size(psdx));

% Fold over spectrum
%disp('line 64 hspec_ss.m');
%pause

psdx=2.*psdx(2:di/2);
f=f(2:di/2);

% Compare Matlab calculated variance to variance calculated from PSD

varx=(std(x))^2;
var2x=sum (psdx)*df;

I=di/2-2;
itv=1;

% smoothing routine

clear itv
N2=max(size(psdx));
I=1;
check2=1;
check(1)=1;
	while check2<N2
	I=I+1;
	check2=check2+I;
	check(I)=check(I-1)+I;
	end
for q=1:I-1
	itv(q)=check(q+1)-check(q);
end

I=I-1;

check=1;

for i=1:I

if i~=I
psdx2(i)=sum (psdx(check:check+(i-1)))/(i);
f2(i)=sum (f(check:check+(i-1)))/(i);
end

if i==I
psdx2(i)=sum (psdx(check:N2))/(N2-check);
f2(i)=sum (f(check:N2))/(N2-check);
end

check=check+(i-1);

end

psdx=psdx2;
f=f2;

% end smoogthing

if y~=0

% Calculate y signal Power Spectral Density

yy=fft(y,N);
pyy=yy.*conj(yy)/N^2;
psdy=pyy/df;
psdy=2*psdy(2:di/2);

vary=(std(y))^2;
var2y=sum (psdy)*df;

% Calculate cross spectrum and separate into co and quad spectrums.

cross=xx.*conj(yy)/N^2;

cr2=cross/df;

co=2*real(cr2(2:di/2));
qu=2*imag(cr2(2:di/2));

% Smoothing Routine and Plots

%if ques2=='y'

check=1;
for i=1:I

if i~=I
psdy2(i)=sum (psdy(check:check+(i-1)))/(i);
co2(i)=sum (co(check:check+(i-1)))/(i);
qu2(i)=sum (qu(check:check+(i-1)))/(i);
end

if i==I
psdy2(i)=sum (psdy(check:N2))/(N2-check);
co2(i)=sum (co(check:N2))/(N2-check);
qu2(i)=sum (qu(check:N2))/(N2-check);
end

check=check+(i-1);

end

psdy=psdy2;
co=co2;
qu=qu2;

%end

f=f(:);
co=co(:);
qu=qu(:);


end

% Make sure each variable is a column vector and output the array.

   f=f(:);
   psdx=psdx(:);
if ~isempty(psdy);   psdy=psdy(:); end
if ~isempty(co);   co=co(:);end
if ~isempty(qu);   qu=qu(:);end



if y==0

   P=[ f psdx];

end

if y~=0
  
   P=[ f psdx psdy co qu];

end
