function Y=surrogate(x,M) % % Y=surrogate(x,M) % % Generates M amplitude adjusted phase shuffled surrogate time series from x. % Useufel for testing the underlying assumption that the null hypothesis consists % of linear dynamics with possibly non-linear, monotonically increasing, % measurement function. % % Required Input Parameters: % % x % Nx1 vector of doubles % % M % 1x1 scalar specifying the number of surrogate time series to % generate. % % Required Output Parameters: % % Y % NxM vector of doubles % % % % References: % %[1] Kaplan, Daniel, and Leon Glass. Understanding nonlinear dynamics. Vol. 19. Springer, 1995. % % % Written by Ikaro Silva, 2014 % Last Modified: November 20, 2014 % Version 1.0 % Since 0.9.8 % % % % % See also MSENTROPY, SURROGATE %endOfHelp % 1. Amp transform original data to Gaussian distribution % 2. Phase randomize #1 % 3. Amp transform #2 to original % Auto-correlation function should be similar but not exact! x=x(:); N=length(x); Y=zeros(N,M); for m=1:M %Step 1 y=randn(N,1); y=amplitudeTransform(x,y,N); %Step 2 y=phaseShuffle(y,N); %Step 3 y=amplitudeTransform(y,x,N); Y(:,m)=y; end %%% Helper functions function target=amplitudeTransform(x,target,N) %Steps: %1. Sort the source by increasing amp %2. Sort target as #1 %3. Swap source amp by target amp %4. Sort #3 by increasing time index of #1 X=[[1:N]' x]; X=sortrows(X,2); target=[X(:,1) sort(target)]; target=sortrows(target,1); target=target(:,2); function y=phaseShuffle(x,N) %%Shuffle spectrum X=fft(x); Y=X; mid=floor(N/2)+ mod(N,2); phi=2*pi*rand(mid-1,1); %Generate random phase Y(2:mid)=abs(X(2:mid)).*cos(phi) + j*abs(X(2:mid)).*sin(phi); if(~mod(N,2)) %Even series has Nyquist in the middle+1 because of DC Y(mid+2:end)=conj(flipud(Y(2:mid))); Y(mid+1)=X(mid+1); else %Odd series is fully symetric except for DC Y(mid+1:end)=conj(flipud(Y(2:mid))); end y=real(ifft(Y));