function varargout = computepsd(Sxx,Svxx,w,range,nfft,Fs,esttype)
  %
  %  Slight modification of original MATLAB's computepsd to include correct
  %  scaling for the variance, i.e var(a*x) = a^2*var(x)
  %
  % VERSION:    $Id: computepsd.m,v 1.1 2009/07/22 15:29:59 miquel Exp $
  %
  % HISTORY:    21-07-09 M Nofrarias
  %                  Creation
  %
  %%%%%%%%%%%%%%%%%
  %
  % COMPUTEPSD  Compute the one-sided or two-sided PSD or Mean-Square.
  %   [Pxx,Pvxx,W,UNITS] = COMPUTEPSD(Sxx,Svxx,W,RANGE,NFFT,Fs,ESTTYPE) where the
  %   inputs and outputs are:
  %
  %   Inputs:
  %    Sxx   - Whole power spectrum [Power]; it can be a vector or a matrix.
  %            For matrices the operation is applied to each column.
  %    W     - Frequency vector in rad/sample or in Hz.
  %    RANGE - Determines if a 'onesided' or a 'twosided' Pxx and Sxx are
  %            returned.
  %    NFFT  - Number of frequency points.
  %    Fs    - Sampling Frequency.
  %    ESTTYPE - A string indicating the estimate type: 'psd', or 'ms' value.
  %
  %   Outputs:
  %    Pxx   - One-sided or two-sided PSD or MEAN-SQUARE (not scaled by Fs)
  %            depending on the input arguments RANGE and TYPE.
  %    W     - Frequency vector 0 to 2*Nyquist or 0 to Nyquist depending on
  %            range, units will be either rad/sample (if Fs is empty) or Hz
  %            (otherwise).
  %    UNITS - Either 'rad/sample' or 'Hz'
  %   Author(s): R. Losada
  %   Copyright 1988-2005 The MathWorks, Inc.
  %   $Revision: 1.1 $  $Date: 2009/07/22 15:29:59 $
 
  if nargin < 7
    iswinfs = false;
    if nargin < 6,
      esttype = 'psd';
    end
  end
  
  w = w(:); % Make sure we always returns a column vector for frequency
  
  % Generate the one-sided spectrum [Power] if so wanted
  if strcmp(range,'onesided'),
    if rem(nfft,2),
      select = 1:(nfft+1)/2;  % ODD
      Sxx_unscaled = Sxx(select,:); % Take only [0,pi] or [0,pi)
      Svxx_unscaled = Svxx(select,:); % Take only [0,pi] or [0,pi)
      Sxx = [Sxx_unscaled(1,:); 2*Sxx_unscaled(2:end,:)];  % Only DC is a unique point and doesn't get doubled
      Svxx = [Svxx_unscaled(1,:); 4*Svxx_unscaled(2:end,:)];     
    else
      select = 1:nfft/2+1;    % EVEN
      Sxx_unscaled = Sxx(select,:); % Take only [0,pi] or [0,pi)
      Svxx_unscaled = Svxx(select,:); % Take only [0,pi] or [0,pi)
      Sxx = [Sxx_unscaled(1,:); 2*Sxx_unscaled(2:end-1,:); Sxx_unscaled(end,:)]; % Don't double unique Nyquist point
      Svxx = [Svxx_unscaled(1,:); 4*Svxx_unscaled(2:end-1,:); Svxx_unscaled(end,:)]; % Don't double unique Nyquist point
    end
    w = w(select);
  end
  
  % Compute the PSD [Power/freq]
  if ~isempty(Fs),
    Pxx = Sxx./Fs; % Scale by the sampling frequency to obtain the psd
    Pvxx = Svxx./Fs^2; 
    units = 'Hz';
  else
    Pxx = Sxx./(2.*pi); % Scale the power spectrum by 2*pi to obtain the psd
    Pvxx = Svxx./(2.*pi)^2;
    units = 'rad/sample';
  end
  
  if strcmpi(esttype,'ms'),
    varargout = {Sxx,Svxx,w,units};  % Mean-square
  else
    varargout = {Pxx,Pvxx,w,units};  % PSD
  end
  
  % [EOF] computepsd.m