using Gadfly
using ColorSchemes
using Statistics
using Printf
using FFTW
set_default_plot_size(20cm, 15cm)
style(background_color=colorant"white");

An astute reader may have noticed something odd when reading the first post on power measurement. The FFT size used was 10,000. Anyone who has learned about the FFT algorithm knows that algorithm is drastically faster when an FFT size that’s a power of 2 is used (eg 2^10 == 1024).

So why use an FFT size of 10,000? Well, there’s a source of loss that wasn’t accounted for in that post and under certain conditions this loss is 0 dB. This loss is called scalloping loss and scalloping loss is zero when the tone in the FFT (10 kHz in this case) falls on an FFT bin boundary.

When computing the frequency domain representation of a signal, the representation is broken up into bins. The sampling rate in the first post was 100 kHz and there were 10,000 bins meaning that each bin represents 10 Hz of bandwidth.

When the tone frequency falls on that 10 Hz bonudary (eg 40,500 Hz), there’s no scalloping loss. However, when the tone frequency doesn’t fall on a boundary (eg 40,509 Hz) the signal is affected.

Scalloping loss comes from the FFT window that’s being used. In this case no window is explicitly applied, which really means that a rectangualar window is being used (since all the samples before and after the signal are 0).

In this post, we’ll examine scalloping loss and how you can compensate for it. First, we’ll setup a signal that’s identical to the one used in the first post. Then, we’ll take 3 different FFTs of that signal and see the affect of scalloping. Next, we develop a method for compensating for scalloping loss. Finally, we plot scalloping loss directly and then the output of the signal with scalloping loss compensated for.

Step 1 - Signal Setup

# Signal setup
duration = 1                        # Signal durations [s]
fs = 100e3                          # Sampling rate, [Hz]
f = 10e3                            # Tone frequency [Hz]
t = 1/fs:1/fs:duration              # Signal time vector [s]
x = exp.(im * 2π * f * t)           # Signal with no noise
noise = randn(ComplexF64, size(x))  # Noise
xₙ = x .+ noise                     # xₙ, signal with noise 
;

Plot the first 100 samples of the signal

n = 100
p_1 = plot(Guide.title("x (complex valued signal)"), Guide.ylabel("Amplitude", orientation=:vertical), Guide.xlabel("Time [s]"))
p_2 = plot(Guide.title("xₙ (complex valued signal)"), Guide.ylabel("Amplitude",  orientation=:vertical), Guide.xlabel("Time [s]"))

push!(p_1, layer(x=t[1:n], y=real(x[1:n]), Geom.line, color=[colorant"blue"]))
push!(p_1, layer(x=t[1:n], y=imag(x[1:n]), Geom.line, color=[colorant"red"]))
push!(p_2, layer(x=t[1:n], y=real(xₙ[1:n]), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=t[1:n], y=imag(xₙ[1:n]), Geom.line, color=[colorant"red"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

Step 2 - Three different FFTs of the same signal

In this section we see three different FFTs of the same signal. The only difference is the FFT Size (10000, 8192 and 4096).

The top plot shows the FFT output over the whole bandwidth and the bottom plot shows a zoom-in of the signal peak.

fft_size = 10000
X = (abs.(fft(xₙ[1:fft_size]))).^2 ./ fs ./ fft_size / (fft_size/length(xₙ))
peak_power_10000 = findmax(10 .* log10.(X))[1]
freqs = FFTW.fftfreq(fft_size, fs)

p_1 = plot()
push!(p_1, Guide.title("Frequency Response w/ FFT Size = $fft_size"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_1, Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), ymin=-80, ymax=10))
push!(p_1, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))

p_2 = plot()
push!(p_2, Guide.title("Zoomed-In Frequency Response = $fft_size"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_2, Coord.cartesian(xmin=8, xmax=12, ymin=-5, ymax=5))
push!(p_2, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=[-50, 50], y=[0, 0], Geom.line, color=[colorant"red"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

fft_size = 8192
X = (abs.(fft(xₙ[1:fft_size]))).^2 ./ fs ./ fft_size / (fft_size/length(xₙ))
peak_power_8192 = findmax(10 .* log10.(X))[1];
freqs = FFTW.fftfreq(fft_size, fs)

p_1 = plot()
push!(p_1, Guide.title("Frequency Response w/ FFT Size = 8192"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_1, Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), ymin=-80, ymax=10))
push!(p_1, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))

p_2 = plot()
push!(p_2, Guide.title("Zoomed-In Frequency Response w/ FFT Size = 8192"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_2, Coord.cartesian(xmin=8, xmax=12, ymin=-10, ymax=10))
push!(p_2, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=[-50, 50], y=[0, 0], Geom.line, color=[colorant"red"]))


stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

fft_size = 4096
X = (abs.(fft(xₙ[1:fft_size]))).^2 ./ fs ./ fft_size / (fft_size/length(xₙ))
freqs = FFTW.fftfreq(fft_size, fs)

p_1 = plot()
push!(p_1, Guide.title("Frequency Response w/ FFT Size = $fft_size"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_1, Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), ymin=-80, ymax=10))
push!(p_1, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))

p_2 = plot()
push!(p_2, Guide.title("Zoomed-In Frequency Response w/ FFT Size = $fft_size"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_2, Coord.cartesian(xmin=8, xmax=12, ymin=-10, ymax=10))
push!(p_2, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=[-50, 50], y=[0, 0], Geom.line, color=[colorant"red"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

These results show that as we change the FFT Size the signal power drops from 0 dB, to -0.5 dB to -2.5 dB.

Here we can directly see the affects of scalloping. We’re only changing the FFT Size and the power of the signal is decreasing, that doesn’t make sense!

@printf "Peak Power for FFT Size of 10000 = %2.3f dB\n" peak_power_10000
@printf "Peak Power for FFT Size of 8192 = %2.3f dB\n" peak_power_8192
@printf "Peak Power for FFT Size of 4096 = %2.3f dB\n" peak_power_4096
;

Peak Power for FFT Size of 10000 = 0.041 dB
Peak Power for FFT Size of 8192 = -0.552 dB
Peak Power for FFT Size of 4096 = -2.543 dB

Step 3 - Compensate for Scalloping Loss

Here we write a function to compensate for scalloping loss. The frequency domain representation of the rectangular window in a sinc() function. We apply that sinc function at the correct frequency to get the scalloping loss that we need to compensate for.

Note that scalloping loss that we compensate for matches the loss in step 2 for each FFT Size.

function sinc(x::Float64)
    sin(x)/x;
end

function compute_scalloping_loss(Fs::Float64, fft_size::Int64, f::Float64)
    res_bw = Fs/fft_size
    a = π/(res_bw) * mod((f - (res_bw/2)), res_bw) - π/2
    if a == 0
        # Avoid divide by 0 in sinc(x)
        b = 0  
    else
        b = 10 * log10(sinc(a) .^ 2)
    end
    b
end
;

sloss_10000 = compute_scalloping_loss(fs, 10000, f)
sloss_8192 = compute_scalloping_loss(fs, 8192, f)
sloss_4096 = compute_scalloping_loss(fs, 4096, f)

@printf "Scalloping loss for FFT Size = 10000 is %2.3f [dB]\n"  sloss_10000
@printf "Scalloping loss for FFT Size = 8192 is %2.3f [dB]\n"  sloss_8192
@printf "Scalloping loss for FFT Size = 4096 is %2.3f [dB]\n"  sloss_4096

Scalloping loss for FFT Size = 10000 is 0.000 [dB]
Scalloping loss for FFT Size = 8192 is -0.579 [dB]
Scalloping loss for FFT Size = 4096 is -2.420 [dB]

Step 4 - Plot Scalloping Loss and Compensating For It

In this last step a tone is swept from 9980 Hz to 10020 Hz in 0.08 Hz increments. The power of the tone is recorded for each increment. The power of the tone is then plotted as a function of frequency.

l = 500
X_comps = zeros(l)
scalloping_loss = zeros(l)
N = 10000
freq = zeros(l)

for (ii, kk) in enumerate(range(-20, 20, length=l))
    F = f + kk                              # Tone frequency 
    x = exp.(im * 2π * F .* t)              # Signal, a tone at frequency F
    X = fft(x[1:N])                         # FFT of the signal, FFT Size is 10k
    X_comp = abs.(X).^2 / N / fs / (N/fs)   # Signal, compensated for FFT Size, sampling rate and FFT bin size
    
    s_loss = compute_scalloping_loss(fs, N, F)  # Scalloping loss
    scalloping_loss[ii] = s_loss
    
    X_comps[ii] = 10*log10(findmax(abs.(X_comp))[1]) - s_loss  # X, compensated for all of the above and scalloping loss
    freq[ii] = F
end

Here we plot the scalloping loss as a function of frequency. Notice that when the FFT size divides the tone frequency evenly (eg 9990 Hz, 10000 Hz, etc) there is no scalloping loss.

plot(x=freq, 
     y=scalloping_loss, 
     Geom.line, 
     Guide.xlabel("Frequency [Hz]"), 
     Guide.ylabel("Scalloping Loss [dB]"),
     Guide.title("Scalloping Loss as a function of Frequency for FFT Size = $N"))

This is also a good place to discuss wherer the name “scalloped” comes from. It turns out that curtains and other textiles with this shape are referred to as “scalloped”.

A scalloped curtain:

In the last plot of this post, the power of the tone is shown as the tone is swept from 9980 Hz to 10020 Hz. All parameters are compensated for in this plot:

FFT Size (10k)
Sampling rate (100 kHz)
FFT Bin Size
Scalloping loss

With all of these items factored in, the power of the tone stays at a constant 0 dB/Hz for all frequencies examined.

plot(x=freq, 
     y=X_norms,
     Geom.line,  
     Coord.cartesian(ymin=-1, ymax=1),
     Guide.xlabel("Frequency [Hz]"), 
     Guide.ylabel("Peak Power [dB/Hz]"),
     Guide.title("Peak Power of Swept Tone after Compensating for Scalloping Loss"))

Recent Posts

May 11, 2022 - Scalloping Loss in the FFT

Step 1 - Signal Setup

Step 2 - Three different FFTs of the same signal

Step 3 - Compensate for Scalloping Loss

Step 4 - Plot Scalloping Loss and Compensating For It