%
% this matlab script demonstrates the use of dHif 
%
Nw=2;  % frequency multiple (Nw < N/2)
N=21;  % signal length

t=0:N-1;             % time axis
w=2*Nw*pi/length(t); % frequency

switch 2 % select case: 0..3
case 0
  x=cos(w*t);        % signal
  H=sin(w*t);        % analytic Hilbert transform
  ixp=-w*sin(w*t);   % analytic derivative
  iyp=w*cos(w*t);    % analytic derivative of Hilbert transform
case 1
  x=sin(w*t);        % signal
  H=-cos(w*t);       % analytic Hilbert transform
  ixp=w*cos(w*t);    % analytic derivative
  iyp=w*sin(w*t);    % analytic derivative of Hilbert transform
case 2
  x=cos(w*t)+sin(w*t);        % signal
  H=sin(w*t)-cos(w*t);        % analytic Hilbert transform
  ixp=-w*sin(w*t)+w*cos(w*t); % analytic derivative
  iyp=w*cos(w*t)+w*sin(w*t);  % analytic derivative of Hilbert transform
case 3 % this case illustrates limits of instanteous frequency (note: Nw < N/4)
  a=0.4;
  x=cos(w*t)+sin(2*w*t)*a;          % signal
  H=sin(w*t)-cos(2*w*t)*a;          % analytic Hilbert transform
  ixp=-w*sin(w*t)+2*w*cos(2*w*t)*a; % analytic derivative
  iyp=w*cos(w*t)+2*w*sin(2*w*t)*a;  % analytic derivative of Hilbert transform
end

[omega,y, xp,yp]=dHif(x);  % compute Hilbert transformed-based inst. freq.

% check results
%xp_err=max(max(abs(xp-ixp)))   % check derivative accuracy
%y_err=max(max(abs(y-H)))       % check Hilbert transform accuracy
%yp_err=max(max(abs(yp-iyp)))   % check derivative accuracy

% plot signals
figure(1)
plot(t,x,'b');
hold on
plot(t,y,'r');
plot(t,H,'k.');
hold off
title('x[n]=b y[h]=r H(x(t))=k.');

% plot spectra
figure(2)
subplot(2,1,1)
plot(real(fft(x)),'b')
hold on
plot(imag(fft(x)),'r')
hold off
subplot(2,1,2)
plot(real(fft(y)),'g')
hold on
plot(real(fft(H)),'k.')
plot(imag(fft(y)),'c')
plot(imag(fft(H)),'b.')
hold off
title('Signal spectra (real=b imag=r)  Hilbert spectra (real=g imag=c)')

% plot derivative signals
figure(3)
plot(t,xp,'b');
hold on
plot(t,ixp,'b.');
plot(t,yp,'r');
plot(t,iyp,'r.');
hold off
title('Derivative signals xp[n]=b yp[h]=r');

% plot Instanteous frequency
figure(4)
plot(t,omega);
axis([0 N 0 max(max(omega))+0.1])
title(['Instanteous frequency versus time ideal=',num2str(w)]);