/*
ht.c Joe Mietus Dec 19 2001
-------------------------------------------------------------------------------
ht: Perform Hilbert transform
Copyright (C) 2002 Joe Mietus
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA.
You may contact the author by e-mail (joe at physionet dot org) or postal mail
(MIT Room E25-505A, Cambridge, MA 02139 USA). For updates to this software,
please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________
Perform Hilbert transform and calculate instantaneous magnitude and frequency.
*/
#include <stdio.h>
#include <math.h>
#define LMAX 200000
#define LFILT 128 /* impulse response filter length; must be even!!! */
main(argc, argv)
int argc;
char *argv[];
{
int i, npt, lfilt;
static double time[LMAX+1], x[LMAX+1];
static double xh[LMAX+1], ampl[LMAX+1], phase[LMAX+1], omega[LMAX+1];
double hilb[LFILT+1], pi, pi2, xt, xht;
pi = M_PI;
pi2 = 2*pi;
lfilt = LFILT;
if (argc != 1) {
usage(argv[0]);
exit(1);
}
for (i=1; i<=lfilt; i++)
hilb[i]=1/((i-lfilt/2)-0.5)/pi;
for (i=1; i<=LMAX && scanf("%lf %lf", &time[i], &x[i])==2; i++) {
xh[i] = 0.0;
ampl[i] = 0.0;
}
npt=i-1;
if (npt < lfilt) {
fprintf(stderr, "%s : insufficient data\n", argv[0]);
exit(2);
}
/* Hilbert transform */
convol(x, xh, hilb, npt, lfilt);
/* Ampl and phase */
for (i=lfilt/2+1; i<=npt-lfilt/2; i++) {
xt = x[i];
xht = xh[i];
ampl[i] = sqrt(xt*xt+xht*xht);
phase[i] = atan2(xht ,xt);
if (phase[i] < phase[i-1])
omega[i] = phase[i]-phase[i-1]+pi2;
else
omega[i] = phase[i]-phase[i-1];
}
for (i=lfilt/2+2; i<=npt-lfilt/2; i++)
printf("%g %.6f %.6f\n", time[i], ampl[i], omega[i]/pi2);
}
convol(source, target, filt, npt, lfilt)
int npt, lfilt;
double *source, *target, *filt;
{
int i, l;
double yt;
for (l=1; l<=npt-lfilt+1; l++) {
yt = 0.0;
for (i=1; i<=lfilt; i++)
yt = yt+source[l+i-1]*filt[lfilt+1-i];
target[l] = yt;
}
/* shifting lfilt/1+1/2 points */
for (i=1; i<=npt-lfilt; i++) {
target[i] = 0.5*(target[i]+target[i+1]);
}
for (i=npt-lfilt; i>=1; i--)
target[i+lfilt/2]=target[i];
/* writing zeros */
for (i=1; i<=lfilt/2; i++) {
target[i] = 0.0;
target[npt+1-i] = 0.0;
}
}
usage(prog)
char * prog;
{
fprintf(stderr, "Usage : %s\n\n", prog);
fprintf(stderr, " Reads 2 columns of input (time and x) and outputs\n");
fprintf(stderr, " time and Hilbert transform ampltudes and frequencies.\n");
}