It sounds like you want to show the range and average (mean?) of the spectrum.
One way to do this is to use Welch's Method. Essentially, you divide the data in the time domain into sections (the length of the section will depend on the frequency resolution you need / want and how much data you have). In practice you would normally use overlapping sections (50% overlap). You calculate the spectrum estimate of each of these sections and average across all of the estimates (average). To obtain the range, take the min and max of the estimates.
Apologies if this is not the case, but often a poor estimator is used when calculating the spectrum estimate (periodograms tend to have high variance for example). There are quite a few ways to calculate the estimates but one of the more optimal ways (in many senses of the word) is to use the multitaper method. There is also a package on cran written by a colleague of mine and former student of the developer of this technique .
The method tapers the data using orthogonal tapers with the weighted average being the final estimate. Due to the averaging, you normally end up with a less variable estimate.
Without needing to tweak a lot,
specEst <- spec.mtm(yourData) #yourData is a time series or vector
will estimate the spectrum with default parameters and plot that estimate.
Regardless, my suggestion would be to use the Welch method to obtain mean and range and use the Multitaper Method to estimate the spectrum of each section.