image.c

speech signal process tools

	{
	    double	rsv, isv;

	    rsv = specrec->re_spec_val[j];

	    if (spectype == SPTYP_DB && fun == FN_LOG
		    || spectype == SPTYP_REAL && fun == FN_SQRT)
		data[i][j] = rsv;
	    else
	    {
		switch (spectype)
		{
		case SPTYP_PWR:
		    break;
		case SPTYP_DB:
		    rsv = pow(10.0,rsv/10.0);
		    break;
		case SPTYP_REAL:
		    rsv = rsv*rsv;
		    break;
		case SPTYP_CPLX:
		    isv = specrec->im_spec_val[j];
		    rsv = rsv*rsv + isv*isv;
		    break;
		}

		switch (fun)
		{
		case NONE:
		    data[i][j] = rsv;
		    break;
		case FN_LOG:
		    data[i][j] = 10.0 * log10(rsv);
		    break;
		case FN_EXP:
		    data[i][j] = exp(rsv);
		    break;
		case FN_SQ:
		    data[i][j] = rsv*rsv;
		    break;
		case FN_SQRT:
		    data[i][j] = sqrt(rsv);
		    break;
		}
	    }
	}
    }
    else
    {
	double	*genrec = (double *) buf;

	for (j = 0; j < nelem; j++)
	    switch (fun)
	    {
	    case NONE:
		data[i][j] = genrec[- 1 + elem_array[j]];
		break;
	    case FN_LOG:
		data[i][j] = log(genrec[- 1 + elem_array[j]]);
		break;
	    case FN_EXP:
		data[i][j] = exp(genrec[- 1 + elem_array[j]]);
		break;
	    case FN_SQ:
		data[i][j] = genrec[- 1 + elem_array[j]]
				* genrec[- 1 + elem_array[j]];
		break;
	    case FN_SQRT:
		data[i][j] = sqrt(genrec[- 1 + elem_array[j]]);
		break;
	    default:
		break;
	    }
    }

    switch (lbl_units)
    {
    case 'p':
	x[i] = tag;
	break;
    case 'r':
	x[i] = i + startrec;
	break;
    case 's':
	x[i] = src_sf_def ? (tag - 1)/src_sf
	    : start_time + (startrec + i - 1)/record_freq;
	break;
    }
}

void
reallo_data()
{
    float   **new_chunk;
    long    i, old_size;

    old_size = data_size;
    data_size += DATA_CHUNK;
    x = (double *) realloc((char *) x,
			   (unsigned) data_size * sizeof(double));
    REQUIRE(x, "can't reallocate space for x coordinate values");

    data = (float **) realloc((char *) data,
			      (unsigned) data_size * sizeof(float *));
    new_chunk = (float **) arr_alloc(2, data_dim, FLOAT, 0);
    REQUIRE(data && new_chunk, "can't extend storage for data array");
    for (i = 0; i < DATA_CHUNK; i++)
	data[old_size + i] = new_chunk[i];
    free((char *) new_chunk);
}

void
plotimage()
{
    int		    u, v;	/* pixel indices */
    long	    i;
    double	    delta_x = xmax - xmin;
    double	    delta_z = zmax - zmin;
    static int	    alloc_size = 0;
    static double   *oldrow, *midrow, *newrow;
				/* "rows" of grayscale densities */
    if (alloc_size == 0)
    {
	oldrow = (double *) malloc_d((unsigned) ncols);
	midrow = (double *) malloc_d((unsigned) ncols);
	newrow = (double *) malloc_d((unsigned) ncols);
	spsassert(oldrow && midrow && newrow,
	    "can't allocate space for 3 rows of interpolated values");
	alloc_size = ncols;

	if (debug_level >= 2)
	    Fprintf(stderr,
		"Allocated space for 3 rows of %ld interpolated values.\n",
		ncols);
    }
    else
    if (alloc_size < ncols)
    {
	oldrow = (double *) realloc((char *) oldrow,
				    (unsigned) ncols * sizeof(double));
	midrow = (double *) realloc((char *) midrow,
				    (unsigned) ncols * sizeof(double));
	newrow = (double *) realloc((char *) newrow,
				    (unsigned) ncols * sizeof(double));
	spsassert(oldrow && midrow && newrow,
	    "can't reallocate space for 3 rows of interpolated values");
	alloc_size = ncols;

	if (debug_level >= 2)
	    Fprintf(stderr,
		"Reallocated space for 3 rows of %ld interpolated values.\n",
		ncols);
    }

    if (debug_level > 1)
	Fprintf(stderr, "In plot_image xmin and xmax: %g, %g\n", xmin, xmax);

    if (Bflag) axes_and_titles();
    dot_init();
    (*dev_initbits)();

    if (Dflag)
    {
	u = 0;
	for (i = 0; i < nrecs; i++)
	{
	    pick_nearest(data[i], y, nelem, midrow, ncols);

	    if (Gflag)
	    {
		double	rowmax = -FLT_MAX;
		double	rowmin;
		
		for (v = 0; v < ncols; v++)
		    if (rowmax < midrow[v]) rowmax = midrow[v];

		gain_adj(pwr_flag, &rowmin, &rowmax);

		for (v = 0; v < ncols; v++)
		    midrow[v] = (midrow[v] - rowmin)/(rowmax - rowmin);
	    }
	    else
	    {
		for (v = 0; v < ncols; v++)
		    midrow[v] = (midrow[v] - zmin) / delta_z;
	    }
	    
	    for ( ;
		 u <= nrows - 1
		 && (i == nrecs - 1
		     || u < (nrows-1) * (0.5*(x[i]+x[i+1]) - xmin) / delta_x);
		 u++
		)
	    {
		dot_row(midrow);
	    }
	}
    }
    else
    {
	interp(data[0], y, nelem, newrow, ncols);
	u = 0;
	for (i = 1; i < nrecs; i++)
	{
	    double	*tmp;

	    tmp = oldrow; oldrow = newrow; newrow = tmp;

	    interp(data[i], y, nelem, newrow, ncols);

	    for ( ;
		 u <= nrows - 1
		 && (i == nrecs - 1
		     || u <= (nrows-1) * (x[i]-xmin) / delta_x);
		 u++
		)
	    {
		double  alpha =
		    (u*delta_x/(nrows-1.0) + xmin - x[i-1]) / (x[i] - x[i-1]);
		double  beta = 1.0 - alpha;
		double  gamma;

		if (Gflag)
		{
		    double	rowmax = -FLT_MAX;
		    double	rowmin;

		    for (v = 0; v < ncols; v++)
		    {
			midrow[v] =
			    alpha*newrow[v] + beta*oldrow[v];
			if (rowmax < midrow[v]) rowmax = midrow[v];
		    }

		    gain_adj(pwr_flag, &rowmin, &rowmax);

		    for (v = 0; v < ncols; v++)
			midrow[v] = (midrow[v] - rowmin)/(rowmax - rowmin);
		}
		else
		{
		    alpha /= delta_z;
		    beta /= delta_z;
		    gamma = zmin / delta_z;

		    for (v = 0; v < ncols; v++)
			midrow[v] =
			    alpha*newrow[v] + beta*oldrow[v] - gamma;
		}

		dot_row(midrow);
	    }
	}
    }

    dot_fin();
    dev_fin();
}

/*
 * Given values dat[j] of a function at points y[j] (j = 0, ..., (nelem-1)),
 * determine values row[v] (v = 0, ..., (ncols-1)) by linear interpolation
 * at ncols other points equispaced from y[0] to y[nelem-1].
 */

void
interp(dat, y, nelem, row, ncols)
    float   *dat;
    double  *y;
    long    nelem;
    double  *row;
    long    ncols;
{
    long    j;
    long    v;
    double  delta_y;

    delta_y = y[nelem - 1] - y[0];
    v = 0;
    for (j = 1; j < nelem; j++)
    {
	for ( ; v <= (ncols-1)*(y[j]-y[0])/delta_y; v++)
	{
	    double  beta =
		(v*delta_y/(ncols-1.0) + y[0] - y[j-1]) / (y[j] - y[j-1]);

	    row[v] = beta*dat[j] + (1.0 - beta)*dat[j-1];
	}
    }
    if (v < ncols)
	row[v] = dat[j];
}

/*
 * Like  interp  above, but the value at a point in the interval
 * (y[j-1], y[j]), instead of being linearly interpolated between dat[j-1]
 * and dat[j], is simply one or the other of those values, depending on
 * which end of the interval is closer.  The interpolating function,
 * instead of being linear on the interval, is constant on each half of
 * the interval with a jump at the center.
 */

void
pick_nearest(dat, y, nelem, row, ncols)
    float   *dat;
    double  *y;
    long    nelem;
    double  *row;
    long    ncols;
{
    long    j;
    long    v;
    double  delta_y;

    delta_y = y[nelem - 1] - y[0];
    v = 0;
    for (j = 1; j < nelem; j++)
	for ( ; v < (ncols-1) * (0.5*(y[j-1]+y[j]) - y[0]) / delta_y; v++)
	    row[v] = dat[j-1];
    for ( ; v < ncols; v++)
	row[v] = dat[nelem-1];
}

int (*
read_cmap(fname, len) )[3]
    char	*fname;
    int		len;
{
    char    	*grange;
    long    	*elements, num_elts;
    FILE	*file;
    int		(*map)[3];	/* colormap entries read from file */
    int		(*imap)[3];	/* interpolated colormap entries */
    int		num_lines;
    unsigned	alloc_size = 16;
    char	line[256];
    long	i, j;
    double  	alpha, ratio;

    fname = strtok(fname, "[");
    grange = strtok((char *) 0, "]");
    elements =  grange ? grange_switch(grange, &num_elts) : NULL;

    if (!fname) return NULL;

    file = fopen(fname, "r");
    if (!file)
    {
	Fprintf(stderr, "%s: can't open colormap file %s.\n", ProgName, fname);
	return NULL;
    }

    map = (int (*)[3]) malloc((unsigned) alloc_size * sizeof(int [3]));
    spsassert(map, "can't allocate space for colormap");

    for (num_lines = 0; fgets(line, 256, file); num_lines++)
    {
	if (num_lines == alloc_size)
	{
	    alloc_size *= 2;
	    map = (int (*)[3])
		  realloc((char *) map,
			  (unsigned) alloc_size * sizeof(int [3]));
	    spsassert(map, "can't reallocate space for colormap");
	}
	sscanf(line, "%d%d%d",
	       &map[num_lines][0], &map[num_lines][1], &map[num_lines][2]);
    }

    if (elements)
    {
	int 	(*submap)[3];
	long	e;

	submap = (int (*)[3]) malloc((unsigned) num_elts * sizeof(int[3]));
	spsassert(submap, "can't allocate space for subset of colormap");

	for (i = 0; i < num_elts && (e = elements[i]) < num_lines; i++)
	{
	    submap[i][0] = map[e][0];
	    submap[i][1] = map[e][1];
	    submap[i][2] = map[e][2];
	}
	if (i == num_elts)
 	    num_lines = num_elts;
	else
	{
	    Fprintf(stderr, "%s: no line %ld in colormap file.\n",
		    ProgName, e);
	    free((char *) submap);
	    submap = NULL;
	    num_lines = 0;
	}

	free((char *) map);
	map = submap;
	free((char *) elements);
    }

    if (debug_level >= 2)
    {
	Fprintf(stderr, "%d colormap entries\n", num_lines);
	for (i = 0; i < num_lines; i++)
	    Fprintf(stderr, "%ld:\t%d %d %d\n",
		    i, map[i][0], map[i][1], map[i][2]);
    }

    if (num_lines < 2)
    {
	if (map) free((char *) map);
	return NULL;
    }

    if (num_lines == len) return map;

    imap = (int (*)[3]) malloc((unsigned) len * sizeof(int [3]));
    spsassert(map, "can't allocate space to interpolate colormap");

    ratio = (double) (num_lines - 1) / (len - 1);

    for (i = 0; i < len; i++)
    {
	alpha = i * ratio;
	j = alpha;
	if (j == num_lines - 1) j--;
	alpha -= j;

	imap[i][0] = 0.5 + alpha * map[j+1][0] + (1.0 - alpha) * map[j][0];
	imap[i][1] = 0.5 + alpha * map[j+1][1] + (1.0 - alpha) * map[j][1];
	imap[i][2] = 0.5 + alpha * map[j+1][2] + (1.0 - alpha) * map[j][2];
    }

    return imap;
}

void
gain_adj(flag, rowmin, rowmax)
    int	    flag;
    double  *rowmin, *rowmax;
{
    if (flag)
    {
	switch (fun)
	{
	case NONE:
	    *rowmin = *rowmax * pow(10.0, gainlolim/10.0);
	    *rowmax = *rowmax * pow(10.0, gainhilim/10.0);
	    if (*rowmax > zmax)
	    {
		*rowmin *= zmax / *rowmax;
		*rowmax = zmax;
	    }
	    if (*rowmin < zmin)
	    {
		*rowmax *= zmin / *rowmin;
		*rowmin = zmin;
	    }
	    break;
	case FN_LOG:
	    *rowmin = *rowmax + gainlolim;
	    *rowmax = *rowmax + gainhilim;
	    if (*rowmax > zmax)
	    {
	    	*rowmin += zmax - *rowmax;
		*rowmax = zmax;
	    }
	    if (*rowmin < zmin)
	    {
		*rowmax += zmin - *rowmin;
		*rowmin = zmin;
	    }
	    break;
	case FN_EXP:
	    *rowmin = pow(*rowmax, pow(10.0, gainlolim/10.0));
	    *rowmax = pow(*rowmax, pow(10.0, gainhilim/10.0));
	    if (*rowmax > zmax)
	    {
		*rowmin = pow(zmax, pow(10.0, gainlolim/10.0
				       - gainhilim/10.0 ));
		*rowmax = zmax;
	    }
	    if (*rowmin < zmin)
	    {
		*rowmax = pow(zmin, pow(10.0, gainhilim/10.0
				       - gainlolim/10.0 ));
		*rowmin = zmin;
	    }
	    break;
	case FN_SQ:
	    *rowmin = *rowmax * pow(10.0, gainlolim/5.0);
	    *rowmax = *rowmax * pow(10.0, gainhilim/5.0);
	    if (*rowmax > zmax)
	    {
		*rowmin *= zmax / *rowmax;
		*rowmax = zmax;
	    }
	    if (*rowmin < zmin)
	    {
		*rowmax *= zmin / *rowmin;
		*rowmin = zmin;
	    }
	    break;
	case FN_SQRT:
	    *rowmin = *rowmax * pow(10.0, gainlolim/20.0);
	    *rowmax = *rowmax * pow(10.0, gainhilim/20.0);
	    if (*rowmax > zmax)
	    {
		*rowmin *= zmax / *rowmax;
		*rowmax = zmax;
	    }
	    if (*rowmin < zmin)
	    {
		*rowmax *= zmin / *rowmin;
		*rowmin = zmin;
		}
	    break;
	}
    }
    else
    {
	*rowmin = *rowmax + gainlolim;
	*rowmax = *rowmax + gainhilim;
	if (*rowmax > zmax)
	{
	    *rowmin += zmax - *rowmax;
	    *rowmax = zmax;
	}
	if (*rowmin < zmin)
	{
	    *rowmax += zmin - *rowmin;
	    *rowmin = zmin;
	}
    }
}
/*
 * for debug printout of long-integer arrays
 */

void
pr_larray(text, n, arr)
    char    *text;
    int     n;
    long    *arr;
{
    int     i;

    Fprintf(stderr, "%s - %d points:\n", text, n);
    for (i = 0; i < n; i++)
    {
        Fprintf(stderr, "%ld ", arr[i]);
        if (i%5 == 4 || i == n - 1) Fprintf(stderr, "\n");
    }
}