www.pudn.com > 医学算法.rar > hf.c

#include	
#include	

/*
	hf - homomorphic filt UCT waveforms

	Carl Crawford
	Purdue University
	W. Lafayette, IN 47907

	Feb. 13, 1981
*/

#define	L	128	/* half length of spectrum */
#define L2      (L + L)
#define	M	8	/* log base 2 of 2L */
#define MN	256	/* maximum n */

int	n = 100;	/* number of samples per projection */
int	k = 1;		/* number of projections */
float	gam = 1.0;	/* gamma */
float	gamp = 1.0;	/* gamma prime */

char    *ifn = "udata";  /* input file name */
char	*ofn = "for003.dat"; /* output file name */
char    *sfn = "spectrum"; /* save magnitude of spectrums */
int     spec;           /* 1=save spectrum */
float	proj[MN];	/* projection */
float   a[L2];          /* real part of fft array */
float   b[L2];          /* imaginary part of fft array */
float   aa[L2];		/* real part of fft array */
float   bb[L2];		/* imaginary part of fft array */
int     tsig[L2];       /* time signal */
FILE    *ifd,*ofd,*sfd,*cfd,*ssfd; /* file descriptors */
int	ssspec;		/* save final spectrum */
float   m0,m1;          /* moments */
int     nofilt;         /* 1=no filter */
int     time;           /* 1=data in time domain */
float   tmp;            /* temp var */
int     b1 = 10;	/* band for division */
int     b2 = 50;        /* band for division */
float	sigma = 3.0;	/* square of standard deviation */
int     div;            /* 1=spectral division */
float   c0,c1;          /* constants for spectral division */
float   df;             /* frequency division */
float   dt = .05;       /* time division */
float   fg[L];          /* f to the power gamma */
float   k11,k12,k22;    /* constants for div */
int     band = 20;      /* part of cepstrum kept */
int     print = 1;      /* 1 = print info */
int     swp;            /* indicator for swap subroutine */
int	raw;		/* 1= raw mode */
int	sw;		/* 1=don't swap */
int	cep;		/* 1=save cepstrum */
char	*cfn= "cepstrum";	/* cepstrum file name */
int	logt;		/* 1=log taken of power spectrum */

char    *table[] = {
		"n",
		"of",
		"if",
		"k",
		"sf",
		"gammap",
		"sigma",
		"b1",
		"b2",
		"band",
		"gamma",
		"dt",
		"cf",
		0
	};

main(argc,argv)
	int	argc;
	char	**argv;
{
	char cc;
	register ii,j,i;
	float fv();

	while(argv++ , --argc){
		if(argv[0][0] == '-')while(cc = *++*argv)switch(cc){

			case 'p':	/* print flag */
				print = 0;
				break;
			case 'd':       /* spectral division */
				div = 1;
				break;
			case 'n':       /* no filtering */
				nofilt = 1;
				break;
			case 't':       /* data in time domain */
				time = 1;
				break;
			case 's':       /* save spectrum */
				spec = 1;
				break;
			case 'S':	/* save final spectrum */
				ssspec = 1;
				break;
			case 'c':	/* save cepstrum */
				cep = 1;
				break;
			case 'r':	/* raw mode */
				raw = 1;
				break;
			case 'w':	/* don't swap */
				sw = 1;
				break;
			default:
				fprintf(stderr,"bad flag: -%c\n",cc);
				exit(1);
			}
		else
			switch(i = comm(*argv)){

			case 1:
				n = atoi(*argv + 2);
				break;
			case 2:
				ofn = *argv + 3;
				break;
			case 3:
				ifn = *argv + 3;
				break;
			case 4:
				k = atoi(*argv + 2);
				break;
			case 5:
				sfn = *argv + 3;
				spec = 1;
				break;
			case 6:
				gamp = atof(*argv + 7);
				break;
			case 7:
				sigma = atof(*argv + 6);
				break;
			case 8:
				b1 = atoi(*argv + 3);
				break;
			case 9:
				b2 = atoi(*argv + 3);
				break;
			case 10:
				band = atoi(*argv + 5);
				break;
			case 11:
				gam = atof(*argv + 6);
				break;
			case 12:
				dt = atof(*argv + 3);
				break;
			case 13:
				cfn = *argv +3;
				break;
		}
	}

	if(n <= 0 || n > MN)err("illegal n","");
	if(ifn[0] == '-' && !ifn[1]){
		ifd = stdin;
	}else{
		if((ifd = fopen(ifn,"r")) == NULL)err("can't open: ",ifn);
	}
	if((ofd = fopen(ofn,"w")) == NULL)err("can't create: ",ofn);
	if(spec)if((sfd = fopen(sfn,"w")) == NULL)err("can't create: ",sfn);
	if(cep)if((cfd = fopen(cfn,"w")) == NULL)err("can't create: ",cfn);
	if(ssspec)if((ssfd = fopen("fspec","w")) == NULL)err("can't create: ","fspec");

	if(print){
		printf("HF DATA\n\n");
		printf("n: %d\n",n);
		printf("k: %d\n",k);
		printf("sigma: %f\n",sigma);
		printf("gamma: %f\n",gam);
		printf("gammap: %f\n",gamp);
		printf("dt: %f\n",dt);
		if(div){
			printf("b1: %d\n",b1);
			printf("b2: %d\n",b2);
		}
		printf("band: %d\n",band);
		printf("input file: %s\n",(ifd == stdin)? "stdin" : ifn);
		printf("output file: %s\n",ofn);
		printf("%s technique\n",div?"division" : "fs");
		if(time)printf("time domain data\n");
		if(nofilt)printf("no homomorphic filtering\n");
		if(spec)printf("spectrum saved in: %s\n",sfn);
		if(cep)printf("cepstrum saved in: %s\n",cfn);
		if(ssspec)printf("final spectrum saved in: fspec\n");
		if(raw)printf("raw mode\n");
		if(sw)printf("no projection swapping\n");
	}


	df = 1.0 / (L2 * dt);
	if(!div){
		b1 = 3;
		b2 = L - 1;
	}
	k11 = k12 = k22 = 0.0;
	for(i=b1;i<=b2;i++){
		fg[i] = pow(i*df,gam);
		k11 += 1;
		k12 += fg[i];
		k22 += fg[i] * fg[i];
	}
	c0 = k11 / (k11 * k22 - k12 * k12) / 2.0;
	c1 = k12 / (k11 * k22 - k12 * k12) / 2.0;

	for(ii=0;ii