DSP.jl power spectra normalization

Specific Domains Signal and Image Processing
1

Can anyone explain the rationale behind the FFT normalization done in DSP.jl, for instance in the periodogram function ?

Here below is the concerned code of the fft2pow! function, which computes the power spectra.
I am concerned on the one-sided spectrum. In the code, the power spectrum is normalized by m1 for the first frequency bin, by m2 for bins 2:n-1 and by m1 or m2 for the last bin depending on whether the size of the window is even or odd:

function fft2pow!(out::Array{T}, s_fft::Vector{Complex{T}}, nfft::Int, r::Real, onesided::Bool, offset::Int=0) where T
m1 = convert(T, 1/r)
n = length(s_fft)
if onesided
m2 = convert(T, 2/r)
out[offset+1] += abs2(s_fft[1])*m1
for i = 2:n-1
@inbounds out[offset+i] += abs2(s_fft[i])*m2
end
out[offset+n] += abs2(s_fft[end])*ifelse(iseven(nfft), m1, m2)
else
…

Looking at the call of fft2pow! in the periodogram function at line 275 of unit periodograms.jl it seems to me that r here above is passed as fs*norm2, where fs is the sampling rate and norm2 is sum(abs2, window), where window is the taper.

I don’t see the interest in taking into account the sampling rate in the normalization of the power spectra. This leads to undesirable behavior. For example, I computed the power spectra for a pure sine wave at an exact discrete Fourier frequency (using a rectangular window) two times, changing only the sampling rate and i obtained different powers at the only non-zero frequency bin.

3 Likes
2

Not sure why this two year old topic popped up at the top of my “latest” sorting so sorry if this is weird necromancy, but I don’t think that that normalization is too out there. If you changed the sampling rate of those two functions as you describe, wouldn’t that lead to two sequences that have different norms? In that sense, accounting for the sampling in the size of the periodogram coefficients (and thus the norm of the sequence of coefficients) seems like maintaining a proper invariance.

For what it’s worth, by the way, the Percival and Walden multitaper book, which I would personally consider a gospel, is very disciplined about this type of normalization. I don’t use DSP.jl and probably never will and I won’t pretend that I am careful about normalizing by the sampling rate in my applications where everything has the same sampling rate, but I could imagine for some applications that it’s important to keep track of.

1 Like
3

@Marco-Congedo, do you see any inconsistency with other codes?
Say with Matlab (some examples posted on the web).

DSP.jl’s periodogram seems to depend on the signal length considered but not on the sampling rate:

DSP.jl code used: 
using DSP, Plots, Printf

plot_font = "Computer Modern";
default(fontfamily=plot_font, framestyle=:box, fg_color_legend = nothing, tickfontsize=7, guidefontsize=8)


########    Test#1:  vary sampling period Δt    ########

# Input:
f0 = 20;                # sinousoid frequency [Hz]
Δt = [1, 2, 4]*1e-3;    # sampling period [s]
fs = 1.0 ./Δt           # sampling frequency [s]
T = 0.2                 # input signal duration [s]
t = [0:Δt:T-Δt for Δt in Δt]
s = [sin.(2Ď€*f0*t) for t in t]

# Compute:
S = [periodogram(s; fs=fs) for (fs,s) in zip(fs,s)]

# Plot:
str = permutedims("Δt = " .* string.(Δt*1e3) .* " ms")
p1 = plot(t,s, xlabel="Time [s]", ylabel="Amplitude", label=str)
p2 = plot(;xlabel="Freq [Hz]",ylabel="Periodogram (dB)", xlims=(0,125), ylims=(-400,40), grid=true)
[plot!(p2, S.freq, DSP.pow2db.(S.power), label=str[i]) for (i,S) in pairs(S)]
f0_ = S[1].freq[argmax(S[1].power)]    # OK = 20 Hz
S0_ = S[1].power[argmax(S[1].power)]   # depends on signal duration T
S0dB = 10*log10(S0_)
str = @sprintf("f0 = %.0f Hz,  S0 = %.0f dB",f0_, S0dB)
annotate!(f0_ + 2, S0_, text(str,:blue,:left,7, plot_font))
plot(p1,p2, legend=:outertop)


########    Test#2:  vary signal length T    ########

# Input:
f0 = 20;                # sinousoid frequency [Hz]
Δt = 2e-3;              # sampling period [s]
fs = 1.0 /Δt            # sampling frequency [s]
T = [1.0, 0.4, 0.2]     # input signal duration [s]
t = [0:Δt:T-Δt for T in T]
s = [sin.(2Ď€*f0*t) for t in t]

# Compute:
S = [periodogram(s; fs=fs) for s in s]

# Plot:
str = permutedims("T = " .* string.(T) .* " s")
p1 = plot(t,s, xlabel="Time [s]", ylabel="Amplitude", label=str)
p2 = plot(;xlabel="Freq [Hz]",ylabel="Periodogram (dB)", xlims=(0,125), ylims=(-400,40), grid=true)
[plot!(p2, S.freq, DSP.pow2db.(S.power), label=str[i]) for (i,S) in pairs(S)]
p3 = deepcopy(p2)
plot!(p3,xlims=(19,22), ylims=(-20,3), grid=true, title="zoom")
f0_ = [S.freq[argmax(S.power)] for S in S]    # OK = 20 Hz
S0_ = [S.power[argmax(S.power)] for S in S]    # depends on signal duration T
S0dB_ = 10*log10.(S0_)
str = [@sprintf("f0=%.0f Hz, S0= %.0f dB",f0, S0dB) for (f0,S0dB) in zip(f0_,S0dB_) ]
[annotate!(f0+0.2, S0dB, text(str,:blue,:left,5, plot_font)) for (f0,S0dB,str) in zip(f0_,S0dB_,str) ]
plot(p1,p2,p3, layout=(1,3),legend=:outertop)
1 Like
4

NB:
The above results look consistent with DSP.jl’s periodogram code comment:

The computed periodogram is normalized so that the area under the periodogram is equal to the uncentered variance (or average power) of the original signal

2 Likes
5

No @rafael.guerra, i didn’t notice any inconsistency with other established code, i was just wondering while writing my code.