1
$\begingroup$

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

enter image description here

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

enter image description here

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.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Do you have an extremely low-off the chart outlier around frequency 1.7? $\endgroup$
    – gung - Reinstate Monica
    Sep 10, 2015 at 3:10
  • $\begingroup$ 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 $\endgroup$
    – user2113499
    Sep 10, 2015 at 3:45
  • $\begingroup$ Try min(spc). What do you get? $\endgroup$
    – gung - Reinstate Monica
    Sep 10, 2015 at 3:53
  • $\begingroup$ I get 1.4e-9 as the min $\endgroup$
    – user2113499
    Sep 10, 2015 at 3:55

1 Answer 1

Reset to default
0
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$
    – user2113499
    Sep 10, 2015 at 16:32
  • $\begingroup$ My pleasure! I hope it turns out to be useful. $\endgroup$
    – driegert
    Sep 10, 2015 at 19:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.