# DSP.jl power spectra normalization

**URL:** <https://discourse.julialang.org/t/dsp-jl-power-spectra-normalization/26162>\
**Category:** Signal and Image Processing\
**Tags:** dsp\
**Created:** [July 8, 2019, 9:41pm UTC](https://discourse.julialang.org/t/dsp-jl-power-spectra-normalization/26162 "2019-07-08T21:41:37Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marco-Congedo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco-congedo/32/7321_2.png) [@Marco-Congedo](https://discourse.julialang.org/u/Marco-Congedo)\
**Post date:** [July 8, 2019, 9:41pm UTC](https://discourse.julialang.org/t/dsp-jl-power-spectra-normalization/26162/1 "2019-07-08T21:41:37Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [July 31, 2021, 10:00pm UTC](https://discourse.julialang.org/t/dsp-jl-power-spectra-normalization/26162/2 "2021-07-31T22:00:16Z")

</div>

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.

---

<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:** [August 1, 2021, 10:48am UTC](https://discourse.julialang.org/t/dsp-jl-power-spectra-normalization/26162/3 "2021-08-01T10:48:25Z")

</div>

@Marco-Congedo, do you see any inconsistency with other codes?  
Say with Matlab (some examples [posted](https://ww2.mathworks.cn/help/signal/ug/measuring-the-power-of-deterministic-periodic-signals.html) on the web).

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

- Variable sampling rate (fixed signal window length):  

- Variable signal window length (fixed sampling rate):  

> **DSP.jl code used:**
>
> ```julia
> 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)
> 
> ```

---

<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:** [August 1, 2021, 1:02pm UTC](https://discourse.julialang.org/t/dsp-jl-power-spectra-normalization/26162/4 "2021-08-01T13:02:12Z")

</div>

**NB:**  
The above results look consistent with [DSP.jl’s periodogram code comment](https://github.com/JuliaDSP/DSP.jl/blob/ea75bfbfdf9f7181695b3024f2491567a4df395e/src/periodograms.jl):

_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_

---

<div class="post-metadata">

**Author:** ![Marco-Congedo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco-congedo/32/7321_2.png) [@Marco-Congedo](https://discourse.julialang.org/u/Marco-Congedo)\
**Post date:** [August 1, 2021, 2:48pm UTC](https://discourse.julialang.org/t/dsp-jl-power-spectra-normalization/26162/5 "2021-08-01T14:48:41Z")

</div>

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