aliasing issue with floating-point resampling ratio

[ericphanson]

When calling resample with a floating-point value, it takes a different codepath than with a rational value. The floating-point version seems more susceptible to aliasing issues.

Here is an example, minimized from a real-world issue, that demonstrates the problem:

using DSP, CairoMakie

function hpf(x, freq; fs)
    return filtfilt(digitalfilter(Highpass(freq; fs), Butterworth(4)), x)
end

function lpf(x, freq; fs)
    return filtfilt(digitalfilter(Lowpass(freq; fs), Butterworth(4)), x)
end

function middle_third(x)
    third = div(length(x), 3)
    return x[third:(2 * third + 1)]
end

function plt_filters(x; fs)
    x = lpf(x, 35.0; fs)
    x = hpf(x, 0.5; fs)
    return x
end

sine_wave(freq_hz) = sin.(2π .* freq_hz .* ts)

n = 45 # seconds
fs = 250
ts = range(0, n; length=fs * n)
data = 0.05 * sine_wave(0.75) + 0.01 * sine_wave(5.0) + 0.025 * sine_wave(10.0) +
       # lots of high frequency noise
       sum(100 * sine_wave(f) for f in 90:125)

resampling_ratio = 1 / 1.00592
resampled = resample(data, resampling_ratio)
ts_resampled = resample(ts, resampling_ratio)

rational_resampled = resample(data, rationalize(resampling_ratio))
rational_ts_resampled = resample(ts, rationalize(resampling_ratio))

# individual axes
fig = let
    fig = Figure()
    colors = Makie.wong_colors()
    ax_kwargs = (; ylabel="Data", limits=(nothing, nothing, -1, 1))

    ax1 = Axis(fig[1, 1]; ax_kwargs...)
    l1 = lines!(ax1, middle_third(ts), middle_third(plt_filters(data; fs)); color=colors[1])

    ax2 = Axis(fig[2, 1]; ax_kwargs...)
    l2 = lines!(ax2, middle_third(ts_resampled), middle_third(plt_filters(resampled; fs));
                color=colors[2])

    ax3 = Axis(fig[3, 1]; xlabel="Time (s)", ax_kwargs...)
    l3 = lines!(ax3, middle_third(rational_ts_resampled),
                middle_third(plt_filters(rational_resampled; fs)); color=colors[3])

    Legend(fig[4, 1], [l1, l2, l3], ["original", "resampled", "rational resampled"];
           orientation=:horizontal)
    fig
end

save("mwe.png", fig)

# combined
fig = let
    fig = Figure()
    colors = Makie.wong_colors()
    ax_kwargs = (; ylabel="Data", limits=(nothing, nothing, -1, 1))
    ax = Axis(fig[1, 1]; xlabel="Time (s)", ax_kwargs...)

    lines!(ax, middle_third(ts_resampled), middle_third(plt_filters(resampled; fs));
           color=colors[2], label="resampled")
    lines!(ax, middle_third(rational_ts_resampled),
           middle_third(plt_filters(rational_resampled; fs)); color=colors[3],
           label="rational resampled")
    lines!(ax, middle_third(ts), middle_third(plt_filters(data; fs)); color=colors[1],
           label="original")

    axislegend(ax)
    fig
end
save("mwe-combined.png", fig)
fig

[wheeheee]

Could you test this with 0.7.9 to see if this might be a regression with the changes in 0.8.0?

[ericphanson]

This actually is on 0.7.9, I will try to run it on 0.8.0 but I need to update my code for changes (no fs argument to Lowpass)

[wheeheee]

Sorry, didn't quite notice. Fingers crossed it magically disappears in 0.8...

[ericphanson]

Unfortunately not, same results on 0.8

DSP v0.8 code

using DSP, CairoMakie

function hpf(x, freq; fs)
    return filtfilt(digitalfilter(Highpass(freq), Butterworth(4); fs), x)
end

function lpf(x, freq; fs)
    return filtfilt(digitalfilter(Lowpass(freq), Butterworth(4); fs), x)
end

function middle_third(x)
    third = div(length(x), 3)
    return x[third:(2 * third + 1)]
end

function plt_filters(x; fs)
    x = lpf(x, 35.0; fs)
    x = hpf(x, 0.5; fs)
    return x
end

sine_wave(freq_hz) = sin.(2π .* freq_hz .* ts)

n = 45 # seconds
fs = 250
ts = range(0, n; length=fs * n)
data = 0.05 * sine_wave(0.75) + 0.01 * sine_wave(5.0) + 0.025 * sine_wave(10.0) +
       # lots of high frequency noise
       sum(100 * sine_wave(f) for f in 90:125)

resampling_ratio = 1 / 1.00592
resampled = resample(data, resampling_ratio)
ts_resampled = resample(ts, resampling_ratio)

rational_resampled = resample(data, rationalize(resampling_ratio))
rational_ts_resampled = resample(ts, rationalize(resampling_ratio))

# individual axes
fig = let
    fig = Figure()
    colors = Makie.wong_colors()
    ax_kwargs = (; ylabel="Data", limits=(nothing, nothing, -1, 1))

    ax1 = Axis(fig[1, 1]; ax_kwargs...)
    l1 = lines!(ax1, middle_third(ts), middle_third(plt_filters(data; fs)); color=colors[1])

    ax2 = Axis(fig[2, 1]; ax_kwargs...)
    l2 = lines!(ax2, middle_third(ts_resampled), middle_third(plt_filters(resampled; fs));
                color=colors[2])

    ax3 = Axis(fig[3, 1]; xlabel="Time (s)", ax_kwargs...)
    l3 = lines!(ax3, middle_third(rational_ts_resampled),
                middle_third(plt_filters(rational_resampled; fs)); color=colors[3])

    Legend(fig[4, 1], [l1, l2, l3], ["original", "resampled", "rational resampled"];
           orientation=:horizontal)
    fig
end

save("mwe.png", fig)

# combined
fig = let
    fig = Figure()
    colors = Makie.wong_colors()
    ax_kwargs = (; ylabel="Data", limits=(nothing, nothing, -1, 1))
    ax = Axis(fig[1, 1]; xlabel="Time (s)", ax_kwargs...)

    lines!(ax, middle_third(ts_resampled), middle_third(plt_filters(resampled; fs));
           color=colors[2], label="resampled")
    lines!(ax, middle_third(rational_ts_resampled),
           middle_third(plt_filters(rational_resampled; fs)); color=colors[3],
           label="rational resampled")
    lines!(ax, middle_third(ts), middle_third(plt_filters(data; fs)); color=colors[1],
           label="original")

    axislegend(ax)
    fig
end
save("mwe-combined.png", fig)
fig

[wheeheee]

I tried to resample using a FIRFilter and different Nϕ, and found that it gets smoother with Nϕ=64 (there are still artifacts, presumably when ϕAccumulator jumps back). But then the resampled signal also becomes shifted / scaled, and after a bit of digging I found that these lines here don't seem to pass on Nϕ correctly...

DSP.jl/src/Filters/stream_filt.jl

Lines 198 to 201 in 06f8776

	function FIRFilter(rate::AbstractFloat, Nϕ::Integer=32)
	h = resample_filter(rate, Nϕ)
	FIRFilter(h, rate)
	end

On L200, FIRFilter looks like it should have an additional argument of Nϕ, but upon fixing it, the smoothing effect previously observed disappears, so that's really not the problem here.

[wheeheee]

Ok, saw wrongly. Indeed, choosing a higher value of Nϕ does reduce artifacts. Currently, that parameter is not exposed in resample, and also as raised in another issue resample_filter isn't really documented either, probably should do that too.

There are still going to be artifacts, but less, if you use something like

resample_phases(s, rate, Nϕ=32) = filt(FIRFilter(resample_filter(rate, Nϕ), rate, Nϕ), s)
resampled = resample_phases(data, resampling_ratio, 128)
ts_resampled = resample_phases(ts, resampling_ratio, 128)

[wheeheee]

You could also play around with rel_bw in resample_filter to see if it can give you what you want. I find setting Nϕ=64 and rel_bw=0.7 looks good for this example at least.

[martinholters]

When working on #596, I din't take a closer look at the algorithm, but from what I recollect, it does upsampling by an integer (Nϕ?) followed by linear interpolation to obtain the proper points "in-between". Neither of those operations should cause those peaks significantly exceeding the input signal's range. So I do think this exposes some kind of bug, which may be circumvented by choosing suitable parameters, but preferably should be fixed, of course.

[wheeheee]

Just to add on, the Nϕ used to create pfb for the rational resampling example here is quite large (6250 vs default 32 for float). So this might not just be a problem with FIRArbitrary resampling...
edit: the graphs shown here are the result of filtering the resampled signals with a bandpass filter...without that, the resampled signals don't seem to exceed the input signal's range.

[ericphanson]

the graphs shown here are the result of filtering the resampled signals with a bandpass filter...without that, the resampled signals don't seem to exceed the input signal's range.

Yeah, fwiw, that's why I was describing it as an aliasing issue; the original signal has a lot of high frequency noise, but it is filtered out cleanly by the bandpass. The resampled version has these unexpected spikes when filtered the same way.

It seems this arbitrary sampling rate resampling is a tricky business. Searching, I came across https://www.mathworks.com/help/dsp/ug/efficient-sample-rate-conversion-between-arbitrary-factors.html and https://www.mathworks.com/help/dsp/ref/designrateconverter.html which seem to be MATLAB's solution here. I don't have a MATLAB license or I'd check how it does on this signal. But maybe a multi-stage approach like that would do better here.

[martinholters]

Ah, right, I hadn't looked at how those plots are generated carefully enough. So the (unfiltered) data contains those spikes at 1 s intervals due to the "high frequency noise", which really is a series of sines 1 Hz apart (90 Hz to 125 Hz). Then the highpass filter should eliminate those, but for the (irrational-)resampled signal, a significant amount remains due to aliasing components below (or just slightly above) the filter's cutoff. It should be noted, though, that the spikes are still attenuated by some 80 dB or so, which isn't great, but not catastrophic either.

A way to more directly visualize the aliasing effects is by looking at the spectrogram of a linear sine sweep. Consider a long sweep from 0 to pi:

sweep = let N=2^24
	[sin(0.5*pi*n^2/N) for n in 1:N]
end

As expected, it's spectrogram consists of a single diagonal and some windowing artifacts:

OTOH, the spectrogram for resample(sweep, resampling_ratio) (with resampling_ratio as above) shows quite some aliasing:

(Note that the colorbar is focused on the background; the original signal is above 30dB, so the SNR is not as bad as one might think at first glance.)

For comparison, rational resampling (resample(sweep, rationalize(resampling_ratio))) does better:

For completeness, I've created the spectrograms with:

function sgramplot(x)
	sgram = spectrogram(x, 4096; window=hanning)
	
	fig = Figure()
	ax = Axis(fig[1,1]; xlabel="ω", ylabel="n")
	hm = heatmap!(ax, freq(sgram), time(sgram), pow2db.(power(sgram)); colorrange=(-130, -30))

	Colorbar(fig[:, end+1], hm)
	return fig
end

So contrary to my first impression, I no longer think this issue exposes a bug. Rather, the choice of algorithm and default parameters leads to a performance/quality trade-off that delivers insufficient quality in this case. So what @wheeheee has proposed might indeed be the solution here, not merely a work-around.

The main take-way here is probably that we should improve the documentation on how to achieve this. From a look at the Mathworks documents, I get the feeling that it might also be worthwhile to replace the linear interpolation we're doing with cubic interpolation, but that would require a deep dive into how FIRArbitrary works, I guess. In the long term, something like designRateConverter would surely be nice, but for now, I guess we have to leave the higher-level design work to the user.

[wheeheee]

I tried out designRateConverter in MATLAB, although I'm not familiar with MATLAB code, and got this out (which doesn't look like the resampled + bandpass-filtered signals above?). It seems that designRateConverter did not use a multi-stage approach, and I didn't know how to coax it into doing so. Still, for what it's worth:

Details

function sw = sine_wave(freq_hz, ts)
    sw = sinpi(2 * freq_hz .* ts);
end

function xt = middle_third(x)
    third = fix(length(x) / 3);
    xt = x(third:(2 * third + 1));
end

n = 45
fs = 250
ts = linspace(0, n, fs * n)
rg = 90:125
lots = arrayfun(@(x) sine_wave(x, ts), rg, UniformOutput=false)
data = 0.05 * sine_wave(0.75, ts) + 0.01 * sine_wave(5.0, ts) ...
     + 0.025 * sine_wave(10.0, ts) + ...
       100 * sum(cat(1, lots{:}))

data = data'

SRC = designRateConverter(InputSampleRate=fs, ...
    OutputSampleRate=fs/1.00592, ...
    Verbose=true, StopbandAttenuation=60)

resampled_data = SRC(data)
resampled_ts = linspace(0, n, length(resampled_data))

plot(middle_third(resampled_ts), middle_third(resampled_data))

filted = bandpass(resampled_data, [0.5 35.0], fs)

plot(middle_third(resampled_ts), middle_third(filted))

The last plot (resampled, bandpass-filtered) looks like this