function [S, max_i, fluxdir, S3tr] = strfnc( trace, freq, maxlag ) %STRFNC Structure function calculation (following Van Atta, 1977) % Jason Kelley NEWAg Lab OSU % Last modified 09mar16 % Questions? Contact % % INPUTS 'trace' data to analyze (Nx1 vector array) % 'freq' sampling frequency (Hz) % 'maxlag' maximum lag time, (seconds) % % OUTPUTS 'S' structure functions S(r)^n and lag r (in seconds) for each row % only calculates 2 3 5 order to save memory, % column order corresponds to SF order, also calculated -S^3(r)/r % format: [r S^2(r) S^3(r) -S^3(r)/r S^5(r)] % 'S3tr' sinusoidal trace of 3rd structure function, used in spectral ramp detection m = length(trace); lags = 1:maxlag*freq; rn = length(lags); S = zeros(rn,5); % old method %{ tic for j = lags r = lags(j); early = trace(1:end-r); later = trace(r+1:end); diffs = later-early; for i = [2 3 5] S(j,i) = sum((diffs).^i)/(m-r); end %structure functions at lag j end % lags j S(:,6) = lags./freq; time = toc; fprintf('loop takes %f s\n', time); %} %tic % begin new method % dev %{ % % for f = lags % filt{lags(f)} = [-1; zeros(lags(f)-1,1); 1]; % end % for f = lags % TC(f,:) = conv(trace,filt{f},'same'); % end % Sc = [sum(TC.^1,) sum(TC.^2) sum(TC.^3) sum(TC.^4) sum(TC.^5)]; % Sc = Sc./repmat(m-lags,1,5); % Sc(:,6) = lags./freq; Sc = zeros(rn,6); %} % lagged differences by convolution filt = [ones(1,rn); -eye(rn)]; % singleton comparators at 1:rn lags e.g. [1 0 0 -1] cT = conv2(trace,filt); % conv filter with trace to get all lags cT = cT(rn+1:end-rn,:); % trim edges. 'same' does not work here for conv2 %time = toc; fprintf('conv takes %f s\n', time); % dev %{ % Sc(:,2) = sum(power(cT,ones(m-rn,rn).*2),1); % time = toc; fprintf('square and sum takes %f s\n', time); % %Sc(:,3) = sum(power(cT,ones(m-rn,rn).*3),1)./(m-(1:rn)); % Sc(:,3) = sum(cT,1)'.^Sc(:,2); % time = toc; fprintf('cumulative sum takes %f s\n', time); % Sc(:,5) = sum(power(cT,ones(m-rn,rn).*5),1)./(m-(1:rn)); % time = toc; fprintf('square and sum takes %f s\n', time); %} %tic cTp(:,:,1) = power(cT,ones(m-rn,rn).*2); % for second order SF cTp(:,:,2) = cTp(:,:,1).*cT; % third order SF cTp(:,:,3) = cTp(:,:,1).*cTp(:,:,2); % fifth order SF w = (m-rn-(1:rn))-1; % unbiased weighting vector 1/(N-1) S(:,2) = sum(cTp(:,:,1),1)./w; % to simplify, column order corresponds to SF order S(:,3) = sum(cTp(:,:,2),1)./w; % i.e. 3rd order SF is S(:,3) S(:,5) = sum(cTp(:,:,3),1)./w; %time = toc; fprintf('square and sum takes %f s\n', time); S(:,1) = lags./freq; % sample N -> dt S(:,4) = -S(:,3)./S(:,1); % ratio used to detect maximum contributing ramps % identify time lag at which S_3(r)/r is maximized, account for flux direction by sign of S^3(r)/r fluxdir = sign(nanmean(S(:,4),1)); [~,max_i] = max(fluxdir.*(S(:,4))); %max_r = S(max_i,1); % actually it is better to output max_i since lookup can occur for max_r % optional output of 3rd order S.F. trace to use in spectral time lag characterization if nargout>3 %[max3sf(k),Sii(k)] = findpeaks(S(:,4),'NPEAKS',1); S3tr = cTp(:,max_i,2); end end %function