function Y=shuffle(fileName)
%
% Generates 50 surrogate time series from series in fileName
% by with AmplitudeAdjustedPhaseShuffle
%Underlying assumption is that the Null Hypothesis consists
%of linear dynamics with possibly non-linear, monotonically increasing,
%measurement function.
%
%  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=load(fileName);
x=x(:);
N=length(x);
M=20;

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

    %Save to file under surr_m
    save(['surr_' num2str(m)],'y')
end




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