energyop.m

function [ey,ex]=energyop(sig,gr)
%% %calculates the energy operator of a signal
%% %input
%1. Raw signal (Vector)
%2. gr (Plot or not plot)
%% %Output
%Energy operator signal (ey)
%Teager operator (ex)
%% %Method
%The Teager Energy Operator is determined as
	%(x(t)) = (dx/dt)^2+ x(t)(d^2x/dt^2) (1.1)
%in the continuous case (where x_ means the rst derivative of x, and x¨ means the second
%derivative), and as
	%[x[n]] = x^2[n] + x[n - 1]x[n + 1] (1.2)
%in the discrete case.
%% Method
%Note that the function is vectorized for optimum processing speed(Keep calm and vectorize)
%Author : Hooman Sedghamiz 
%% hoose792@student.liu.se
%%
if nargin<2
    gr=0;
end
sig=sig(:);
%% (x(t)) = (dx/dt)^2+ x(t)(d^2x/dt^2) 
%Operator 1
y=diff(sig);
y=[0;y];
squ=y(2:length(y)-1).^2;
oddi=y(1:length(y)-2);
eveni=y(3:length(y));
ey=squ - (oddi.*eveni);
%% [x[n]] = x^2[n] - x[n - 1]x[n + 1] 
%operator ex
squ1=sig(2:length(sig)-1).^2;
oddi1=sig(1:length(sig)-2);
eveni1=sig(3:length(sig));
ex=squ1 - (oddi1.*eveni1);
ex = [ex(1); ex; ex(length(sig)-2)]; %make it the same length
%% plots
if gr
figure,ax(1)=subplot(211);plot((sig/max(sig))-mean(sig/max(sig)),'b'),
hold on,
plot((ey/max(ey))-mean(ey/max(ey)),'Linewidth',2,'LineStyle','--','color','r'),
axis tight;
hleg1=legend('Original Signal','Energy Operator');
set(hleg1,'Location','NorthWest')
ax(2)=subplot(212);plot((sig/max(sig))-mean(sig/max(sig)),'b'),
hold on,
plot((ex/max(ex))-mean(ex/max(ex)),'Linewidth',2,'LineStyle','--','color','g'),
hleg2=legend('Original Signal','Teager Energy');
set(hleg2,'Location','NorthWest')
axis tight,
zoom on;
linkaxes(ax,'x');
end