0
$\begingroup$

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.

enter image description here

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.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Baddioes
    Commented Apr 11 at 17:41

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.