Fuglevand.m

Go to the documentation of this file.

namespace models{
namespace motorunit{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class Fuglevand
  :public handle {
    public:

                fsamp = 1000;
len = 10;
                RR = 30;
        fmin = 8;
fmax = 35;
        sd_ISI = 0.1;
                RP = 100;
RT = 3;
        lengthWin = 1000;
        Resolution = 100;

        
    public:

        pende;

        
    private:

        r;

        fSeed;

        
    public: /* ( Dependent ) */

        Seed;

        
    public: /* ( Dependent ) */

        Fuglevand() {
            this.pende= (this.fmax-this.fmin)/60; /*  gain g_e in the paper */


            this.Seed= 1;
        }
        function forces = getForces(colvec<double> mu,T) {
            
            if nargin < 3
                T = 2000;
                if nargin < 2
                    mu = [.1 5];
                end
            end

            /*  Add window of single peak to simulation time */
            T = T + this.lengthWin;

            firing = zeros(1,T);
            forces = firing;

            k = 2;
            gain = 1;

            /*  recruitment threshold excitation */
            rte = exp(log(this.Resolution)*mu(1))*this.RR/this.Resolution;
            exDrive = mu(2)*this.RR/10;

            if rte <= exDrive*1.001

                /*  peak twitch force */
                P = exp(log(this.Resolution)*mu(1))*this.RP/this.Resolution;
                Tfac = 90*(1./P)^(log(this.RT)/log(this.RP));

                /*  instantaneous discharge frequency */
                finst = this.pende*(exDrive-rte)+this.fmin;
                finst = min(finst, this.fmax);

                /*  previous firing time
                 * rand is a pseudorandom value drawn from the standard uniform distribution on the open interval (0,1). */
                istprec = this.r.rand/finst;

                istatt = 0;  /*  current firing time */

                while(istatt+0.02 < T/1000)
                    /*  variability due to the jitter
                     * standard normal distri with mean 0 and standard deviation sd_ISI
                     * rand is a pseudorandom value drawn from the standard normal distribution */
                    var_jitter = this.sd_ISI*this.r.randn*0;

                    /*  moment of activation (time) */
                    istatt = istprec+(1+var_jitter)/finst;
                    /*  check that the firing is not too close to the previous one */
                    if(istatt < istprec+0.02)
                        istatt = istprec+0.02;
                    end
                    firing(round(istatt*1000)) = 1;
                    Ti_ISI = Tfac/(istatt-istprec)/1000;
                    if(Ti_ISI < 0.4)
                        gain(k) = 1.0;/* #ok */

                    else
                        gain(k) = (1-exp(-2*Ti_ISI^3))/Ti_ISI/((1-exp(-2*0.4^3))/0.4);/* #ok */

                    end
                    k=k+1;
                    istprec = istatt;
                end /*  while */

                hlp = [0:round(this.lengthWin/1000*this.fsamp)]/this.fsamp/(Tfac/1000);
                winCoeff = P.*(hlp.*exp(1-hlp));
                winCoeff = [0 winCoeff(winCoeff>1.e-4)];

                spikes = find(firing==1);
                for i=1:size(spikes,2)
                    if(spikes(i)+size(winCoeff,2) <= length(forces))
                        idx = spikes(i) + (1:size(winCoeff,2));
                        forces(idx) = forces(idx)+gain(i)*winCoeff;
                    end
                end
            end
            /*  Subtract window of single peak from force */
            forces = forces(1:end-this.lengthWin);
        }


        
#if 0 //mtoc++: 'set.Seed'
function  Seed(value) {
            this.r= RandStream(" mt19937ar "," Seed ",value);
            this.fSeed= value;
        }

#endif


        
#if 0 //mtoc++: 'get.Seed'
function value = Seed() {
            value = this.fSeed;
        }

#endif

    
    public: /* ( Static ) */ /* ( Dependent ) */

        static function  test_ForceResponse() {
            f = models.motorunit.Fuglevand;
            [x,y] = meshgrid(0:.05:1,0:.4:10);
            p = [x(:) y(:)]^t;
            force = zeros(1,numel(x));
            for k=1:size(p,2)
                fmu = f.getForces(p(:,k));
                force(k) = max(fmu);
            end
            force = reshape(force,size(x,1),[]);
            surf(x,y,force," FaceColor "," interp "," EdgeColor "," interp ");
        }

    
};
}
}