Obliczanie spektralnej entropii w MATLAB

Jak obliczyć Entropię spektralną sygnału w MATLAB? Znam podstawowe kroki, ale byłoby miło, gdyby ktoś mógł pomóc,

1. Oblicz widmo mocy sygnału za pomocą polecenia FFT w MATLAB.
2. Oblicz gęstość widmową mocy za pomocą widma mocy lub dowolnej innej techniki.
3. Normalizuj gęstość widmową mocy pomiędzy $[0, 1]$, dzięki czemu można ją traktować jako funkcję gęstości prawdopodobieństwa $p_i$.
4. Oblicz Entropię $H(s) = -\sum p_i\log_2\left(p_i\right)$

Technicznie nie jest to forum MATLAB, ale mogę wyjaśnić ci bardziej szczegółowo następujące kroki: Załóżmy, że twój sygnał wejściowy to $x[n]$, a jego DFT to $X(f)$. Do rzeczywistych sygnałów można użyć jednostronnego DFT, ponieważ druga połowa byłaby zbędna, gdy spojrzeć na jego gęstość widmową mocy. (PSD).

Po obliczeniu DFT sygnału PSD jest po prostu $|X(f)|^2$. Oznacza to, że musisz wziąć absolutną wielkość wyniku DFT, do kwadratu.

Musisz teraz znormalizować PSD, aby można go było postrzegać jako funkcję gęstości prawdopodobieństwa (PDF). Tak więc znormalizowany PSD (nazwijmy to$PSD_n$) będzie po prostu:

$$PSD_n(f) = \frac{PSD(f)}{\sum_{f=\frac{-fs}{2}}^{f = \frac{f_s}{2}} PSD(f)}$$

Wreszcie, twoja entropia spektralna będzie:

$$E = -\sum_{f=\frac{-fs}{2}}^{f = \frac{f_s}{2}} PSD_n(f) log_2[PSD_n(f)]$$

Po prostu tu robię

Mój kod źródłowy:

    [x, Fs, nbits] = wavread('ederwander.wav'); 


    winSize = 2048;

    n_samples = length(x);


    %50% overlap or 0 to not use overlap
    OverlapStep = 50;

    if OverlapStep > 0

        Overlap = floor((OverlapStep*winSize) / 100); 
        nFrames=floor(n_samples/Overlap)-1; 
    else
        Overlap= winSize;
        nFrames=floor(n_samples/Overlap)-1;
    end

    Entropy = zeros(nFrames,1);

    k=1;
    inc=1;

    while ( (k+winSize-1) <= n_samples )

        FrameSignal = x(k:k+winSize-1);

        v = FrameSignal .* hann(length(FrameSignal));           

        N = length(v);

        Y=fft(v);

        % Compute the Power Spectrum
        sqrtPyy = ((sqrt(abs(Y).*abs(Y)) * 2 )/N);
        sqrtPyy = sqrtPyy(1:winSize/2);



       %Normalization
       d=sqrtPyy(:);
       d=d/sum(d+ 1e-12);

       %Entropy Calculation
       logd = log2(d + 1e-12);
       Entropy(inc) = -sum(d.*logd)/log2(length(d));


       k=k+Overlap;
       inc=inc+1;
end

Ten kod źródłowy wykonuje obliczenia spektralnej entropii z każdego bloku w ramce ...