Logarithmic Smoothing - A Tutorial and Two Package Ideas

LostInMusic · July 27, 2022, 3:22pm

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

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 - basically a nice logspace API
GitHub - Firionus/NonuniformResampling1D.jl: Resamples uniformly sampled 1D data at nonuniform locations - 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

baggepinnen · July 27, 2022, 5:02pm

Have you considered plotting something like the Welch periodigram rather than the basic FFT (standard periodigram)?

LostInMusic · July 27, 2022, 7:59pm

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:

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) 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) But enough of me rambling for tonight , thanks for commenting

baggepinnen · July 28, 2022, 7:49am

Cool, thanks for sharing

this reminds me of the Wigner-Ville spectrogram.

JFViola · November 10, 2023, 4:02pm

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.

rafael.guerra · November 10, 2023, 10:21pm

FWIW, a simple smoothing alternative using Loess.jl is provided below that seems to produce results similar to those of the package example.

CODE

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:

LostInMusic · November 15, 2023, 2:50pm

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

rafael.guerra · November 15, 2023, 4:21pm

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

LostInMusic · November 16, 2023, 3:00am

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:

CODE

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:

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

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