Estimation of NegEntropy in MATLAB

7 次查看（过去 30 天）
Marimuthu Ananthavelu
回答： William Rose 2021-5-28
I am trying to evaluate the different ICA algorithms using MATLAB. To do that, one of the measure which I use is to estimate the non-gaussianity using NegEntropy. I am trying to find a function which can do this on time series EEG data. Would appreciate your help.
Thanks
Kind regards,
Mari

请先登录，再回答此问题。

回答（2 个）

William Rose
William Rose 2021-5-27
编辑：William Rose 2021-5-27
Ran in:
@Marimuthu Ananthavelu and @Smriti Thakur,
pentropy(), in the SIgnal Processing Toolbox, computes the entropy of a signal's power spectrum, or of its spectrogram (i.e. time-varying spectrum). By default, it divides the signal into segments and computes a spectrum and an entropy for each segment. This can be useful in analysis of speech signals. If you want to estimate the entropy (H) of the entire signal, without segmentation, specify 'Instantaneous',false, as shown below:
y=(rand(1,256)-.5)*sqrt(12); %uniform random noise with mean=0, st.dev.=1
H=pentropy(y,1,'Instantaneous',false)
H = 0.9965
The entropy returned is normalized, by default, to the entropy of a white noise signal of equal length. You can test pentropy() by calculating H for a Gaussian white noise signal. In the following example, 100 Gaussian random signals are generated, each of length N=256. The entropy of each is computed and saved. At the end, the mean, standard deviation, minimum, amd maximum entropies are displayed.
N=256; H=zeros(1,100);
for i=1:100,H(i)=pentropy(randn(1,N),1,'Instantaneous',false);end
disp([mean(H),std(H),min(H),max(H)]);
0.9974 0.0017 0.9884 0.9995
As the signal duration becomes long, H approaches unity. See example below, which is like above, but with N=32768 this time.
N=32768; H=zeros(1,100);
for i=1:100,H(i)=pentropy(randn(1,N),1,'Instantaneous',false);end
disp([mean(H),std(H),min(H),max(H)]);
0.9998 0.0000 0.9997 0.9999
A value of H less than unity represents non-Gaussianity. The smaller H is, the more non-Gaussian y is.
pentropy() is not sensitive to some important forms of non-Gaussianity. The spectral entropy of uniform random noise (in the first example above) is almost identical to Gaussian random noise. This should not happen, because in theory, a Gaussian signal has maximum entropy.
The difficulty in controlling how the spectrum is computed inside pentropy() is another potential disadvatage of pentropy(). When I turned off the automatic scaling, and varied the signal duration, the entropy of the supposedly unscaled signal did not change as I expected. I think this is because the power spectrum routine called by pentropy() has default settings for the number of frequencies in the spectrum. See examples below.
N=4096;
H=pentropy(randn(1,N),1,'Instantaneous',false,'Scaled',false)
H = 9.9921
I expected H=log2(N)=12.
N=256;
H=pentropy(randn(1,N),1,'Instantaneous',false,'Scaled',false)
H = 9.9780
I expected H=log2(N)=8.
The shortcomings of pentropy() show that a time-domain entropy function is needed. This function will construct the histogram of the time-domain values of the signal, then normalize the histogram so it is a probability density, then compute the entropy of the density. A search of Matlab file exchange for negentropy returns two hits, but neither one looks like what we want.
William Rose
William Rose 2021-5-28
@Marimuthu Ananthavelu and @Smriti Thakur,
Here is code to estimate the entropy and negentropy of time domain signal, y. In this example, y is uniform random noise, but y could be any signal. The units for Hent(y) and J(y) are bits.
N=1024; %signal length
y=sqrt(12)*(rand(1,N)-.5); %uniform random noise, mean=0, stdev=1
h=histogram(y,'Normalization','pdf'); %estimate f(x)=the p.d.f.
Hent=-h.BinWidth*sum(h.Values(h.Values>0).*log(h.Values(h.Values>0))/log(2)); %compute entropy
J=log(std(y))/log(2)+2.0471-Hent; %compute negentropy
fprintf('Uniform random noise: Hent=%.4f, J=%.4f\n',Hent,J); %display results
The justification for the code above is in the attached document.

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Estimation 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by