% function P = get_rand_seq(n_seq,n_trials,reward_range,corr_tol,fun,include_fun)
%
% This function can be used to create multiple sequences that look like a
% random walk, but that
% (a) pass through the whole range rather than doing so by chance
% and
% (b) all have the same frequency spectrum (i.e. the same rate of change).
% The correlation range can be set to a specific value denoting how much
% these sequences are allowed to correlate.
% Please check the sequences, because just by chance specific unwanted
% patterns can emerge (such as random fluctuation with a specific shape).
%
% Please note: This function only checks for correlations across trials for
% all frequencies together, not for correlations across trials for a given
% frequency. This means it is possible that specific frequencies alone are
% highly correlated. For example, an observer may not perceive specific
% changes, but only very slow changes (equivalent to a low-pass filter).
% In this case, the sequences may be perceived as being correlated although
% they are not.
%
% Input:
% n_seq: Number of sequences
% n_trials: Number of trials for each sequence (I'd recommend creating
% shorter blocks that all have the same sequence length, e.g. when
% you have 200 trials, maybe create four sequence pairs of 50 trials)
% reward_range: This is a 1x2 vector denoting the range of rewards we
% will get later for our sequence (e.g. [0.2 0.8]). Please note that
% we use Fourier transformation, i.e. that these ranges might be
% exceeded. We control for that by slightly reducing the input range.
% If you care about exactly reaching the extreme values, then
% normalize your results (which will however influence the frequency
% spectrum). (default: [0.05 0.95])
% corr_tol: This is a 1x2 vector denoting the range in which the
% sequences are allowed to be correlated (default: [-0.05 0.05]).
% fun (optional): Per default, all values within reward_range are equally
% likely (i.e. a linear function). If you want to change that, you
% can manually enter the function yourself. In this case,
% reward_range and n_trials will not play a role.
% [default: linspace(reward_range(1),reward_range(2),n_trials) ]
% include_fun (optional): 0 or 1. Whether to include the sequence in fun (default: 0)
%
% Example:
% P = get_rand_seq(2,64,[0.2 0.8])
% plot(P), axis([1 64 0 1])
% % check power spectrum difference
% pw1 = abs(fft(P(:,1))); pw2 = abs(fft(P(:,2)));
% % Should be close to zero
% range(pw1-pw2)
%
% 2014/03 Martin Hebart
function P = get_rand_seq(n_seq,n_trials,reward_range,corr_tol,fun,include_fun)
max_iter = 100000; % number of iterations for search (per sequence)
% set defaults
if ~exist('reward_range','var') || isempty(reward_range)
reward_range = [0.05 0.95];
end
if ~exist('corr_tol','var') || isempty(corr_tol)
corr_tol = [-0.05 0.05];
end
if ~exist('include_fun','var') || isempty(include_fun)
include_fun = 0;
end
run_checks; % see function at end
% adjust reward range
m = mean(reward_range);
reward_range = 0.95*(reward_range-m)+m;
if ~exist('fun','var') || isempty(fun)
x = linspace(reward_range(1),reward_range(2),n_trials)'; % default value
else
x = fun;
end
% run fourier transformation, get power and phase
fftx = fft(x);
powerx = abs(fftx);
phasex = angle(fftx);
if include_fun % should we include the original function?
P(:,1) = fun;
else
% carry out phase scrambling and include scrambled into output
randphase = angle(fft(rand(size(x))));
% this below is the same as above, but shifting all phases by +/- pi/2 providing us with uncorrelated sequences within each frequency band
% randphase(1) = phasex(1); randphase(2:end) = pi*round(rand(size(randphase(2:end))))-pi/2;
scrphase = phasex + randphase;
P(:,1) = real(ifft(powerx.*exp(sqrt(-1)*(scrphase))));
end
for i_seq = 2:n_seq
ct = 0; % counter for bailing out
ind = tril(true(i_seq),-1); % index for comparison of correlations
while 1
ct = ct+1;
% do same as above, keep power
randphase = angle(fft(rand(size(x))));
% randphase(1) = phasex(1); randphase(2:end) = pi*round(rand(size(randphase(2:end))))-pi/2;
scrphase = phasex + randphase;
p = real(ifft(powerx.*exp(sqrt(-1)*(scrphase))));
% correlate all sequences
r = corr([P p]);
% check if all is in right range, otherwise repeat
if all(r(ind)>corr_tol(1)) && all(r(ind) % update P
break
end
% check if number of maximum iterations is exceeded
if ct > max_iter
fprintf('With given parameters couldn''t find %i sequences that match.\nPlease adjust parameters or edit function and change variable max_iter!\n',n_seq)
P = [];
return
end
end
end
if any(P(:)>1)
warning('GET_RAND_SEQ:LARGER1','Some entries were larger than 1 and will be truncated to 0. This may slightly affect the frequency, mean, and std')
P(P(:)>1) = 1;
end
if any(P(:)<0)
warning('GET_RAND_SEQ:SMALLER0','Some entries were smaller than 0 and will be truncated to 0. This may slightly affect the frequency, mean, and std')
P(P(:)<0) = 0;
end
function run_checks
if n_seq < 2
error('Minimum number of sequences is 2')
end
if any(reward_range<0|reward_range>1)
error('Variable reward_range must be in range 0 to 1')
end
if reward_range(1) == reward_range(2)
error('Range of variable reward_range must be larger than 0')
end
if reward_range(1)>reward_range(2)
reward_range = reward_range([2 1]);
end
if corr_tol(1)>corr_tol(2)
corr_tol = corr_tol([2 1]);
end
if exist('fun','var') && ~isempty(fun)
if size(fun,1)==1
fun = fun';
end
if ~isempty(n_trials) && n_trials ~= length(fun)
warning('RUN_CHECKS:N_TRIALS_IGNORED','Input variable ''n_trials'' ignored. Using length of ''fun'' as number of trials.')
end
if reward_range(1) ~= min(fun) || reward_range(2) ~= max(fun)
warning('RUN_CHECKS:REWARD_RANGE_IGNORED','Input variable ''reward_range'' ignored. Using range of ''fun'' as range.')
end
end
end
end