function [Spec,f,SpecConf] = spectrum2(varargin)
%SPECTRUM Power spectrum estimate of one or two data sequences.
%   P=SPECTRUM2(X,NFFT,NOVERLAP,WIND) estimates the Power Spectral Density of
%   signal vector(s) X using Welch's averaged periodogram method. X must be a
%   set of column vectors. They are divided into overlapping sections, each of 
%   which is detrended and windowed by the WINDOW parameter, then zero padded 
%   to length NFFT. The magnitude squared of the length NFFT DFTs of the sec-
%   tions are averaged to form Pxx. P is a two times SIZE(X,2) column matrix 
%   P = [Pxx Pxxc]; the second half of columns Pxxc are the 95% confidence 
%   intervals. The number of rows of P is NFFT/2+1 for NFFT even, (NFFT+1)/2 
%   for NFFT odd, or NFFT if the signal X is complex. If you specify a scalar 
%   for WINDOW, a Hanning window of that length is used.
%
%   [P,F] = SPECTRUM2(X,NFFT,NOVERLAP,WINDOW,Fs) given a sampling frequency
%   Fs returns a vector of frequencies the same length as Pxx at which the 
%   PSD is estimated.  PLOT(F,P(:,1:end/2)) plots the power spectrum estimate 
%   versus true frequency.
%
%   [P, F] = SPECTRUM2(X,NFFT,NOVERLAP,WINDOW,Fs,Pr) where Pr is a scalar 
%   between 0 and 1, overrides the default 95% confidence interval and 
%   returns the Pr*100% confidence interval for Pxx instead.
%
%   SPECTRUM2(X) with no output arguments plots the PSD in the current
%   figure window, with confidence intervals.
%
%   The default values for the parameters are NFFT = 256 (or SIZE(X,1),
%   whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2, 
%   and Pr = .95. You can obtain a default parameter by leaving it out
%   or inserting an empty matrix [], e.g. SPECTRUM(X,[],128).
%
%   P = SPECTRUM2(X,Y) performs spectral analysis of the two times SIZE(X,2)
%   sequences X and Y using the Welch method. SPECTRUM returns the 8 times 
%   SIZE(X,2)column array
%      P = [Pxx Pyy Pxy Txy Cxy Pxxc Pyyc Pxyc]
%   where
%      Pxx  = X-vector power spectral density
%      Pyy  = Y-vector power spectral density
%      Pxy  = Cross spectral density
%      Txy  = Complex transfer function from X to Y = Pxy./Pxx
%      Cxy  = Coherence function between X and Y = (abs(Pxy).^2)./(Pxx.*Pyy)
%      Pxxc,Pyyc,Pxyc = Confidence range.
%   All input and output options are otherwise exactly the same as for the
%   single input case.
%   
%   SPECTRUM2(X,Y) with no output arguments will plot Pxx, Pyy, abs(Txy),
%   angle(Txy) and Cxy in sequence, pausing between plots.
%
%   SPECTRUM2(X,...,DFLAG), where DFLAG can be 'linear', 'mean' or 'none', 
%   specifies a detrending mode for the prewindowed sections of X (and Y).
%   DFLAG can take the place of any parameter in the parameter list
%   (besides X) as long as it is last, e.g. SPECTRUM(X,'none');
%
%   See also SPECTRUM, PSD, CSD, TFE, COHERE, SPECGRAM, SPECPLOT, DETREND, 
%     PMTM, PMUSIC.
%   ETFE, SPA, and ARX in the Identification Toolbox.

%   The units on the power spectra Pxx and Pyy are such that, using
%   Parseval's theorem: 
%
%        SUM(Pxx)/LENGTH(Pxx) = SUM(X.^2)/size(x,1) = COV(X)
%
%   The RMS value of the signal is the square root of this.
%   If the input signal is in Volts as a function of time, then
%   the units on Pxx are Volts^2*seconds = Volt^2/Hz.
%
%   Here are the covariance, RMS, and spectral amplitude values of
%   some common functions:
%         Function   Cov=SUM(Pxx)/LENGTH(Pxx)   RMS        Pxx
%         a*sin(w*t)        a^2/2            a/sqrt(2)   a^2*LENGTH(Pxx)/4
%Normal:  a*rand(t)         a^2              a           a^2
%Uniform: a*rand(t)         a^2/12           a/sqrt(12)  a^2/12
%   
%   For example, a pure sine wave with amplitude A has an RMS value
%   of A/sqrt(2), so A = SQRT(2*SUM(Pxx)/LENGTH(Pxx)).
%
%   See Page 556, A.V. Oppenheim and R.W. Schafer, Digital Signal
%   Processing, Prentice-Hall, 1975.

error(nargchk(1,8,nargin))
[msg,x,y,nfft,noverlap,window,Fs,p,dflag]=specchk2(varargin);
error(msg)


if isempty(p),
   p = .95;   % default confidence interval even if not asked for
end

n = size(x,1);	  % Number of data points
ns = size(x,2);   % Number of signals
nwind = length(window);
if n < nwind    % zero-pad x (and y) if length less than the window length
   x(nwind,:)=0;  n=nwind;
   if ~isempty(y), y(nwind)=0;  end
end

k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
index = 1:nwind;
KMU = k*norm(window)^2; % Normalizing scale factor ==> asymptotically unbiased
% KMU = k*sum(window)^2;% alt. Nrmlzng scale factor ==> peaks are about right

if (isempty(y)) % Single sequence case.
   Pxx = zeros(nfft,ns); Pxx2 = zeros(nfft,ns);
   for i=1:k
      for l=1:ns
         if strcmp(dflag,'linear')
            xw = window.*detrend(x(index,l));
         elseif strcmp(dflag,'none')
            xw = window.*(x(index,l));
         else
            xw = window.*detrend(x(index,l),0);
         end
         Xx = abs(fft(xw,nfft)).^2;
         Pxx(:,l) = Pxx(:,l) + Xx;
         Pxx2(:,l) = Pxx2(:,l) + abs(Xx).^2;
      end
      index = index + (nwind - noverlap);
   end
   % Select first half
   if ~any(any(imag(x)~=0)),   % if x and y are not complex
      if rem(nfft,2),    % nfft odd
         select = [1:(nfft+1)/2];
      else
         select = [1:nfft/2+1];   % include DC AND Nyquist
      end
   else
     select = 1:nfft;
   end
   Pxx = Pxx(select,:);
   Pxx2 = Pxx2(select,:);
   cPxx = zeros(size(Pxx));
   if k > 1
      c = (k.*Pxx2-abs(Pxx).^2)./(k-1);
      c = max(c,zeros(size(Pxx)));
      cPxx = sqrt(c);
   end
   ff = sqrt(2)*erfinv(p);  % Equal-tails.
   Pxxc = ff.*cPxx/KMU;
   P = Pxx/KMU;
   Pc = Pxxc;
else
   Pxx = zeros(nfft,ns); % Dual sequence case.
   Pxy = Pxx; Pxx2 = Pxx; Pxy2 = Pxx;
   Pyy = zeros(nfft,1); Pyy2 = Pyy; 
   for i=1:k
      if strcmp(dflag,'linear')
         yw = window.*detrend(y(index));
      elseif strcmp(dflag,'none')
         yw = window.*(y(index));
      else
         yw = window.*detrend(y(index),0);
      end
      Yy = fft(yw,nfft);
      Yy2 = abs(Yy).^2;
      Pyy = Pyy + Yy2;
      Pyy2 = Pyy2 + abs(Yy2).^2;
      for l=1:ns
         if strcmp(dflag,'linear')
            xw = window.*detrend(x(index,l));
         elseif strcmp(dflag,'none')
            xw = window.*(x(index,l));
         else
            xw = window.*detrend(x(index,l),0);
         end
         Xx = fft(xw,nfft);
         Xx2 = abs(Xx).^2;
         Pxx(:,l) = Pxx(:,l) + Xx2;
         Pxx2(:,l) = Pxx2(:,l) + abs(Xx2).^2;
         Xy  = Yy .* conj(Xx);
         Pxy(:,l) = Pxy(:,l) + Xy;
         Pxy2(:,l) = Pxy2(:,l) + Xy .* conj(Xy);
      end
      index = index + (nwind - noverlap);
   end
   % Select first half
   if ~any(any(imag([x y])~=0)),   % if x and y are not complex
      if rem(nfft,2),    % nfft odd
         select = [1:(nfft+1)/2];
      else
         select = [1:nfft/2+1];   % include DC AND Nyquist
      end
   else
      select = 1:nfft;
   end
   Pxx = Pxx(select,:);
   Pxy = Pxy(select,:);
   Pxx2 = Pxx2(select,:);
   Pxy2 = Pxy2(select,:);
   Pyy = Pyy(select);
   Pyy2 = Pyy2(select);
   cPxx = zeros(size(Pxx));
   cPyy = zeros(size(Pyy));
   cPxy = cPxx;
   if k > 1
      c = max((k.*Pxx2-abs(Pxx).^2)./(k-1),zeros(size(Pxx)));
      cPxx = sqrt(c);
      c = max((k.*Pyy2-abs(Pyy).^2)./(k-1),zeros(size(Pyy)));
      cPyy = sqrt(c);
      c = max((k.*Pxy2-abs(Pxy).^2)./(k-1),zeros(size(Pxx)));
      cPxy = sqrt(c);
   end
   Txy = Pxy./Pxx;
   Cxy = (abs(Pxy).^2)./(Pxx.*repmat(Pyy,1,size(Pxx,2)));
   ff = sqrt(2)*erfinv(p);  % Equal-tails.
   Pxx = Pxx/KMU;
   Pyy = Pyy/KMU;
   Pxy = Pxy/KMU;
   Pxxc = ff.*cPxx/KMU;
   Pxyc = ff.*cPxy/KMU;
   Pyyc = ff.*cPyy/KMU;
   P = [Pxx Pyy Pxy Txy Cxy];
   Pc = [Pxxc Pyyc Pxyc];
end
freq_vector = (select - 1)'*Fs/nfft;
   
if nargout == 0,   % do plots
   newplot;
   if Fs==2, xl='Frequency'; else, xl = 'f  [Hz]'; end
   nplot=1+(1-isempty(y))*4;
   
   subplot(nplot,1,1);
   c = [max(Pxx-Pxxc,0)  Pxx+Pxxc];
   c = c.*(c>0);
 
   h=semilogy(freq_vector,Pxx,...
	      freq_vector,c(:,1:size(c,2)/2),'--',...
	      freq_vector,c(:,size(c,2)/2+1:end),'--');
   title('\bf X Power Spectral Density')
   ylabel('P_x')
   xlabel(xl)
   if length(h)>3
      s={};
      for k=1:length(h)/3
         c=get(h(k),'Color');
         set(h(k+length(h)/3),'Color',c);
         set(h(k+2*length(h)/3),'Color',c);
	 s{k}=['x_' num2str(k)];
      end
      legend(h(1:length(h)/3),s);
   end
   if (isempty(y)),   % single sequence case
      return
   end
   
   subplot(nplot,1,2);
   c = [max(Pyy-Pyyc,0)  Pyy+Pyyc];
   c = c.*(c>0);
   h=semilogy(freq_vector,Pyy,...
	      freq_vector,c(:,1),'--',...
	      freq_vector,c(:,2),'--');
   if size(Pxx,2)>1
      for k=1:length(h)/3
         c=get(h(k),'Color');
         set(h(k+length(h)/3),'Color',c);
         set(h(k+2*length(h)/3),'Color',c);
      end
   end
   
   title('\bf Y Power Spectral Density')
   ylabel('P_y')
   xlabel(xl)
   
   subplot(nplot,1,3);
   semilogy(freq_vector,abs(Txy));
   title('\bf Transfer function magnitude')
   ylabel('T_{xy}^{(m)}')
   xlabel(xl)

   subplot(nplot,1,4);
   plot(freq_vector,(angle(Txy))), ...
   title('\bf Transfer function phase')
   ylabel('T_{xy}^{(p)}')
   xlabel(xl)
   
   subplot(nplot,1,5);
   plot(freq_vector,Cxy);
   title('\bf Coherence')
   ylabel('C_{xy}')
   xlabel(xl)
   if exist ('niceaxes') ~= 0
      niceaxes(findall(gcf,'Type','axes'));
   end
elseif nargout ==1, 
   Spec = P;
elseif nargout ==2, 
   Spec = P;
   f = freq_vector;
elseif nargout ==3, 
   Spec = P;
   f = freq_vector;
   SpecConf = Pc;
end

function [msg,x,y,nfft,noverlap,window,Fs,p,dflag] = specchk2(P)
%SPECCHK Helper function for SPECTRUM
%   SPECCHK(P) takes the cell array P and uses each cell as 
%   an input argument.  Assumes P has between 1 and 7 elements.

msg = [];
if size(P{1},1)<=1
   if max(size(P{1}))==1
      msg = 'Input data must be a vector, not a scalar.';
   else
      msg = 'Requires column vector input.';    
   end
   x=[];
   y=[];
elseif (length(P)>1),
   if (all(size(P{1},1)==size(P{2})) & (size(P{1},1)>1) ) | ...
         size(P{2},1)>1,   % 0ne signal or 2 present?
      % two signals, x and y, present
      x = P{1}; y = P{2}; 
      % shift parameters one left
      P(1) = [];
   else 
      % only one signal, x, present
      x = P{1}; y = []; 
   end
else  % length(P) == 1
   % only one signal, x, present
   x = P{1}; y = []; 
end

% now x and y are defined; let's get the rest

if length(P) == 1
   nfft = min(size(x,1),256);
   window = hanning(nfft);
   noverlap = 0;
   Fs = 2;
   p = [];
   dflag = 'linear';
elseif length(P) == 2
   if isempty(P{2}),   dflag = 'linear'; nfft = min(size(x,1),256); 
   elseif isstr(P{2}), dflag = P{2};     nfft = min(size(x,1),256); 
   else                dflag = 'linear'; nfft = P{2};   end
   window = hanning(nfft);
   noverlap = 0;
   Fs = 2;
   p = [];
elseif length(P) == 3
   if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2};     end
   if isempty(P{3}),   dflag = 'linear'; noverlap = 0;
   elseif isstr(P{3}), dflag = P{3};     noverlap = 0;
   else                dflag = 'linear'; noverlap = P{3}; end
   window = hanning(nfft);
   Fs = 2;
   p = [];
elseif length(P) == 4
   if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2};     end
   if isstr(P{4})
      dflag = P{4};
      window = hanning(nfft);
   else
      dflag = 'linear';
      window = P{4};  window = window(:);   % force window to be a column
      if length(window) == 1, window = hanning(window); end
      if isempty(window), window = hanning(nfft); end
   end
   if isempty(P{3}), noverlap = 0;  else noverlap=P{3}; end
   Fs = 2;
   p = [];
elseif length(P) == 5
   if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2};     end
   window = P{4};  window = window(:);   % force window to be a column
   if length(window) == 1, window = hanning(window); end
   if isempty(window), window = hanning(nfft); end
   if isempty(P{3}), noverlap = 0;  else noverlap=P{3}; end
   if isstr(P{5})
      dflag = P{5};
      Fs = 2;
   else
      dflag = 'linear';
      if isempty(P{5}), Fs = 2; else Fs = P{5}; end
   end
   p = [];
elseif length(P) == 6
   if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2};     end
   window = P{4};  window = window(:);   % force window to be a column
   if length(window) == 1, window = hanning(window); end
   if isempty(window), window = hanning(nfft); end
   if isempty(P{3}), noverlap = 0;  else noverlap=P{3}; end
   if isempty(P{5}), Fs = 2;     else    Fs = P{5}; end
   if isstr(P{6})
      dflag = P{6};
      p = [];
   else
      dflag = 'linear';
      if isempty(P{6}), p = .95;    else    p = P{6}; end
   end
elseif length(P) == 7
   if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2};     end
   window = P{4};  window = window(:);   % force window to be a column
   if length(window) == 1, window = hanning(window); end
   if isempty(window), window = hanning(nfft); end
   if isempty(P{3}), noverlap = 0;  else noverlap=P{3}; end
   if isempty(P{5}), Fs = 2;     else    Fs = P{5}; end
   if isempty(P{6}), p = .95;    else    p = P{6}; end
   if isstr(P{7})
      dflag = P{7};
   else
      msg = 'DFLAG parameter must be a string.'; return
   end
end

% NOW do error checking
if isempty(msg)
   if (nfft<length(window)), 
      msg = 'Requires window''s length to be no greater than the FFT length.';
   end
   if (noverlap >= length(window)),
      msg = 'Requires NOVERLAP to be strictly less than the window length.';
   end
   if (nfft ~= abs(round(nfft)))|(noverlap ~= abs(round(noverlap))),
      msg = 'Requires positive integer values for NFFT and NOVERLAP.';
   end
   if ~isempty(p),
      if (prod(size(p))>1)|(p(1,1)>1)|(p(1,1)<0),
         msg = 'Requires confidence parameter to be a scalar between 0 and 1.';
      end
   end
   if (min(size(y))~=1)&(~isempty(y)),
      msg = 'Requires column vector input as second signal.';
   end
   if (size(x,1)~=length(y))&(~isempty(y)),
      msg = 'Requires X and Y to have the same number of rows.';
   end
end