Just for completeness, I confirmed that this worked, and you know for sure what you’re getting in terms of the fit:
using ApproxFun, LsqFit, Plots
function model(x,p)
f = Fun(Chebyshev(Interval(0,2*pi)),p)
f.(x)
end
xvals = collect(0:.1:(2*pi))
yvals = sin.(xvals)
yvalsnoise = yvals .+ randn(length(xvals)) * 0.1
thefit = curve_fit(model,xvals,yvalsnoise,zeros(5))
f = Fun(Chebyshev(Interval(0,2*pi)),thefit.param)
plot(xvals,yvalsnoise)
plot!(xvals,f.(xvals))
You might find FastChebInterp easier for this. It can do Chebyshev regression from an arbitrary set of points to a given degree. (If you can evaluate at Chebyshev points it’s even better.). It’s lower level than ApproxFun, but that might be a plus here.
Of course, you have to think carefully about whether polynomials are a good fit here (pun intended), or if some specialized basis would be better for your problem, e.g. to capture the exponential tail. And of course you could also use polynomials after some change of variables or rescaling, or use more specialized polynomials like Gauss-Laguerre.
You could fit an LTI model directly using frequency-domain data, here are a number of methods implemented. I’d try subspaceid first. It sometimes requires some tweaking of the internal model order r to yield a good result. If you have time-domain data, you can of course use any other method from CSI as well.
Background information: I am currently looking at wind spectra. So there is no
simple dynamical system behind it, only mother earth…
They are usually characterized by the mean and the turbulence intensity. I want
to add one or more parameters to make it easier to compare different wind time
series in the frequency domain.
You can find some more info about this topic for example here:
There is never a simple linear system behind anything in practice, yet we use linear models to great effect all over the place 
The spectrum you show does not look very complicated, in fact, it looks like the frequency response of a low order linear system, so I’d imagine that white noise filtered through a linear system could produce a very similar spectrum. Whether or not that would be useful is another question, maybe not, but if you had a parametric model of a linear system that matches your data, you’d have an easy way to get the desired bandwidths.
This is giving indeed a much more streightforward approximation:
using PythonPlot, FastChebInterp
f_ex = 0.2
N = 10
Ts = 0.5
tmax = (N-1) * Ts
t = 0:Ts:tmax
signal = sin.(2π * f_ex .* t) # sin (2π f t)
f = chebregression(t, signal, 5)
plot(t, signal, label = "Signal")
t1=0.0:0.01:4.5
plot(t1, f.(t1), label = "Approximation")
legend()
grid(true)
Thanks for the hint!
I successfully applied Chebychev interpolation on the spectrum, and it works, but I needed a pretty high order, in this case 18:
So probably the suggestion from @baggepinnen is the better approach which would need less parameters…
You may want to discard the frequencies below 0.001 , then a polynom of 3th or 4th order could be sufficient, at least for your purposes. But sure, it is better to have some idea of a physical model to fit.
I added root finding with Roots.jl and now I also get the -3db and -10db frequencies calculated:
Main weakness of this approach is that I do not get error margins for the result…
If I need that, I might try one of the many other approaches suggested in this thread.
Thank you all for your expertise!
It would probably be possible to augment FastChebInterp to give error bars from chebregression, since it’s just a least-square fit in a particular basis.
It’d be interesting to compare this approach to Welch, if you have the time and the inclination…
Agreed, or at least below .0001 that flat region isn’t telling you anything much but to make a polynomial stay flat for a long distance would require more terms.
Which brings us back to Turing…
If I find the time, I will try… I mean, at some point I have to write a paper, and Welch has the advantage that I could cite a reference…
I would try a smoothing spline. There are a number of packages for this in Julia. I tend to using SmoothingSplines.jl.
For more ideas a number of people posted here.
You seem to have a solution that works, but I thought I would mention Savitsky-Golay as an option. I have found it to work very well.
https://www.juliabloggers.com/savitzky-golay-filters-julia/
edit: Use it to smooth out the noise before you do a fit.
Interesting! Thanks for sharing!