Home > atmlab > math > movingstd.m

movingstd

PURPOSE ^

movingstd: efficient windowed standard deviation of a time series

SYNOPSIS ^

function s = movingstd(x,k,windowmode)

DESCRIPTION ^

 movingstd: efficient windowed standard deviation of a time series
 usage: s = movingstd(x,k,windowmode)

 Movingstd uses filter to compute the standard deviation, using
 the trick of std = sqrt((sum(x.^2) - n*xbar.^2)/(n-1)).
 Beware that this formula can suffer from numerical problems for
 data which is large in magnitude.

 At the ends of the series, when filter would generate spurious
 results otherwise, the standard deviations are corrected by
 the use of shorter window lengths.

 arguments: (input)
  x   - vector containing time series data

  k   - size of the moving window to use (see windowmode)
        All windowmodes adjust the window width near the ends of
        the series as necessary.

        k must be an integer, at least 1 for a 'central' window,
        and at least 2 for 'forward' or 'backward'

  windowmode - (OPTIONAL) flag, denotes the type of moving window used
        DEFAULT: 'central'

        windowmode = 'central' --> use a sliding window centered on
            each point in the series. The window will have total width
            of 2*k+1 points, thus k points on each side.
        
        windowmode = 'backward' --> use a sliding window that uses the
            current point and looks back over a total of k points.
        
        windowmode = 'forward' --> use a sliding window that uses the
            current point and looks forward over a total of k points.

        Any simple contraction of the above options is valid, even
        as short as a single character 'c', 'b', or 'f'. Case is
        ignored.

 arguments: (output)
  s   - vector containing the windowed standard deviation.
        length(s) == length(x)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

movingstd.m

SOURCE CODE ^

0001 function s = movingstd(x,k,windowmode)
0002 % movingstd: efficient windowed standard deviation of a time series
0003 % usage: s = movingstd(x,k,windowmode)
0004 %
0005 % Movingstd uses filter to compute the standard deviation, using
0006 % the trick of std = sqrt((sum(x.^2) - n*xbar.^2)/(n-1)).
0007 % Beware that this formula can suffer from numerical problems for
0008 % data which is large in magnitude.
0009 %
0010 % At the ends of the series, when filter would generate spurious
0011 % results otherwise, the standard deviations are corrected by
0012 % the use of shorter window lengths.
0013 %
0014 % arguments: (input)
0015 %  x   - vector containing time series data
0016 %
0017 %  k   - size of the moving window to use (see windowmode)
0018 %        All windowmodes adjust the window width near the ends of
0019 %        the series as necessary.
0020 %
0021 %        k must be an integer, at least 1 for a 'central' window,
0022 %        and at least 2 for 'forward' or 'backward'
0023 %
0024 %  windowmode - (OPTIONAL) flag, denotes the type of moving window used
0025 %        DEFAULT: 'central'
0026 %
0027 %        windowmode = 'central' --> use a sliding window centered on
0028 %            each point in the series. The window will have total width
0029 %            of 2*k+1 points, thus k points on each side.
0030 %
0031 %        windowmode = 'backward' --> use a sliding window that uses the
0032 %            current point and looks back over a total of k points.
0033 %
0034 %        windowmode = 'forward' --> use a sliding window that uses the
0035 %            current point and looks forward over a total of k points.
0036 %
0037 %        Any simple contraction of the above options is valid, even
0038 %        as short as a single character 'c', 'b', or 'f'. Case is
0039 %        ignored.
0040 %
0041 % arguments: (output)
0042 %  s   - vector containing the windowed standard deviation.
0043 %        length(s) == length(x)
0044 
0045 % check for a windowmode
0046 if (nargin<3) || isempty(windowmode)
0047   % supply the default:
0048   windowmode = 'central';
0049 elseif ~ischar(windowmode)
0050   error 'If supplied, windowmode must be a character flag.'
0051 end
0052 % check for a valid shortening.
0053 valid = {'central' 'forward' 'backward'};
0054 windowmode = lower(windowmode);
0055 ind = strmatch(windowmode,valid);
0056 if isempty(ind)
0057   error 'Windowmode must be a character flag: ''c'', ''b'', or ''f''.'
0058 else
0059   windowmode = valid{ind};
0060 end
0061 
0062 % length of the time series
0063 n = length(x);
0064 
0065 % check for valid k
0066 if (nargin<2) || isempty(k) || (rem(k,1)~=0)
0067   error 'k was not provided or not an integer.'
0068 end
0069 switch windowmode
0070   case 'central'
0071     if k<1
0072       error 'k must be at least 1 for windowmode = ''central''.'
0073     end
0074     if n<(2*k+1)
0075       error 'k is too large for this short of a series and this windowmode.'
0076     end
0077   otherwise
0078     if k<2
0079       error 'k must be at least 2 for windowmode = ''forward'' or ''backward''.'
0080     end
0081     if (n<k)
0082       error 'k is too large for this short of a series.'
0083     end
0084 end
0085 
0086 % Improve the numerical analysis by subtracting off the series mean
0087 % this has no effect on the standard deviation.
0088 x = x - mean(x);
0089 
0090 % we will need the squared elements
0091 x2 = x.^2;
0092 
0093 % split into the three windowmode cases for simplicity
0094 A = 1;
0095 switch windowmode
0096   case 'central'
0097     B = ones(1,2*k+1);
0098     s = sqrt((filter(B,A,x2) - (filter(B,A,x).^2)*(1/(2*k+1)))/(2*k));
0099     s(k:(n-k)) = s((2*k):end);
0100   case 'forward'
0101     B = ones(1,k);
0102     s = sqrt((filter(B,A,x2) - (filter(B,A,x).^2)*(1/k))/(k-1));
0103     s(1:(n-k+1)) = s(k:end);
0104   case 'backward'
0105     B = ones(1,k);
0106     s = sqrt((filter(B,A,x2) - (filter(B,A,x).^2)*(1/k))/(k-1));
0107 end
0108 
0109 % special case the ends as appropriate
0110 switch windowmode
0111   case 'central'
0112     % repairs are needed at both ends
0113     for i = 1:k
0114       s(i) = std(x(1:(k+i)));
0115       s(n-k+i) = std(x((n-2*k+i):n));
0116     end
0117   case 'forward'
0118     % the last k elements must be repaired
0119     for i = (k-1):-1:1
0120       s(n-i+1) = std(x((n-i+1):n));
0121     end
0122   case 'backward'
0123     % the first k elements must be repaired
0124     for i = 1:(k-1)
0125       s(i) = std(x(1:i));
0126     end
0127 end
0128

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005