# Comparing custom Welch periodogram with Julia's Welch periodogram

I wrote my own PSD struct in Julia for practice purposes. It involves two type of constructors: a "direct" PSD and a PSD based on Welch' method (the one that matters to us). I provide the code for the curious, but it is not essential to my question:

mutable struct PSD
freq::Vector
spectrum::Vector

# Direct PSD
function PSD(x::Vector, sampling_rate::Int)
N = length(x)
hann = hanning(N) # Hanning window
x = x .* hann
ft = abs2.(fft(x))
ft = ft[1:(div(N, 2) + 1)] # Make one sided
freq = [i for i in 0:( length(ft) - 1)] .* sampling_rate/N
normalization = 2/(sum(hann.^2))
spectrum = ft * normalization
new(freq, spectrum)
end

# Welch's method
function PSD(x::Vector, fs::Int, L::Int, overlap::AbstractFloat)
hann = hanning(L)
segs = overlaps(x, L, overlap)
map!(x -> x .* hann, segs, segs) # Apply Hanning window
map!(x -> abs2.(fft(x)), segs, segs) # Compute |H(f)|²
map!(x -> x[1:(div(L, 2) + 1)], segs, segs) # Make one sided
map!(x -> 2 .* x ./ ( sum(hann.^2)  ), segs, segs) # (2 * |H(f)|²) / (∑ wᵢ²)
w = sum(segs) ./ length(segs)
freq = [i for i in 0:(length(segs[1])-1)] .* fs/L
new(freq, w)
end
end


I used EEG data to compare my own Welch method with the one provided by the DSP.jl package. There is a lot of agreement, but there is one specific part of the spectrum which looks totally different. Here's an image comparing the two.

The blue one is my custom Welch method; the orange one is the one from the DSP package. As you can see, around frequency $$100$$ the start to diverge, and in particular there is a big downward "spike" in my custom estimation which doesn't exist in DSP's estimation.

Of course, both estimations used the same data (an EEG signal with sampling rate 500).

What could explain this difference? Is this something that is known or should I contact the developers of the package? Thanks in advance.

• Have you looked at the source code? I would debug and try to follow line by line in both your code and the source code and compare outputs for a set of data to see where yours differs from the source implementation. Commented Apr 11 at 17:41