function [y, ymean] = bands2spectdens(Lp, fo, N, faxis)
    // Convert band spectrum to spectral density
    // Usage:
    //       y = bands2spectdens(Lp, fo, N, faxis)
    // where
    //       Lp:    Vector of band levels in dB
    //       fo:    Vector of centre frequencies of bands 
    //              in Hz
    //       N:     Number of frequencies where the spectral
    //              density will be given. If N is omitted, 
    //              only the bounds of the intervals where
    //              the spectral density is linear are given  
    //       faxis: Type of frequency distribution ("log",
    //              "lin")              
    //       y:     Two-row matrix whose first row contains 
    //              frequencies in Hz and whose second row 
    //              contains the corresponding spectral 
    //              densities in Pa^2/Hz 
    //       ymaen: Vector of mean values of the spectral
    //              density over each band
    //
    // NOTES: 1) The spectral density is given as a continuous
    //           piecewise linear function with the same total
    //           energy in each band as in the given bands.
    //        2) The mean spectral density over adjacent bands
    //           must be different 
    //
    // --------------------------------
    // Author: Federico Miyara
    // Date:   2021-03-09     

    if 1==2
        // Example for testing purposes only --DELETE--
        Lp = [20   40 50  60  65  62   55   58   50   45    35];
        fo = [16 31.5 63 125 250 500 1000 2000 4000 8000 16000];
        N = 100;
        faxis = "log";
    end
    
    // Number of bands
    n = length(Lp); 

    // Lower and upper limits of bands
    fs(1:n-1) = sqrt(fo(1:n-1).*fo(2:n));
    fs(n) = fo(n)*fo(n)/fs(n-1);
    fi(2:n) = fs(1:n-1);
    fi(1) = fo(1)*fo(1)/fi(2);
    
    // Average spectral density in each band
    Pref = 20e-6;
    p2 = Pref^2*10.^(Lp/10)./(fs - fi);
    
    if argn(1)==2 then
        ymean = p2;
    end

    // Vector of values of spectral density at the
    // boundaries between bands
    q = zeros(1, n+1);
    
    // First band
    if p2(1)<p2(2)
        q(1) = p2(1)/2;
    elseif p2(1)>p2(2)
        q(1) = p2(1) + (p2(1) - p2(2))/2;
    end
    
    // Inner bands
    //q(2:n) = (p2(1:n-1) + p2(2:n))/2; 
    q(2:n) = sqrt((p2(1:n-1) .* p2(2:n))); 
        
    // Last band
    if p2(n-1)>p2(n)
        q(n+1) = p2(n)/2;
    elseif p2(n-1)<p2(n)
        q(n+1) = p2(n) + (p2(n) - p2(n-1))/2;
    end

    // Interpolating points (r,u)
    r = [];
    u = []; 
       
    // First band
    if p2(1)>p2(2)
        // Density is maximum at the first band 
        ro = (fi(1) + fs(1))/2;
        r = [r, fi(1), ro];
        u = [u, q(1), p2(1)]
    elseif p2(1)<p2(2)
        // Density is minimum at the first band         
        ro = fs(1)*(q(2)-p2(1)); 
        ro = ro + fi(1)*(p2(1)-q(1));
        ro = ro / (q(2)-q(1));
        r = [r, fi(1), ro   ];
        u = [u, q(1),  p2(1)];
    end

    // Inner bands
    for k=2:n-1
        // Monotone increasing
        if p2(k-1)<p2(k) & p2(k)<p2(k+1)
            ro = fs(k)*(q(k+1) - p2(k)) + fi(k)*(p2(k) - q(k));
            ro = ro / (q(k+1) - q(k));
            r = [r, fi(k), ro];
            u = [u, q(k), p2(k)];
        end
        // Monotone decreasing
        if p2(k-1)>p2(k) & p2(k)>p2(k+1)
            ro = fs(k)*(q(k+1) - p2(k)) + fi(k)*(p2(k) - q(k));
            ro = ro / (q(k+1) - q(k));
            r = [r, fi(k), ro];
            u = [u, q(k), p2(k)];
        end
        // Maximum
        if p2(k-1)<p2(k) & p2(k)>p2(k+1)
            ro = (p2(k) - q(k)) / (2*p2(k) - q(k) - q(k+1));
            ro = fi(k) + ro*(fs(k) - fi(k));
            r = [r, fi(k), ro];
            M(k) = (p2(k) - q(k)) * (q(k+1) - q(k));
            M(k) = M(k) / (2*p2(k) - q(k) - q(k+1));     
            M(k) =  2*p2(k) - q(k+1) + M(k);
            u = [u, q(k), M(k)];        
        end
        // Minimum
        if p2(k-1)>p2(k) & p2(k)<p2(k+1)
            B = fs(k)*q(k+1)/q(k) - fi(k);
            B = B - 2*p2(k)*(fs(k) - fi(k))/q(k);
            C = fi(k)*fs(k)*q(k+1)/q(k);    
            ro1 = -B + sqrt(B^2 - 4*C);        
            ro2 = -B - sqrt(B^2 - 4*C);
            r = [r, fi(k), ro1,     ro2];
            u = [u, q(k),  p2(k)/2, p2(k)/2];
        end
    end
    
    // Last band    
    if p2(n)>p2(n-1)
        // Density is maximum at the last band 
        ro = (fi(n) + fs(n))/2;
        r = [r, fi(n), ro,    fs(n)];
        u = [u, q(n),  p2(n), q(n+1)]
    elseif p2(n)<p2(n-1)
        // Density is minimum at the last band         
        ro = fs(n) * (q(n+1) - p2(n)); 
        ro = ro + fi(n) * (p2(n) - q(n));
        ro = ro / (q(n+1) - q(n));
        r = [r, fi(n), ro,    fs(n) ];
        u = [u, q(n),  p2(n), q(n+1)];
    end
    
    // Decide number of frequencies from arguments
    // and compute the interpolated spectral density
    if argn(2)<3
        // No number of frequencies is provided. Use
        // the minimum set of frequencies defining
        // the piecewise linear approximation
        y = [r; u];
    else
        // The number of frequencies is provided 
        if argn(2)==4
            // The type of frequencye distribution is 
            // given  
            if faxis=="log" 
                // If it is explicitly logarithmic, 
                // create a logarithmivally istributed 
                // frequency set
                ff = logspace(log10(fi(1)), log10(fs(n)), N);
            else
                // In any other case create a linearly
                // distributed frequency set 
                ff = linspace(fi(1), fs(n), N);
            end
        elseif argn(2)==3
            // If there is no explicit mention of how
            // should be the frequencies distributed,
            // distribute them linearly  
            ff = linspace(fi(1), fs(n), N);
        end
        // Interpolate linearly
        // NOTE: since logspace doesn't provide an exact
        //       match of the first and last points,
        //       the extrapolation feature of interp1 
        //       must be used to prevent NaN's.
        y = [ff; interp1(r, u, ff, "linear", "extrap")];
    end
    
 
    if 1==2
        // Piecewise constant interpolation for 
        // comparison
        rr = [];
        uu = [];
        for k=1:n
            rr = [rr, fi(k), fs(k)];
            uu = [uu, p2(k), p2(k)];
        end
        
        figure(1);
        gcf().background = 8;
        clf()
        plot("ln", r, real(u))    
        plot("ln", rr, real(uu),"r") 

        figure(2);
        gcf().background = 8;
        clf()
        plot("ln",r,10*log10(real(u)))    
        plot("ln", rr, 10*log10(real(uu)),"r")
    end
               

    
endfunction
    
