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)

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