--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/computepsd.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,82 @@
+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