# Logarithmic Smoothing - A Tutorial and Two Package Ideas

**URL:** https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877
**Category:** Signal and Image Processing
**Tags:** package, tutorials, smoothing
**Created:** [July 27, 2022, 3:22pm UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877 "2022-07-27T15:22:12Z")
**Posts on this page:** 9
**Page:** 1

<div class="post-metadata">

### Author: ![LostInMusic](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lostinmusic/32/7486_2.png) [@LostInMusic](https://discourse.julialang.org/u/LostInMusic)
#### Post date: [July 27, 2022, 3:22pm UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/1 "2022-07-27T15:22:12Z")

</div>

As a person interested in acoustics, I often plot FFT results on double-logarithmic scales. The results are often pretty noisy and smoothing is typically applied, for example with a width of 1/12 octaves. Unfortunately, so far there has been no easy way to perform smoothing and resampling on a logarithmic axis in Julia.

I’ve now written a little demo/tutorial notebook to change that: [GitHub - Firionus/logarithmic\_smoothing\_julia\_demo: How to perform logarithmic smoothing in Julia](https://github.com/Firionus/logarithmic_smoothing_julia_demo)

If you also need log smoothing, I hope you like the way it works in the tutorial. I’ve certainly implemented it badly way too many times before and now finally found the time to do it properly. If you see any mistakes or have questions, please leave a comment 😃

As for the broader appeal, the tutorial uses two unregistered packages I wrote:

- [GitHub - Firionus/NonlinearSequences.jl: logarithmically and power spaced values](https://github.com/Firionus/NonlinearSequences.jl) - basically a nice logspace API
- [GitHub - Firionus/NonuniformResampling1D.jl: Resamples uniformly sampled 1D data at nonuniform locations](https://github.com/Firionus/NonuniformResampling1D.jl) - a simple window-based resampler that allows arbitrary output coordinates and adjusts the smoothing width based on the distance between output coordinates

Before registering these packages I’d like to take in some feedback. Is their functionality something you can see yourself using? Do you think they should be registered? Do you see any obvious problems with their code? If so, let me know in the comments. I’d love to hear from the community ❤

---

<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 27, 2022, 5:02pm UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/2 "2022-07-27T17:02:44Z")

</div>

Have you considered plotting something like the [Welch periodigram](https://docs.juliadsp.org/stable/periodograms/#DSP.Periodograms.welch_pgram) rather than the basic FFT (standard periodigram)?

---

<div class="post-metadata">

### Author: ![LostInMusic](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lostinmusic/32/7486_2.png) [@LostInMusic](https://discourse.julialang.org/u/LostInMusic)
#### Post date: [July 27, 2022, 7:59pm UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/3 "2022-07-27T19:59:44Z")

</div>

The Welch method also smooths on the linear frequency axis like most other algorithms, so it isn’t applicable to this problem. Maybe this quick plot serves to visualize that:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/4/44451f7c6645cf488fd7a4271fa72f7aca587de8.jpeg)

With the Welch method, either the high frequencies are still noisy (Welch 50sp) or the resolution at low frequencies is way too low (Welch 250sp). The log-smoothed version has similar detail throughout the frequency range when plotted against a logarithmic frequency scale.  
Additionally, in acoustics most of the time you’re looking at the FFT of an impulse response, so energy concentration in the time domain is quite high. But the Welch method seems designed more for noise-like rather than click-like signals.

Since you brought up a method that windows in the time domain: Instead of performing a single FFT and then smoothing, you can perform a DFT where you apply a different window for each frequency and the width is proportional to the inverse of the frequency. Variants of this exist in many acoustic measurements softwares, e.g. “frequency dependent window” in REW and Acourate, “Adaptive Window” in MLSSA ([http://www.mlssa.com/pdf/DRA-MLSSA-Manual-10WI-8.pdf](http://www.mlssa.com/pdf/DRA-MLSSA-Manual-10WI-8.pdf)) or “Dual Gate” in ARTA.  
The result is also log-smoothed, but reduces the influence of parts outside the center of the impulse response for high frequencies. This is thought to better resemble human hearing, and there are good reasons to believe that. But I haven’t seen any proper empiric research on that yet. (Paper on that topic: [https://www.researchgate.net/publication/241628057\_Generalized\_Fractional-Octave\_Smoothing\_of\_Audio\_and\_Acoustic\_Responses](https://www.researchgate.net/publication/241628057_Generalized_Fractional-Octave_Smoothing_of_Audio_and_Acoustic_Responses)) But enough of me rambling for tonight 🙈, thanks for commenting 👍

---

<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 28, 2022, 7:49am UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/4 "2022-07-28T07:49:41Z")

</div>

Cool, thanks for sharing 🙂

> [@LostInMusic](#):
>
> you can perform a DFT where you apply a different window for each frequency and the width is proportional to the inverse of the frequency.

this reminds me of the [Wigner-Ville](https://en.wikipedia.org/wiki/Wigner_distribution_function) spectrogram.

---

<div class="post-metadata">

### Author: ![JFViola](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jfviola/32/203766_2.png) [@JFViola](https://discourse.julialang.org/u/JFViola)
#### Post date: [November 10, 2023, 4:02pm UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/5 "2023-11-10T16:02:37Z")

</div>

Hi - a bit late to the party, but i wanted to encourage you to further pursue this project! Logsmoothing is also quite important in other disciplines like astronomy & climate science (estimating climate variability from PSDs), so if there are no packages yet, they would be very helpful. I am new/ not really using Julia yet, but esp. with climate science there is the movement towards Julia, so it might be fitting.  
Cheers.

---

<div class="post-metadata">

### Author: ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [November 10, 2023, 10:21pm UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/6 "2023-11-10T22:21:56Z")

</div>

FWIW, a simple smoothing alternative using `Loess.jl` is provided below that seems to produce results similar to those of the [package example](https://github.com/Firionus/logarithmic_smoothing_julia_demo/blob/main/example_picture.png).

> **CODE**
>
> ```julia
> using Loess, DelimitedFiles, Plots
> 
> # https://github.com/Firionus/logarithmic_smoothing_julia_demo/blob/main/example_frequency_response.csv
> 
> file = raw"C:\Users\jrafa\OneDrive\Julia_Code\DATASETS\example_frequency_response.csv"
> 
> f, a, p = eachcol(readdlm(file, ',', skipstart=2)) # skip first row with 0 frequency
> 
> N = length(f)
> n = min(2000, N) # number of spectra points to retain
> e0 = exp(log(N)/n)
> ix = collect(round(Int, e0^i) for i in 0:n)
> 
> model = loess(f[ix], a[ix], span=0.016)
> la = Loess.predict(model, f[ix])
> 
> plot(f, a, xaxis=:log, xlims=(20, 20e3), ylims=(40, 120), c=:blues)
> plot!(f[ix], la, c=:blue)
> 
> ```

Result:

 ![Loess_log_spectrum](https://global.discourse-cdn.com/julialang/original/3X/6/7/67e13faeb7760dfb253f51ca7e797447be8a12e8.png)

---

<div class="post-metadata">

### Author: ![LostInMusic](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lostinmusic/32/7486_2.png) [@LostInMusic](https://discourse.julialang.org/u/LostInMusic)
#### Post date: [November 15, 2023, 2:50pm UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/7 "2023-11-15T14:50:06Z")

</div>

Hey, that’s a nice demo of how to approximate it 😊, and almost certainly with better performance than my package.

But notice the higher variability in your plot compared to my example around 4 - 10 kHz? That’s because you throw out data randomly, which again creates noise that isn’t there when you consider all data points.  
That’s why I wrote my own package for this, and unfortunately I don’t see a way to use loess.jl with unevenly spaced data points. I don’t think there’s a way to “weigh” the data points.

But thanks for sharing nonetheless, I appreciate it 👍

---

<div class="post-metadata">

### Author: ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [November 15, 2023, 4:21pm UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/8 "2023-11-15T16:21:44Z")

</div>

> [@LostInMusic](#):
>
> But notice the higher variability in your plot compared to my example around 4 - 10 kHz?

You can turn the `span` keyword knob in `loess()`. For example, try `span=0.02` in the code provided above.

---

<div class="post-metadata">

### Author: ![LostInMusic](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lostinmusic/32/7486_2.png) [@LostInMusic](https://discourse.julialang.org/u/LostInMusic)
#### Post date: [November 16, 2023, 3:00am UTC](https://discourse.julialang.org/t/logarithmic-smoothing-a-tutorial-and-two-package-ideas/84877/9 "2023-11-16T03:00:29Z")

</div>

I ran your code on my own computer, and it started looking even more similar to my example! Then I upped the number of retained points to N, and wow, that result is similar:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/1/5129e4fe369754bfbbec6696fcffcc2fee6a59b7.png)

> **CODE**
>
> ```plaintext
> using Loess, DelimitedFiles, Plots
> 
> # https://github.com/Firionus/logarithmic_smoothing_julia_demo/blob/main/example_frequency_response.csv
> 
> file = raw"example_frequency_response.csv"
> 
> f, a, p = eachcol(readdlm(file, ',', skipstart=2)) # skip first row with 0 frequency
> 
> N = length(f)
> n = N#min(2000, N) # number of spectra points to retain
> e0 = exp(log(N)/n)
> ix = collect(round(Int, e0^i) for i in 0:n)
> 
> print(length(ix))
> 
> model = loess(f[ix], a[ix], span=0.017)
> la = Loess.predict(model, f[ix])
> 
> plot(f, a, xaxis=:log, xlims=(20, 20e3), ylims=(40, 120), c=:blues, alpha=1, label="raw data")
> plot!(f[ix], la, c=:blue, linestyle=:dot, label="loess w/ repeat hack")
> 
> using NonuniformResampling1D, NonlinearSequences
> 
> resolution = 1/12 # octaves
> oversampling = 16 # results in `oversampling/resolution` many points per octave
> flog = octspace(20, 20e3, resolution/oversampling)
> fax = range(f[1], f[end], length=length(f))
> ylog = nuresample(fax, a, flog, hann_window(1*oversampling))
> 
> plot!(flog, ylog, c=:orange, linestyle=:dot, label="nuresample")
> 
> ```

That’s interesting! Your variant actually uses repetition of frequency/amplitude values at lower frequencies and leaves out values at higher frequencies:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/e/2eca33c956b6f92c2dc5d2e62a44b827bb141747.png)

The result is nearest-neighbor behavior at low frequencies, where NonunformResampling1D automatically uses upsampling to maintain smooth shapes:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/1/c1d76cc3a96132e962d98489dc4eb7a18ee4a1c3.png)

Your approach is a bit hacky, IMO, but the results clearly are very similar to my approach, with much less code 👍
