% MATLAB implementation of the PVModel
% Present implementation is text-based only
% Currently we can predict for
% 	* Stern, Trahiotis, and Tao paper
%	* Pure tones
% 	* AM tones

% Written by Richard Stern, with help from Gerald Ng

% This version localized triple-AM stimuli
% Last revised 3/19/97 - rms

% clear

global As f0 strflag threshold;

tic

% In the present version, the time increments for t and tau
% must be identical (represented by tinc)

% All times are in msec; all frequencies are in kHz
% NOTE: This version currently assumes periodicity along the t
% axis .... i.e. the total number of t values must represent
% an integer number of periods.

% ********* Initializing variables ***********

[finc,fcvec,numfc] = f_init;   % initialization of freq variables

[tinc,t,tau,numt,numtau] = t_init;  % intializing temporal vars.

[strnum,strflag] = str_init;   % initializing straightness vars.
%[strnum] = str_init;   % initializing straightness vars.

% [As,f0,fm,modindex,stimset,ts] = stim_ini;

g=sqrt(20);

% Initializing indices for Lttau, etc.

rrLttau = zeros(numt-1,numtau);
rlLttau = zeros(numt-1,numtau);
indexrl = zeros(numt-1,numtau);
indexrr = zeros(numt-1,numtau);

for m = 1:numtau,
  for n = 1:numt-1, %leave out the last sample
    indexrr(n,m) = 1+rem(round((2*t(n)+tau(m))/tinc),(numt-1));
    indexrl(n,m) = 1+rem(round((2*t(n)-tau(m))/tinc),(numt-1));

    if(indexrr(n,m) <= 0)
      indexrr(n,m) = indexrr(n,m)+numt-1;
      end
     if(indexrl(n,m) <= 0)
     indexrl(n,m) = indexrl(n,m)+numt-1;
     end
    end
  end



% Setting up sluggishness filter

freqslug = .005; % cutoff for sluggishness filter
h = exp(-2*pi*freqslug*t);
h = h(1,1:numt-1); % discard last point

% The characteristic frequency loop

% Generate characteristic frequencies to be used

Ltaufc = zeros(numfc,numtau);
tempLttaufc = zeros((numt-1),numtau);


% Call all stimulus-specific parameters here
[As,f0,fm,modindex,ts,stimset] = stim_init;


for i = 1:numfc,
  Lttaufc = zeros((numt-1),numtau);
  fc = fcvec(i)


   
  sttcomp;	% This is the subroutine that works on the stimulus-specific information.
  			% It is indended that stimulus/experiment-specific changes be primarily
			% in this routine (plus the initialization files).

  Lttaufc = rlLttau .* rrLttau;

% mesh(tau,t,rlLttau)
% mesh(tau,t,rrLttau)
% mesh(tau,t,Lttau)
% view(0,90)

% saving Lttaufc to files according to frequency blocks
% indexed from 1 to numfc

  eval(['Lttaufc' int2str(i) '=Lttaufc;']);
  eval(['save Lttaufc' int2str(i) ' Lttaufc' int2str(i) ';']);

end	% Primary frequency loop

% Straightness weighting
% Performed by adding in log domain
% This segment also performs temporal averaging
% by convolving with the "sluggishness filter"

% Lttaufc
if strflag > 0

firstf=ceil(strnum/2);            % edges of straightness frequencies
lastf=numfc-floor(strnum/2);
temptave = zeros(numt-1,numtau);   % initialization for temporal avg.

  for it = firstf:lastf,

    if it==firstf

      for j = 1:strnum,     % load initial Lttaufc's
        eval(['load Lttaufc' int2str(j) ';']);

        if j <= floor(strnum/2)  % perform sluggishness only for frequencies
          for k = 1:numtau,      % that don't need to be updated
            temptave(:,k) = ...
              (circconv((eval(['Lttaufc' int2str(j) '(:,k)']))',h))';
            end    % sluggishness convolution

          % plot(tau,Ltave(1,:))
          Ltaufc(j,:)=temptave(1,:);
          end   

        eval(['Lttaufc' int2str(j) '=log(Lttaufc' int2str(j) '+eps);']);
        end    % initial loading of frequency blocks
      end    % initial case

    eval(['tempLttaufc=Lttaufc' int2str(it) ';']);

    for l = 1:floor(strnum/2),
      eval(['tempLttaufc = tempLttaufc+Lttaufc' int2str(it+l) ';']);
      eval(['tempLttaufc = tempLttaufc+Lttaufc' int2str(it-l) ';']);
      end   %  straightness weighting loop
  
    tempLttaufc = exp(tempLttaufc);

    for k=1:numtau,            % temporal averaging for present fblock
      temptave(:,k) = (circconv(tempLttaufc(:,k)',h)');
      end        % sluggishness convolution

    % plot(tau,Ltave(1,:));
    Ltaufc(it,:)=temptave(1,:);

    eval(['clear Lttaufc' int2str(it-floor(strnum/2)) ';']); % save space
    if it < lastf
      eval(['load Lttaufc' int2str(it+1+floor(strnum/2)) ';']); % save space
      eval(['Lttaufc' int2str(it+1+floor(strnum/2)) '=log(Lttaufc' ...
        int2str(it+1+floor(strnum/2)) '+eps);']);
    else
      eval(['clear Lttaufc' int2str(it) ';']);
      for q = (lastf+1):numfc,     % temporal averaging for remaining blocks
        eval(['Lttaufc' int2str(q) '=exp(Lttaufc' int2str(q) ');']);
        for k=1:numtau,
          temptave(:,k) = ...
            (circconv((eval(['Lttaufc' int2str(q) '(:,k)']))',h))';
          end    % sluggishness convolution

        % plot(tau,Ltave(1,:));
        Ltaufc(q,:)=temptave(1,:);
        eval(['clear Lttaufc' int2str(q) ';']);
        end   % for loop to process remaining blocks
      end    % if else (reached the last frequency)
    end   % for loop to perform straightness operation
  
  else    % strflag == 0

    temptave = zeros(numt-1,numtau);   % initialization for temporal avg.

    for i = 1:numfc,               % sluggishness convolution
      eval(['load Lttaufc' int2str(i) ';']);
      for k=1:numtau,
        temptave(:,k) = (circconv((eval(['Lttaufc' int2str(i) '(:,k)']))',h))';
        end
      eval(['clear Lttaufc' int2str(i) ';']);
      Ltaufc(i,:)=temptave(1,:); 

      end

    end     %if strflag

ptaufc = ptau(tau,fcvec);
exp_p = zeros(1,numfc);
tau_2d = zeros(numfc,numtau);  % for 2d matrix multiplication
  for j=1:numfc,
    tau_2d(j,:)=tau;
    end

exp_p = sum((tau_2d.*Ltaufc.*ptaufc)')./sum((Ltaufc.*ptaufc)');
position = mean(exp_p);		% Un-normalized predicted lateral position 
[ts position]
  
toc

