cleaning up the code

otoolej · Apr 18, 2014 · 51cc57e · 51cc57e
1 parent e61235a
commit 51cc57e
Show file tree

Hide file tree

Showing 10 changed files with 292 additions and 124 deletions.
diff --git a/README.md b/README.md
@@ -1,43 +1,103 @@
-# Methods to Estimate Instantaneous Energy (including the nonlinear energy operator)
+Estimating Instantaneous Energy 
+================================
 
-
-Collection of M-files (computer code) to implements instantaneous energy measures, as
+Collection of M-files (computer code) to implement instantaneous energy measures, including the "nonlinear energy operator", as
 described in [[1]](#references).
-Requires Matlab or Octave (programming environments).
+Requires Matlab or Octave programming environments.
 
 
 # Contents
-- [quick start](#quick-start)
-- [requirements](#requirements)
-- [test computer setup](#test-computer-setup)
-- [licence](#licence)
-- [references](#references)
-- [contact](#contact)
+* [overview](#overview)
+* [quick start](#quick-start)
+* [requirements](#requirements)
+* [test computer setup](#test-computer-setup)
+* [licence](#licence)
+* [references](#references)
+* [contact](#contact)
+
+
+# overview
+Implements methods to estimate frequency-weighted instantaneous energy.  Implements the
+Teager–Kaiser operator, often referred to as the nonlinear energy operator, and a similar
+frequency-weight operator proposed in reference [[1]](#references). The Teager–Kaiser operator is
+simply defined, for discrete signal x(n), as
+```
+Ψ[x(n)] = x²(n) - x(n+1)x(n-1)
+```  
+and the proposed energy measure is defined as
+```
+ Γ[x(n)] = y²(n) + H[y(n)]²
+```
+where y(n) is the derivative of x(n), estimated using the central-finite difference
+equation y(n)=[x(n+1)-x(n-1)]/2, and H[·] is the discrete Hilbert transform of $x(n)$.  Reference [[1]](#references) contains more details.
 
 
 # quick start
-TO DO....
 
-### files
+The following example generates the Teager–Kaiser operator and the proposed
+envelope–derivative operator for a test signal (sum of two sinusoidals signals),
+cut-and-paste the following code into Matlab (or Octave):
+```matlab
+  % generate two sinusoidal signals:
+  N=256; n=0:N-1;
+  w1=pi/(N/32); ph1=-pi+2*pi*rand(1,1);  a1=1.3;
+  w2=pi/(N/8); ph2=-pi+2*pi*rand(1,1);  a2=3.1;
+  x1=a1.*cos(w1.*n + ph1);  x2=a2.*cos(w2.*n + ph2);
+  x=x1+x2;
+
+  % compute instantaneous energy:
+  x_env_diff=cal_freqweighted_energy(x,1,'envelope_diff');
+  x_teager  =cal_freqweighted_energy(x,1,'teager');
+
+  % plot:
+  figure(1); clf; 
+  subplot(211); hold all; plot(x); ylabel('amplitude');
+  subplot(212); hold all; plot(x_env_diff,'-'); plot(x_teager,'--');
+  ylabel('energy');
+  legend('envelope-derivative','Teager-Kaiser');
+```
+
+
+## properties
+To test the properties of the operators with some example signals, call
+`properties_test_Hilbert_NLEO(number)` with argument `number` either 0,1,2,3, or 4. For
+example, 
+```matlab
+  >> properties_test_Hilbert_NLEO(3);
+```
+calls the function with frequency modulated signal with instantaneous frequency law of
+0.1+0.3sin(tπ/N).
+
+## noise Analysis
+To compare the bias for each method, run the function
+```matlab
+  >> bias_of_estimators;
+```
+which computes the mean-value (and therefore an approximation to the Expectation operator)
+of 10,000 iterations of white Gaussian noise. This then produces Fig. 2 in the ‘Noise
+Analysis’ section of [[1]](#references).
+
+
+# files
 All Matlab files (.m files) have a description and an example in the header. To read this
 header, type `help <filename.m>` in Matlab.  Directory structure is as follows: 
 ```
 .
-├── bias_of_estimators.m
-├── cal_freqweighted_energy.m
-├── compare_nleo_methods.m
-├── discrete_Hilbert_diff_operator.m
-├── do_bandpass_filtering.m
-├── general_nleo.m
-├── nleo_parameters.m
-├── pics/                              # directory for figures
-├── plot_eeg_examples.m
-├── properties_test_Hilbert_NLEO.m
-└── script_test_eeg_data.m
+├── bias_of_estimators.m              # noise analysis: estimate bias with WGN
+├── cal_freqweighted_energy.m         # select method to estimate instantaneous energy
+├── discrete_Hilbert_diff_operator.m  # proposed envelope–derivative operator
+├── do_bandpass_filtering.m           # simply band-pass filtering
+├── general_nleo.m                    # general Nonlinear Energy Operator (Plotkin–Swamy)
+├── nleo_parameters.m                 # set parameters here (directions etc.)
+├── pics/                             # directory for figures
+├── properties_test_Hilbert_NLEO.m    # test properties of the different operator
+└── requires_EEG_data		          # directory containing files for EEG analysis
+    ├── compare_nleo_methods.m		  # these files require EEG data to run.
+    ├── plot_eeg_examples.m	
+    └── script_test_eeg_data.m
 ```
 
 
-
 # requirements
 Either Matlab (R2012 or newer,
 [Mathworks website](http://www.mathworks.co.uk/products/matlab/)) or Octave (v3.6 or
@@ -48,7 +108,7 @@ newer, [Octave website](http://www.gnu.org/software/octave/index.html), with the
 
 # test computer setup
 - hardware:  Intel(R) Xeon(R) CPU E5-1603 0 @ 2.80GHz; 8GB memory.
-- operating system: Ubuntu GNU/Linux x86_64 distribution (Raring, 13.04), with Linux kernel 3.5.0-28-generic 
+- operating system: Ubuntu GNU/Linux x86_64 distribution (Raring, 13.04), with Linux kernel 3.11.0-19-generic
 - software: Octave 3.6.4 (using Gnuplot 4.6 patchlevel 1), with 'octave-signal' toolbox and Matlab (R2009b, R2012a, and R2013a)
 
 ---
@@ -88,6 +148,8 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 # references
 
+1. J.M. O' Toole and N.J. Stevenson, “Assessing instantaneous energy in the EEG: a
+non-negative, frequency-weighted energy operator”, under review, 2014
 
 ---
 

diff --git a/bias_of_estimators.m b/bias_of_estimators.m
@@ -1,10 +1,11 @@
 %-------------------------------------------------------------------------------
-% bias_of_estimators: assess whether estimators are biased to noise are not
+% bias_of_estimators: Estimate bias of operators by approximating expectation 
+% operator. Using 10,000 iterations of white Gaussian noise.
 %
-% Syntax: []=bias_of_estimators()
+% Syntax: []=bias_of_estimators(PRINT_PLOTS)
 %
 % Inputs: 
-%      - 
+%      PRINT_PLOTS - print or not, either 0 or 1; 
 %
 % Outputs: 
 %     [] - 
@@ -16,55 +17,72 @@
 % John M. O' Toole, University College Cork
 % Started: 28-03-2014
 %
-% last update: Time-stamp: <2014-04-07 00:45:52 (otoolej)>
+% last update: Time-stamp: <2014-04-18 11:25:37 (otoolej)>
 %-------------------------------------------------------------------------------
 function []=bias_of_estimators(PRINT_PLOTS)
 if(nargin<1 || isempty(PRINT_PLOTS)), PRINT_PLOTS=0; end
 
+
+%---------------------------------------------------------------------
+% 0. set parameters
+%---------------------------------------------------------------------
 nleo_parameters;
 FONT_NAME='Arial';
 FONT_SIZE=12;
 
-
+% number of iterations (L) and signal length (N) of white Gaussian noise:
 N=64; L=10000;
 xx=randn(L,N);
 
+% set to 1 to take absolute value of the operator:
 ABS_NLEO=0;
+% set to 1 to up-sampling signal before analysis:
 RESAMPLE_X=1;
 
 Nup=N*4; Nup2=ceil(N*6.33);
 nplain_teag=zeros(L,Nup);
 nenv_diff=zeros(L,Nup);
 nplain_got=zeros(L,Nup2);
 
+
+%---------------------------------------------------------------------
+% 1. iterate analysis and estimate Expectation operator by taking the 
+%    mean value
+%---------------------------------------------------------------------
 for n=1:L
+
+    % up-sample by a factor of 4 for the Teager--Kaiser and envelope-derivate 
+    % operators:
     if(RESAMPLE_X)
         xxf=resample(xx(n,:),ceil(4*N),N);
     else
         xxf=xx(n,:);
     end
     xxf=xxf./std(xxf);
 
+    % compute the Teager--Kaiser and envelope-derivate operators:
     nenv_diff(n,:)=cal_freqweighted_energy(xxf,1,'envelope_diff');
     nplain_teag(n,:)=cal_freqweighted_energy(xxf,1,'teager');
-
-% $$$     nplain_teag(n,:)=plain_nleo(xxf,1,0,[],1);    
-% $$$     nenv_diff(n,:)=discrete_Hilbert_diff_operator(xxf,1,0,[]);
 
+
+    % up-sample by a factor of 6.33 for the Agarwal--Gotman operator:
     if(RESAMPLE_X)
         xxf=resample(xx(n,:),ceil(6.33*N),N);
         xxf=xxf./std(xxf);        
     end
 
+    % compute the Agarwal--Gotman operator:    
     nplain_got(n,:)=cal_freqweighted_energy(xxf,1,'agarwal');
 end
-
 if(ABS_NLEO)
     nplain_teag=abs(nplain_teag);
     nplain_got=abs(nplain_got);        
 end
 
 
+%---------------------------------------------------------------------
+% 2. PLOT the expectation values (i.e. mean-values)
+%---------------------------------------------------------------------
 Nteag=size(nplain_teag,2);
 npart=4:Nteag-4;
 Ngot=size(nplain_got,2);
@@ -101,14 +119,14 @@
 end
 
 
-return;
-figure(2); clf; hold all;
-subplot(221); plot( nplain_teag.' );
-subplot(222); plot( nplain_got.' );
-subplot(223); plot( nenv_diff.' );
-subplot(224); plot( xxf );
-
-figure(3); clf; hold all;
-pfft(xxf,1);
+% $$$ return;
+% $$$ figure(2); clf; hold all;
+% $$$ subplot(221); plot( nplain_teag.' );
+% $$$ subplot(222); plot( nplain_got.' );
+% $$$ subplot(223); plot( nenv_diff.' );
+% $$$ subplot(224); plot( xxf );
+% $$$ 
+% $$$ figure(3); clf; hold all;
+% $$$ pfft(xxf,1);
 
 
diff --git a/cal_freqweighted_energy.m b/cal_freqweighted_energy.m
@@ -7,22 +7,41 @@
 % Inputs: 
 %     x           - input signal
 %     Fs          - sampling frequency
-%     method      - either 'teager','arg',...
+%     method      - either 'teager','agarwal', 'envelope_diff', 'palmu', 'abs_teager',
+%                   or 'env_only'
 %     wlength_ma  - length (in samples) of output moving-average filter; set to [] to
 %                   turn off
 %     bandpass_filter_params - band-pass filter values; set to [] to turn off
 %
 % Outputs: 
-%     [x_nleo] - 
+%     [x_nleo] - output operator
 %
 % Example:
-%     
 %
+%     % generate two sinusoidal signals:
+%     N=256; n=0:N-1;
+%     w1=pi/(N/32); ph1=-pi+2*pi*rand(1,1);  a1=1.3;
+%     w2=pi/(N/8); ph2=-pi+2*pi*rand(1,1);  a2=3.1;
+%     x1=a1.*cos(w1.*n + ph1);  x2=a2.*cos(w2.*n + ph2);
+%     x=x1+x2;
+%
+%     % compute instantaneous energy:
+%     x_env_diff=cal_freqweighted_energy(x,1,'envelope_diff');
+%     x_teager  =cal_freqweighted_energy(x,1,'teager');
+%
+%     % plot:
+%     figure(1); clf; 
+%     subplot(211); plot(x); ylabel('amplitude');
+%     subplot(212); hold all; plot(x_env_diff,'-'); plot(x_teager,'--');
+%     ylabel('energy');
+%     legend('envelope-derivative','Teager-Kaiser');
+%
+
 
 % John M. O' Toole, University College Cork
 % Started: 04-04-2014
 %
-% last update: Time-stamp: <2014-04-06 02:36:45 (otoolej)>
+% last update: Time-stamp: <2014-04-18 12:23:59 (otoolej)>
 %-------------------------------------------------------------------------------
 function [x_nleo]=cal_freqweighted_energy(x,Fs,method,wlength_ma, ...
                                                  bandpass_filter_params)

diff --git a/discrete_Hilbert_diff_operator.m b/discrete_Hilbert_diff_operator.m
@@ -4,24 +4,37 @@
 % transform (of length-N) for Hilbert transform
 %
 % Γ[x(n)] = y(n)² + H[y(n)]², where y(n) is the derivative of x(n) and H[] is the
-% Hilber transform
+% Hilbert transform
 %
 % Syntax: [x_nleo]=discrete_Hilbert_diff_operator(x)
 %
 % Inputs: 
 %     x - input signal
 %
 % Outputs: 
-%     [x_nleo] - 
+%     x_nleo - envelope-derivative operator (i.e. energy measure)
 %
 % Example:
+%     % generate test signals:
+%     N=256; n=0:N-1;
+%     w1=pi/(N/32); ph1=-pi+2*pi*rand(1,1);  a1=1.3;
+%     x1=a1.*cos(w1.*n + ph1);  
+%
+%     % compute instantaneous energy:
+%     x_env_diff=discrete_Hilbert_diff_operator(x1);
+%
+%     % plot:
+%     figure(1); clf; 
+%     subplot(211); plot(x1); ylabel('amplitude');
+%     subplot(212); plot(x_env_diff,'-'); ylim([0 0.5])
+%     ylabel('energy');
 %     
 %
 
 % John M. O' Toole, University College Cork
 % Started: 03-04-2014
 %
-% last update: Time-stamp: <2014-04-04 16:35:12 (otoolej)>
+% last update: Time-stamp: <2014-04-18 13:09:57 (otoolej)>
 %-------------------------------------------------------------------------------
 function [x_nleo]=discrete_Hilbert_diff_operator(x)
 
@@ -56,6 +69,9 @@
 x_nleo=x_nleo(1:Nstart);
 
 
+% Octave includes imaginary quantization noise:
+x_nleo=real(x_nleo);
+
 %---------------------------------------------------------------------
 % DEBUG: testing and compare 
 %---------------------------------------------------------------------

diff --git a/general_nleo.m b/general_nleo.m
@@ -5,14 +5,28 @@
 % Syntax: x_nleo=general_nleo(x,l,p,q,s)
 %
 % Inputs: 
-%     x,l,p,q,s - 
+%     x       - input signal
+%     l,p,q,s - parameters for the operator; integer values and must satisfy: 
+%                 l+p=q+s and [l,p]≠[q,s]
 %
 % Outputs: 
-%     x_nleo - 
+%     x_nleo - output nonlinear energy operator
 %
 % Example:
-%     
+%     % generate test signals:
+%     N=256; n=0:N-1;
+%     w1=pi/(N/32); ph1=-pi+2*pi*rand(1,1);  a1=1.3;
+%     x1=a1.*cos(w1.*n + ph1);  
 %
+%     % compute instantaneous energy:
+%     x_nleo=general_nleo(x1,1,2,0,3);
+%
+%     % plot:
+%     figure(1); clf; 
+%     subplot(211); plot(x1); ylabel('amplitude');
+%     subplot(212); plot(x_nleo,'-'); ylim([0 0.5])
+%     ylabel('energy');
+
 
 % John M. O' Toole, University College Cork
 % Started: 28-01-2014