function [ fct_fd std_dev_vector y2cMap ] = create_fd_objects( t, y, lambda_fct, lambda_stderr, degree )
    %
    % This function creates a functional data object to approximate the raw
    % data. Puts knots at each x values.
    %
    % Inputs parameters
    % -----------------
    % t    vector of t values to fit with a FDA object
    % y    vector of y values to fit with a FDA object
    % lambda_fct  smoothing parameter for the function object
    % lambda_stderr   smoothing parameter for the standard error object
    % degree - In which derivative are we looking for maxima and minima. To
    %          find max in the function, use degree = 0, and to find spurts (that is,
    %          max in velocity), use degree = 1
    %
    % Output values
    % -------------
    % fct_fd          functional object representing a smooth fit of the input data
    % std_dev_vector  vector of standard error measures at the x values 
    % y2cMap          FDA matrix mapping the data to the curve coefficients  

    rng = [t(1), t(size(t,1))];  % range of t values covered in the raw data

    %  -----------  Create fd objects   ----------------------------
    %  We penalize the highest derivative we will be using + 2.
    %  If we seek spurts and plateaus (i.e., degree = 1), we are interested
    %  in zero crossings in the jerk (3rd derivative).
    %  If we seek maxima in the function itself (degree = 0), we look at 
    %  acceleration zero crossings (2nd derivative)
    %
    penalty_degree = 4 + degree;
    norder         = penalty_degree + 2;
    nbasis         = length(t) + penalty_degree;
    fd_basis       = create_bspline_basis(rng, nbasis, norder, t);

    %  set up roughness penalty smoothing function smooth_basis

    Lfd            = int2Lfd(penalty_degree);
    fct_fdPar      = fdPar(fd_basis, Lfd, lambda_fct);

    %  smooth the data
    fct_fd         = smooth_basis(t, y, fct_fdPar);
    
    %  ----------------------------------------------------------
    %  Estimate standard error of measurement for velocity and 
    %    acceleration, re-smooth using the reciprocal of variance
    %    of estimate as a weight, and display results
    %  ----------------------------------------------------------

    %  set up function smooth_pos

    norderse       = 3 + degree;
    t_count = size(t, 1);
    nbasisse       = t_count + norderse - 2;
    stderrbasis    = create_bspline_basis(rng, nbasisse, norderse, t);
    Wfd0           = fd(zeros(nbasisse,1),stderrbasis);  %  initial value for Wfd
    %wtvec          = ones(t_count,1);  %  weight vector

    fd_fit         = eval_fd(t, fct_fd);
    fd_res         = y - fd_fit;   %  residuals
    fd_resmnsqr    = mean(fd_res.^2,2);  %  mean squared residuals

    %  positively smooth the mean squared residuals

    Lfdobj         = 1;  %  smooth toward a constant
    fct_fdPar      = fdPar(Wfd0, Lfdobj, lambda_stderr);

    Wfd            = smooth_pos(t, fd_resmnsqr, fct_fdPar);

    %  compute the variance and standard error of measurements
    fd_var         = eval_pos(t, Wfd);
    std_dev_vector = sqrt(fd_var);

    %  update weight vector for smoothing data
    wtvec          = 1./fd_var;
    wtvec          = wtvec./mean(wtvec);

    %  set up new smooth of the data using this weight vector
    Lfd            = int2Lfd(4);
    fct_fdpar      = fdPar(fd_basis, Lfd, lambda_stderr);

    %  smooth the data 
    [fct_fd, df, gcv, coef, SSE, penmat, y2cMap] = ...
        smooth_basis(t, y, fct_fdpar, wtvec);

return;