# Adding specific trend line to spectra frequency graph in R

I have a plot of spectra vs frequency and I am trying to add a specific line through the data and what I have right now is

plot(freq, spc, log='xy', type='l')
y.loess <- loess(spc ~ freq, span=0.8, data.frame(x=freq, y=spc))
y.predict <- predict(y.loess, data.frame(x=freq)
lines(freq,y.predict, col='red')


This gives me the following

The black part of the graph is correct and what I need but the red line is incorrect what I need should look something like

I thought loess would work but it's not quite what I am going for. How do I add a line to my data to make it look like the second picture? I feel like nls will work somehow but I do not know how to add the formula for it.

• Do you have an extremely low-off the chart outlier around frequency 1.7? Sep 10, 2015 at 3:10
• I am not sure why the loess line is doing that. I plotted it with a a ylim from 1e-12 to 1e-3 and there were no other points lower than what is shown up above Sep 10, 2015 at 3:45
• Try min(spc). What do you get? Sep 10, 2015 at 3:53
• I get 1.4e-9 as the min Sep 10, 2015 at 3:55

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.

• Thank you for your suggestion. After looking more into Welch's method I believe that is exactly what I need. I didn't know about method. Thank you for introducing it! Sep 10, 2015 at 16:32
• My pleasure! I hope it turns out to be useful. Sep 10, 2015 at 19:32