function [pxx,fxx,f] = cpsd(x,fs,nfft,fasesw)
%function [pxx,fxx,f] = cpsd(x,fs,nfft,fasesw)
% Estima 'nfft' puntos de la PSD de 'x' muestreada a 'fs' Hz
% cpsd(x,fs,nfft) Traza la PSD en escala logaritmica
% cpsd(x,fs,nfft,fasesw) Si fasesw==1 también se traza la fase
% [pxx,fxx,f] = cpsd(x,fs,nfft) Devuelve la PSD en 'pxx',
%   la fase en 'fxx' y el vector de frecuencias en 'f'
  if nargin < 3
    disp('Formato: cpsd(x,fs,nfft)');
  else
  if nargin == 3, fasesw = 0; end
% Cálculo de la PSD y FASE
    npts = length(x);
    w=hanning(nfft);
    h = fft(x(1:nfft).*w,nfft);
    psd = abs(h).^2;
    fase = unwrap(angle(h));
    n = nfft;
    np = 1;
    while n <= npts - nfft;
      h = fft(x(n+1:n+nfft).*w,nfft);
      psd = psd + abs(h).^2;
      fase = fase + unwrap(angle(h));
      n = n + nfft/2;
      np = np + 1;
    end
    psd  = psd(1:nfft/2+1) / (norm(w)^2 * np);
    fase = fase(1:nfft/2+1) / np;
    fvec = fs*(0:nfft/2)/nfft;

    if nargout == 0
      psd = 10*log10(psd*norm(w)^2/sum(w)^2);
      fase = fase / pi;
      if fasesw == 1
        subplot(2,1,1)
          plot(fvec,psd);
          set(gca,'XLim',[0,fs/2]);
          xlabel('Frecuencia (Hz)');
          ylabel('PSD (dB)');
        subplot(2,1,2);
          plot(fvec,fase);
          set(gca,'XLim',[0 fs/2]);
          xlabel('Frecuencia (Hz)');
          ylabel('Fase (xPI)');
      else
        plot(fvec,psd);
        set(gca,'XLim',[0,fs/2]);
        xlabel('Frecuencia (Hz)');
        ylabel('PSD (dB)');
      end
    else
      pxx = psd / length(psd);
      fxx = fase;
      f   = fvec;
    end
  end
end
