Spectral Entropy Calculation in MATLAB

Question

How do I calculate the Spectral Entropy of a signal in MATLAB ? I know the basic steps but it would be nice if someone can help,

1. Calculate the power spectrum of the signal using FFT command in MATLAB.
2. Calculate the Power Spectral Density using the power spectrum or using any other technique.
3. Normalize the Power Spectral Density between $[0, 1]$, so that it can be treated as a probability density function $p_i $.
4. Calculate the Entropy $H(s) = -\sum p_i\log_2\left(p_i\right)$

Answer 1

Technically this is not a MATLAB-esque forum, but I can explain the steps in more detail for you: Suppose your input signal is $x[n]$, and its DFT is $X(f)$. For real signals you may use the one-sided DFT, since the other half would be redundant when you look at its Power Spectral Density. (PSD).

Once you compute the DFT of your signal, the PSD is simply $|X(f)|^2$. That is, you need to take the absolute magnitude of your DFT result, squared.

You now need to normalize the PSD such that it can be viewed as a Probability Density Function, (PDF). Thus, a normalized PSD, (let us call it $PSD_n$) will simply be:

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

Finally, your spectral entropy will be:

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

Answer 2

I Just do here

My source code:

    [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

This source code does Spectral Entropy calculation from every framed block ...