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));