# Julia vs GNU Octave for plot fitting / finding peaks

**URL:** <https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237>\
**Category:** Signal and Image Processing\
**Tags:** data, curve-fitting\
**Created:** [June 29, 2020, 8:44am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237 "2020-06-29T08:44:57Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![doronbehar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/doronbehar/32/16004_2.png) [@doronbehar](https://discourse.julialang.org/u/doronbehar)\
**Post date:** [June 29, 2020, 8:44am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/1 "2020-06-29T08:44:57Z")

</div>

I’m currently using GNU octave to fit a simple cosine with a phase to many data sets stored in text files. For most data sets, the fit works well, but for some others, it doesn’t. I’m wondering whether I’d be able to get better results with Julia. I’ll explain:

The fit function I’m using, as I guess almost any other curve fit function, requires an initial guess for the parameters, as in [LsqFit.jl](https://github.com/JuliaNLSolvers/LsqFit.jl#existing-functionality). In GNU octave, in order to guess the frequency good enough, I rely upon [`findpeaks`](https://octave.sourceforge.io/signal/function/findpeaks.html) to calculate the average distance between the maximum peaks. Here’s an example of 2 measurements which `findpeaks` hasn’t found successfully the peaks in one of them, but it did so good enough for the other:

 ![bad-peaks](https://global.discourse-cdn.com/julialang/original/3X/5/8/58bdff87a8c2f1e10d765b28144048fbd91a8213.png)

There are a bit more details into the way I used `findpeaks`, but that’s irrelevant because this is a Julia Forum and not an octave forum.

Anyway, I’m no Julia user yet but I’ve heard good things about it, so I’m seeking out to hear about the experience of other users from the community, doing similar tasks to this. To be specific, `findpeaks` is my main disappointment out of GNU octave, and I wonder if somewhere in Julia’s community there’s a better `findpeaks` function.

---

<div class="post-metadata">

**Author:** ![MatFi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matfi/32/10002_2.png) [@MatFi](https://discourse.julialang.org/u/MatFi)\
**Post date:** [June 29, 2020, 9:17am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/2 "2020-06-29T09:17:30Z")

</div>

> [@doronbehar](#):
>
> I’m currently using GNU octave to fit a simple cosine with a phase to many data sets stored in text files. For most data sets, the fit works well, but for some others, it doesn’t. I’m wondering whether I’d be able to get better results with Julia. I’ll explain:

I would say that you would get the best / most reliable / performing results in this case through [FFT](https://en.wikipedia.org/wiki/Fast_Fourier_transform).

---

<div class="post-metadata">

**Author:** ![oheil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oheil/32/220745_2.png) [@oheil](https://discourse.julialang.org/u/oheil)\
**Post date:** [June 29, 2020, 10:21am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/3 "2020-06-29T10:21:37Z")

</div>

Hi @doronbehar and welcome.

To specifiy what has been said:

```julia
using Pkg; Pkg.add("FFTW"); Pkg.add("Plots")

```

```julia
using FFTW
using Plots

# Number of points
N=100
# Start time , end time
t0 = 0
tmax=4pi
# Sample period
Ts = (tmax-t0)/N
# time coordinate
t = t0:Ts:tmax

# signal
signal = sin.(t)

# Fourier Transform
F = fft(signal) |> fftshift
freqs = fftfreq(length(t), 1.0/Ts) |> fftshift

#plot
time_domain = scatter(t, signal, title = "Signal");
freq_domain = plot(freqs, abs.(F), title = "Spectrum", xlim=(-1, +1));
plot(time_domain, freq_domain, layout = 2)

#frequencies
f=freqs[findall(abs.(F) .> 0.1)]
#periods
p=1.0 ./ freqs[findall(abs.(F) .> 0.1)] 

```

It is adapted from here:  
[https://stackoverflow.com/questions/56030394/how-to-visualize-fft-of-a-signal-in-julia](https://stackoverflow.com/questions/56030394/how-to-visualize-fft-of-a-signal-in-julia)  
and works well down to N=25 samples for above time scale.

---

<div class="post-metadata">

**Author:** ![doronbehar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/doronbehar/32/16004_2.png) [@doronbehar](https://discourse.julialang.org/u/doronbehar)\
**Post date:** [June 29, 2020, 10:50am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/4 "2020-06-29T10:50:35Z")

</div>

Thank you @oheil for the code example and the welcome. But please note that I’m not looking for an answer to “how to do this”. I’ve also managed to write an equivalent code in Octave, and get a good approximation (also measured by hand) for the frequency but my fits are still disappointing.

What I mean is, that I’m interested in your _experiences_ with Julia in the wild - where you - Julia users had to fit numerous data sets to the same function, with more parameters then in your ideal code example, perhaps for a function such as:

![](https://global.discourse-cdn.com/julialang/original/3X/d/a/da005d5969f40f9de31449f901a4199b1f454885.gif)

Where A\_0 is an offset that also must be discovered by the curve fit.

---

<div class="post-metadata">

**Author:** ![oheil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oheil/32/220745_2.png) [@oheil](https://discourse.julialang.org/u/oheil)\
**Post date:** [June 29, 2020, 10:52am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/5 "2020-06-29T10:52:31Z")

</div>

In this case I missed the point 😁 No real world experience here.

---

<div class="post-metadata">

**Author:** ![platawiec](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/platawiec/32/31914_2.png) [@platawiec](https://discourse.julialang.org/u/platawiec)\
**Post date:** [June 29, 2020, 1:58pm UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/6 "2020-06-29T13:58:33Z")

</div>

Apologies if this isn’t quite what you’re looking for, but if what you are trying to do is fit data specifically to sinusoids, then you might want to use a routine like [harminv](https://github.com/NanoComp/harminv). Unfortunately this is a C library and there is no pure-Julia version, though if you know a bit of C you should be able to easily make a wrapper.

To this specific issue of `findpeaks`, it seems like there are three options, following [this conversation](https://github.com/JuliaDSP/DSP.jl/pull/369):

- [This fork of DSP.jl](https://github.com/tungli/DSP.jl), due to be merged soon as [this PR](https://github.com/JuliaDSP/DSP.jl/pull/369)
- [Peaks.jl](https://github.com/halleysfifthinc/Peaks.jl)
- [`Images.jl` `findlocalmaxima`](https://juliaimages.org/stable/function_reference/#Images.findlocalmaxima)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 29, 2020, 1:59pm UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/7 "2020-06-29T13:59:10Z")

</div>

> [@MatFi](#):
>
> I would say that you would get the best / most reliable / performing results in this case through [FFT](https://en.wikipedia.org/wiki/Fast_Fourier_transform).

You could use an FFT to get an initial guess, but I would then pass it into a least-square fitting routine (or similar) to refine the guess.

If you know _a priori_ that your data consists of a single cosine (plus noise?), then you can do better (possibly _far better_) than a Fourier transform (whose resolution is limited by the sampling time).

There are also other methods, e.g. based on variants of Prony’s method, to extract amplitude/frequency/phase from signals that are known _a priori_ to consist of a sequence of (possibly decaying) sinusoids (e.g. [GitHub - NanoComp/harminv: harmonic inversion algorithm of Mandelshtam: decompose signal into sum of decaying sinusoids](https://github.com/NanoComp/harminv) … it would be nice to re-implement this algorithm in Julia).

---

<div class="post-metadata">

**Author:** ![doronbehar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/doronbehar/32/16004_2.png) [@doronbehar](https://discourse.julialang.org/u/doronbehar)\
**Post date:** [June 29, 2020, 3:04pm UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/8 "2020-06-29T15:04:57Z")

</div>

Thank you all for the links and suggestions. I guess I’ll have to find out by myself how better is Julia vs GNU octave for this task. I’ll probably have to try out all of the options you all linked iteratively for all of my data sets, and do the comparison by my self.

If anyone’s interested, here’s an example data file (more [here](https://gitlab.com/doronbehar/physics1m-lab-harmonic-oscillator/-/tree/c38a2148cda9cec541540d650518c0bb8e5d2144/DRIVEN-DAMPED/measurements)):

[https://gist.github.com/f465318757df65010d26d85c30fd5f3e](https://gist.github.com/f465318757df65010d26d85c30fd5f3e)

And to show case how awful is Octave’s `findpeaks`, note how there are circles around the “peaks” - at almost every point:

 ![bad-peaks](https://global.discourse-cdn.com/julialang/original/3X/8/f/8fbe6e343e9f441759e26887b65dbc891f57093d.png)

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [June 29, 2020, 4:49pm UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/9 "2020-06-29T16:49:19Z")

</div>

> [@stevengj](#):
>
> If you know _a priori_ that your data consists of a single cosine (plus noise?), then you can do better (possibly _far better_) than a Fourier transform (whose resolution is limited by the sampling time).

In case of **a single tone** it is not limited as you can always pad by zeros.  
Padding by zeros is almost a perfect interpolation of the case of a single Sine.  
It won’t work, resolution wise, in case of multiple Sines / Cosines. As in this case you’re limited, in case of a bad SIR case, by the side lobes / width of main lobe which are governed by the time interval (The ratio between them isn’t, it is a function of the _Window_).

Regarding computational efficiency it is indeed sometimes better to have initial guess on a coarse DFT grid and then get better accuracy using Non Linear Least Squares methods.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [June 29, 2020, 6:05pm UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/10 "2020-06-29T18:05:05Z")

</div>

it looks to me like your data is actually constant for a while, and then steps to a new value (there’s a flat blue line involved). This induces infinite curvature at each step location, and makes it look like every location is a peak. If this is just an artifact of the plot, nevermind. But if this is because you have a finite resolution like a 4 bit or even 8 bit A/D converter then you might benefit by converting to floating point representation, and then doing a convolution with a smoothing kernel to remove this infinite curvature. For example

```julia
using DSP

newdata = conv(olddata,[sin(pi*i/6) for i from 0:6])
```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 30, 2020, 8:22am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/11 "2020-06-30T08:22:34Z")

</div>

> [@doronbehar](#):
>
> I guess I’ll have to find out by myself how better is Julia vs GNU octave for this task.

FWIW, I think most of this is about picking the right algorithm, which you can do equally well in Julia and Octave

There is nothing inherently magical in Julia that will make it better for this task, aside from a general preference for programming in Julia (of applicable).

---

<div class="post-metadata">

**Author:** ![doronbehar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/doronbehar/32/16004_2.png) [@doronbehar](https://discourse.julialang.org/u/doronbehar)\
**Post date:** [June 30, 2020, 11:47am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/12 "2020-06-30T11:47:16Z")

</div>

> [@dlakelan](#):
>
> This induces infinite curvature at each step location, and makes it look like every location is a peak.

That’s an excellent observation, you have also correctly guessed that my sensor is of “low” resolution, which is correct. I read about convolution on Wikipedia and I didn’t get much how is that related to “smoothing” my data. However, thinking about it I realised smoothing is exactly what I need to do before I start looking for peaks, and I’ve found:

- `https://www.mathworks.com/help/matlab/data_analysis/convolution-filter-to-smooth-data.html` Which uses a certain convolution function but I don’t know how to apply it to my data.
- [Function Reference: regdatasmooth](https://octave.sourceforge.io/data-smoothing/function/regdatasmooth.html) Which seem to make the life of `findpeaks` easier, but it takes dreadfully a long time.

> [@Tamas\_Papp](#):
>
> FWIW, I think most of this is about picking the right algorithm, which you can do equally well in Julia and Octave. There is nothing inherently magical in Julia that will make it better for this task, aside from a general preference for programming in Julia (of applicable).

Yea I guess you are right, I was mostly disappointed by Octave’s findpeaks, but reading @dlakelan’s comment, I guess I can’t blame it. Smoothing my data with [`regdatasmooth`](https://octave.sourceforge.io/data-smoothing/function/regdatasmooth.html) before using `findpeaks` gives better results, but still not perfect for all of my data sets - false peaks of the smoothed data are still detected:

 ![false-smooth-peaks](https://global.discourse-cdn.com/julialang/original/3X/6/0/60fa3886b8ba8a06e17b0d8979b71f47040e9bb0.png)

Writing a better algorithm that will make the call to `findpeaks` perfect, is the heart of my generic fit procedure. And it should probably take a long time, which I don’t currently have, though I do sense there’s some hope now, given the new knowledge I gained thanks to you all.

One thing that’s still missing for me when exploring Julia from above, is a _dedicated_ equivalent to this `smooth` function which is available on Matlab / Octave:

- `https://www.mathworks.com/help/matlab/ref/smoothdata.html`
- `https://octave.sourceforge.io/data-smoothing/function/regdatasmooth.html`

I understand that @dlakelan’s suggestion using `conv` might work, but I was wondering what could have helped me guess that 2nd argument to `conv`?

```julia
[sin(pi*i/6) for i from 0:6]

```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 30, 2020, 12:18pm UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/13 "2020-06-30T12:18:14Z")

</div>

> [@doronbehar](#):
>
> _dedicated_ equivalent to this `smooth` function

If you are looking for Tikhonov regularization, I think that

> **[GitHub - JuliaImageRecon/RegularizedLeastSquares.jl](https://github.com/JuliaImageRecon/RegularizedLeastSquares.jl)**
>
> Contribute to JuliaImageRecon/RegularizedLeastSquares.jl development by creating an account on GitHub.

has it implemented (and some other packages). And there are of course a ton of other smoothing algorithms implemented in Julia packages.

Incidentally, Matlab and Octave seem to be using different algorithms here. I am not sure that picking an arbitrary (if widely used) algorithm and calling it “smoothing” is a productive approach, there are zillion ways to “smooth” data, depending on the use case.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [June 30, 2020, 1:24pm UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/14 "2020-06-30T13:24:32Z")

</div>

> [@doronbehar](#):
>
> I understand that @dlakelan’s suggestion using `conv` might work, but I was wondering what could have helped me guess that 2nd argument to `conv` ?
> 
> ```julia
> [sin(pi*i/6) for i from 0:6]
> 
> ```

Basic knowledge of signal processing would have helped. Convolution with any non-negative kernel has smoothing qualities. There are tons of them; some have fancy names, some are claimed to be optimal under certain conditions. If you just choose one which is highest in the middle, symmetric, and strictly decreasing towards the edges, it will likely be good enough for your needs. Usually you would want to tune its width to your data though.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [June 30, 2020, 3:02pm UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/15 "2020-06-30T15:02:14Z")

</div>

> [@doronbehar](#):
>
> what could have helped me guess that 2nd argument to `conv` ?

A convolution replaces each data point with a weighted average of its neighboring points. Richard Hammings book “Digital Filters” is an excellent introduction, and very inexpensive.

I just chose that as an example, but it’s called the Hann kernel and does reasonably well at hammering down high frequency noise, which round-off error is. Yes you might need to widen it. Hammings book will be invaluable.

@doronbehar the convolution of the derivative of a kernel is the same as the negative of the derivative of the convolution… So you could do:

N=10 ## changing this allows you to try out different widths  
conv(data, [-cos(pi\*i/N) for i from 0:N])

and then find where the resulting function changes sign, and it should get you a long way towards finding the peaks. You’ll see that the -cos() function forms a half cycle going from -1 to 1… so the convolution is a weighted sum of differences between function values at different separations in time. (f(x+a) - f(x-a)) kind of thing. Since you’re just looking for where the function crosses 0 because that’s where the “smooth derivative” goes to zero and hence a peak is… this should be enough.

---

<div class="post-metadata">

**Author:** ![doronbehar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/doronbehar/32/16004_2.png) [@doronbehar](https://discourse.julialang.org/u/doronbehar)\
**Post date:** [July 1, 2020, 7:25am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/16 "2020-07-01T07:25:21Z")

</div>

Thanks again to you all for the links and suggestions. I conclude by now, that I’ll need to learn signal processing thoroughly at some point.

One last question: Does Julia return a vector of the same size when you use `conv`? On Octave, it doesn’t - if I use as a 2nd argument what I suspect is the matlab equivalent of:

```julia
[sin(pi*i/6) for i from 0:6]

```

I.e `sin((0:6)*pi/6)`, I get a vector with a size larger by 6, and then I don’t know how to plot it or do anything else with it, because the time vector is of the same size as before.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 1, 2020, 7:33am UTC](https://discourse.julialang.org/t/julia-vs-gnu-octave-for-plot-fitting-finding-peaks/42237/17 "2020-07-01T07:33:34Z")

</div>

You can easily find the frequency using, e.g., the ESPRIT method

```julia
using DSP
y = sin.(0:0.1:10) .+ 0.1 .* randn.()
esprit(y, 10, 2)

julia> esprit(y, 10, 2)
2-element Array{Float64,1}:
 -0.015982377052346328
  0.015982377052346328

julia> 0.1/(2pi)
0.015915494309189534

```

After that, it’s a simple least-squares problem to find the coefficients in the model

k\_s \sin(\omega t) + k\_c \cos(\omega t)