Building a spectral envelop of FFT'd audio I'm trying to figure out how to create a good spectral envelope of a signal.

Basically I am, at present, taking a windowed section of audio applying an FFT and generating a bar chart style representation of the magnitude of the spectrum.  

However to add to this I'd really like to have a good spectral envelope overlaid over the top.  Unfortunately I just can't seem to come up with a good solution.

I tried a simple peak picking algorithm where I find every peak in the spectrum (I literally just look for points where the current value is greater than the values in the buckets either side).  This returns me a set of peaks and I then use a catmull-rom interpolation to draw the line over the peaks.

This follows the spectrum reasonably as you can see in the following screenshot:



However its not ideal.  In that image for example you can see how it doesn't follow the spectrum on the lefthand side.

When you zoom out:



You see that the envelope is pretty jaggy.  Its not terrible but is this really the best way to build a spectral envelope?

Is there a better way of picking the ideal peaks so that it truly envelopes the spectrum?  If so can anyone point me towards algorithms?

Also is catmull-rom interpolation the best way of doing it?  Can anyone suggest a better interpolation method and show me how to implement it?

The question is language agnostic but I'm writing in C++.

Edit: Ok I found an implementation of Auto-Regressive modelling.

I then implemented the method as follows:

    std::vector&lt; float &gt; coefficients( 3 );
    AutoRegression( &amp;mSpecBuffer.front(), kFFTSizeDiv2, 3, &amp;coefficients.front(), MAXENTROPY );

    mPeakBuffer.clear();

    mPeakBuffer.push_back( Peak( 0, mSpecBuffer[0] ) );
    mPeakBuffer.push_back( Peak( 1, mSpecBuffer[1] ) );
    mPeakBuffer.push_back( Peak( 2, mSpecBuffer[2] ) );

    int x       = 3;
    int xMax    = kFFTSizeDiv2;
    while( x &lt; xMax )
    {
        const float k3  = (mPeakBuffer.end() - 3)-&gt;peakHeight;
        const float k2  = (mPeakBuffer.end() - 2)-&gt;peakHeight;
        const float k1  = (mPeakBuffer.end() - 1)-&gt;peakHeight;

        const float k   = (k1 * coefficients[0]) + (k2 * coefficients[1]) + (k3 * coefficients[2]);

        mPeakBuffer.push_back( Peak( x, k ) );
        x++;
    }


This gives me the following result:



Which to my eye looks a HELL of a lot better.

Zoomed right out it still looks pretty damned good:



So I just want to check, if this is correct?  If so I'll work along this path further :)
Two techniques commonly used for recovering the spectral envelope of an audio signal:


Low-order AR modeling. Estimate your AR coefficients, then plot the magnitude of the frequency response of the estimated all-pole filter.
Cepstral analysis. Take the logarithm of the magnitude of the FFT output. Take the FFT, truncate to keep only the first coefficients, take the inverse FFT, take the exponential. This is equivalent to applying a low-pass filter on the spectrum (with a log-scale magnitude).