Home > fcn > calculateInverse.m

calculateInverse

PURPOSE ^

SYNOPSIS ^

function [out]=calculateInverse(motor,mode,mode_error,mode_dim)

DESCRIPTION ^

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [out]=calculateInverse(motor,mode,mode_error,mode_dim)  
0002   N=size(mode_error,2);
0003   InverseParameters=DIVA('MotorCortex','InverseParameters');
0004   
0005   for ninput=1:N,                  
0006     % use grid* notation for linearization of
0007     % articulatory/modal trajectories.  Form Jacobian from
0008     % linear forms and perform inverse operation on Jacobian
0009     
0010     nMotor = size(motor,1);
0011     nMode  = size(mode,1);
0012     
0013     current_motor = motor(:,ninput);
0014     current_mode  = mode(:,ninput);
0015     
0016     grid_motor = current_motor*ones(1,nMotor)+...
0017         InverseParameters.DELTA*diag(ones(1,nMotor));
0018     grid_mode   = current_mode(:,ones(1,nMotor));
0019 
0020     % Determine the Jacobian
0021     if(nMode == DIVA('AuditoryCortex','nFormants'))
0022       temp = DIVA('VocalTract','motor',grid_motor);
0023       xMode=temp{1};      
0024     elseif(nMode == DIVA('Sensory','nSensory'))
0025       temp = DIVA('VocalTract','motor',grid_motor);
0026       xMode = DIVA('Sensory','configuration',temp{2});
0027     end
0028     
0029     J1 = (xMode-grid_mode)./InverseParameters.DELTA;
0030     
0031     grid_motor = current_motor*ones(1,nMotor)-...
0032         InverseParameters.DELTA*diag(ones(1,nMotor));
0033     grid_mode  = current_mode(:,ones(1,nMotor));
0034     
0035     if(nMode == DIVA('AuditoryCortex','nFormants'))
0036       temp = DIVA('VocalTract','motor',grid_motor);
0037       xMode=temp{1};      
0038     elseif(nMode == DIVA('Sensory','nSensory'))
0039       temp = DIVA('VocalTract','motor',grid_motor);
0040       xMode = DIVA('Sensory','configuration',temp{2});
0041     end    
0042 
0043     J2 = (xMode-grid_mode)./-InverseParameters.DELTA;
0044     
0045     J = 0.5*(J1+J2);
0046     
0047     % Determine the inverse of the Jacobian
0048     %   Various tricks can be played here related to regularization
0049     %   This step is crucial for the inverse transformation and various
0050     %   different techniques can be applied to create a better solution
0051     %   At present the straightforward pseudoinverse is applied.
0052     
0053     J(find(isnan(J(:)))) = 0;
0054     [U,S,V]=svd(J,0);
0055     
0056     Jinv = pinv(J'*J+InverseParameters.LAMBDA*...
0057                 eye(size(J,2)))*J';
0058     
0059     Jnull=V(:,nMode+1:nMotor); 
0060     Jnull=Jnull*pinv(Jnull'*Jnull)*Jnull';    
0061     
0062 %     out=InverseParameters.ALPHA*...
0063 %         Jinv*mode_error - InverseParams.BETA*...
0064 %         Jnull*repmat(motor(:,ninput),1,N);
0065 
0066 out(:,ninput) = InverseParameters.ALPHA*Jinv*mode_error(:,ninput) -...
0067     InverseParameters.BETA*Jnull*motor(:,ninput);
0068 
0069 %     for ntarget = 1:N,
0070 %       % Jacobian formula
0071 %       out(:,ntarget) = InverseParameters.ALPHA*...
0072 %           Jinv*mode_error(:,ntarget) - InverseParameters.BETA*...
0073 %           Jnull*motor(:,ninput); % was Jnull*motor -> presumes one input/output
0074 %      end
0075   end

Generated on Tue 27-Mar-2007 12:06:24 by m2html © 2003