--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/welchscale.m	Tue Dec 06 11:09:25 2011 +0100
@@ -0,0 +1,54 @@
+% WELCHSCALE scales the output of welch to be in the required units
+%
+% $Id: welchscale.m,v 1.1 2011/12/01 09:41:22 hewitson Exp $
+%
+
+function [yy, dyy, info] = welchscale(xx, dxx, win, fs, norm, inunits)
+  
+  nfft = length(win);
+  S1   = sum(win);
+  S2   = sum(win.^2);
+  enbw = fs * S2 / (S1*S1);
+  
+  if isempty(norm)
+    norm = 'None';
+  end
+  switch lower(norm)
+    case 'asd'
+      yy = sqrt(xx);
+      if isempty(dxx)
+        dyy = dxx;
+      else
+        dyy = 1./2./sqrt(xx) .* dxx;
+      end
+      info.units = inunits ./ unit('Hz^0.5');
+    case 'psd'
+      yy = xx;
+      dyy = dxx;
+      info.units = inunits.^2/unit('Hz');
+    case 'as'
+      yy = sqrt(xx * enbw);
+      if isempty(dxx)
+        dyy = dxx;
+      else
+        dyy = 1./2./sqrt(xx) .* dxx * enbw;
+      end
+      info.units = inunits;
+    case 'ps'
+      yy = xx * enbw;
+      dyy = dxx * enbw;
+      info.units = inunits.^2;
+    case 'none'
+      yy = xx;
+      dyy = dxx;
+      info.units = inunits;
+    otherwise
+      error('Unknown normalisation');
+  end
+  
+  info.nfft = nfft;
+  info.enbw = enbw;
+  info.norm = norm;
+  
+end
+