function y = fmod(x, cf)
    // Frequency modulation of an arbitrary signal
    // 
    // Usage:
    //       y = fmod(x, cf)
    //
    // where
    //       x:  signal vector (carrier)
    //       cf: vector containing the instantaneous
    //           frequency expansion factor 
    //       y:  frequency modulated otput vector
    //
    // This function works in such a way that any constant 
    // frequency component in the signal is frequency
    // modulated by an instantaneous factor cf(k), i.e.,
    // if f is the frequency of such component, after the
    // function has been applied the instantaneous 
    // frequency is cf(k)*f. In other words, the signal
    // is time-warped so to have such modulation. In order 
    // that y has the same duration as x, cf must average 1.
    // If cf is consistently < 1 then only part of the
    // original signal will be output (up to a duration
    // equal to the original one). If cf is consistently
    // > 1 then the whole signal will be output but the
    // duration will be smaller than the original one.
    //  
    // Technical description
    //
    // In order to accomplish the frequency modulation
    // the phase for a given frequency f is computed as 
    // the numerical version of the integral of the 
    // angular frequency:
    //
    //                    k
    // phi(k) = 2*pi*f * SUM(cf(n))/Fs
    //                   n=1
    //
    // If this is then divided by 2*pi*f/Fs we get the
    // phase expresed in terms of the sample index scale:
    //
    //       k
    // km = SUM(cf(n))
    //      n=1
    //
    // As this is independent of frequency, it can be
    // used for the whole signal. This means that, in
    // principle, the modulated signal could be obtained
    // by retrieving the values of the original signal
    // at these values of the index. However, km might 
    // be non-integer. For instance,  if cf were
    // constantly equal to 0.5 (which would correspond 
    // to lowering the whole signal one octave) then
    // kd = 0.5*k. This means that the original signal
    // should be sampled once in-between the original
    // samples. Since these values are unknown, it is
    // necessary to interpolate. 
    //
    // NOTE: It is possible to specify cf such that it
    // has some negative  values. This means an inversion 
    // of the time axis, which can be interpreted as a
    // locally retrograde reading of x. However, in the
    // event that the cumulative sum of cf becomes < 1,
    // we are back at the beginning of x. As there is no
    // information to extrapolate x towards indexes < 1,
    // the corresponding values of y will be assigned
    // NaN's  (not-a-number's). To prevent this error,
    // in general avoid negative values in cf. For the
    // same reason, cf(1) must be >= 1. 
    // 
    // -------------------------------
    // Author: Federico Miyara
    // Date:   2010-06-05
    //         2021-01-25

    // Compute phase in sample index scale 
    km = cumsum(cf);

    // Ensure that no attempt will be made to interpolate
    // further than the original duration of the signal.
    // To that end, the last index for which km does not 
    // go beyond the end of x is detected and km is then
    // truncated
    kmax = sum(km<length(x)); 
    km = km(1:kmax);

    // Interpolate
    y = interp1(1:length(x), x, km, "spline"); 

endfunction
