# How to approximate a noisy spectrum
**URL:** https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519
**Category:** General Usage
**Tags:** question
**Created:** [May 28, 2023, 4:45pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519 "2023-05-28T16:45:00Z")
**Posts on this page:** 16
**Page:** 2
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [May 28, 2023, 9:53pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/21 "2023-05-28T21:53:01Z")
Just for completeness, I confirmed that this worked, and you know for sure what you’re getting in terms of the fit:
```julia
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))
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 28, 2023, 10:14pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/22 "2023-05-28T22:14:23Z")
> [@ufechner7](#):
>
> But I have only a vector of frequency, amplitude pairs, and not a function. ApproxFun needs a function as input.
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.
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [May 29, 2023, 4:32am UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/23 "2023-05-29T04:32:34Z")
You could fit an LTI model directly using frequency-domain data, [here are a number of methods implemented](https://baggepinnen.github.io/ControlSystemIdentification.jl/dev/freq/). I’d try [`subspaceid`](https://baggepinnen.github.io/ControlSystemIdentification.jl/dev/freq/#Statespace) 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.
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [May 29, 2023, 7:56am UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/24 "2023-05-29T07:56:28Z")
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:
> **[Study of Dynamic Response Characteristics of the Wind Turbine Based on...](https://www.mdpi.com/2076-3417/9/12/2392)**
>
> The present research envisages a method for calculating the dynamic responses of the wind turbines under typhoon. The measured power spectrum and inverse Fourier transform are used to generate the fluctuating wind field in the eyewall of the typhoon....
> **[nwpsaf-kn-tr-008.pdf](https://nwp-saf.eumetsat.int/publications/tech_reports/nwpsaf-kn-tr-008.pdf)**
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [May 29, 2023, 8:06am UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/25 "2023-05-29T08:06:57Z")
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.
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [May 29, 2023, 10:43am UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/26 "2023-05-29T10:43:03Z")
> [@stevengj](#):
>
> You might find FastChebInterp easier for this.
This is giving indeed a much more streightforward approximation:
```julia
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!
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [May 29, 2023, 11:25am UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/27 "2023-05-29T11:25:10Z")
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…
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [May 29, 2023, 12:19pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/28 "2023-05-29T12:19:40Z")
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.
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [May 29, 2023, 12:36pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/29 "2023-05-29T12:36:09Z")
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!
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 29, 2023, 1:03pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/30 "2023-05-29T13:03:07Z")
> [@ufechner7](#):
>
> Main weakness of this approach is that I do not get error margins for the result…
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.
---
### Author:  [@mbaz](https://discourse.julialang.org/u/mbaz)
#### Post date: [May 29, 2023, 1:07pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/31 "2023-05-29T13:07:29Z")
It’d be interesting to compare this approach to Welch, if you have the time and the inclination…
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [May 29, 2023, 2:49pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/32 "2023-05-29T14:49:12Z")
> [@Eben60](#):
>
> You may want to discard the frequencies below 0.001 ,
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.
> [@ufechner7](#):
>
> Main weakness of this approach is that I do not get error margins for the result…
Which brings us back to Turing…
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [May 29, 2023, 3:00pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/33 "2023-05-29T15:00:22Z")
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…
---
### Author:  [@Jake](https://discourse.julialang.org/u/Jake)
#### Post date: [May 29, 2023, 6:12pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/34 "2023-05-29T18:12:37Z")
I would try a smoothing spline. There are a number of packages for this in Julia. I tend to using [SmoothingSplines.jl](https://juliapackages.com/p/smoothingsplines).
For more ideas a number of people posted [here](https://discourse.julialang.org/t/smoothing-interpolation-for-tachometer-data/56117).
---
### Author:  [@purplishrock](https://discourse.julialang.org/u/purplishrock)
#### Post date: [May 29, 2023, 6:25pm UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/35 "2023-05-29T18:25:56Z")
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/](https://www.juliabloggers.com/savitzky-golay-filters-julia/)
edit: Use it to smooth out the noise before you do a fit.
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [June 14, 2023, 11:48am UTC](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519/36 "2023-06-14T11:48:04Z")
Interesting! Thanks for sharing!
[Previous page](https://discourse.julialang.org/t/how-to-approximate-a-noisy-spectrum/99519.md?page=1)