spectrogram.c

/*
** Copyright (C) 2007-2015 Erik de Castro Lopo <erikd@mega-nerd.com>
**
** 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 or version 3 of the
** License.
**
** 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, see <http://www.gnu.org/licenses/>.
*/

/*
**	Generate a spectrogram as a PNG file from a given sound file.
*/

/*
**	Todo:
**      - Decouple height of image from FFT length. FFT length should be
*         greater than height and then interpolated to height.
**      - Make magnitude to colour mapper allow abitrary scaling (ie cmdline
**        arg).
**      - Better cmdline arg parsing and flexibility.
**      - Add option to do log frequency scale.
*/

#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>
#include <limits.h>
#include <assert.h>

#include <cairo.h>
#include <fftw3.h>

#include <sndfile.h>

#include "window.h"
#include "common.h"

#define TICK_LEN			6
#define	BORDER_LINE_WIDTH	1.8

#define	TITLE_FONT_SIZE		20.0
#define	NORMAL_FONT_SIZE	12.0

#define	LEFT_BORDER			70.0
#define	TOP_BORDER			30.0
#define	RIGHT_BORDER		75.0
#define	BOTTOM_BORDER		40.0

#define	SPEC_FLOOR_DB		-180.0


typedef struct
{	const char *sndfilepath, *pngfilepath, *filename ;
	int width, height ;
	bool border, log_freq, gray_scale ;
	enum WINDOW_FUNCTION window_function ;
	double spec_floor_db ;
} RENDER ;

/* Interpolate between samples to eliminate the
 * 1-sample-period jitter due to time quantization?
 */
static bool precise_timing = 0;

typedef struct
{	int left, top, width, height ;
} RECT ;

static void
get_colour_map_value (float value, double spec_floor_db, unsigned char colour [3], bool gray_scale)
{	static unsigned char map [][3] =
	{	/* These values were originally calculated for a dynamic range of 180dB. */
		{	255,	255,	255	},	/* -0dB */
		{	240,	254,	216	},	/* -10dB */
		{	242,	251,	185	},	/* -20dB */
		{	253,	245,	143	},	/* -30dB */
		{	253,	200,	102	},	/* -40dB */
		{	252,	144,	66	},	/* -50dB */
		{	252,	75,		32	},	/* -60dB */
		{	237,	28,		41	},	/* -70dB */
		{	214,	3,		64	},	/* -80dB */
		{	183,	3,		101	},	/* -90dB */
		{	157,	3,		122	},	/* -100dB */
		{	122,	3,		126	},	/* -110dB */
		{	80,		2,		110	},	/* -120dB */
		{	45,		2,		89	},	/* -130dB */
		{	19,		2,		70	},	/* -140dB */
		{	1,		3,		53	},	/* -150dB */
		{	1,		3,		37	},	/* -160dB */
		{	1,		2,		19	},	/* -170dB */
		{	0,		0,		0	},	/* -180dB */
	} ;

	float rem ;
	int indx ;

	if (gray_scale)
	{	/* "value" is a negative value in decibels.
		 * black (0,0,0) is for <= -180.0,
		 * white (255,255,255) is for >= 0.0 leaving 254 values
		 * which should cover the intervening range evenly.
		 * (value/spec_floor_db) is between 0.0 and 1.0 (not inclusive)
		 * because both value and spec_floor_db are negative.
		 * (v/s) * 254.0 is from 0.000001 to 253.9999999 and
		 * we want 1 to 254, so ceil((v/s) * 254) gives us that,
		 * converted to 254 to 1 by subtracting from 255. */
		int c ; /* The pixel value */

		if (value <= spec_floor_db)
			c = 0 ;
		else
		{	c = 255 - lrintf (ceil ((value / spec_floor_db) * 254.0)) ;

			if (c < 1 || c > 254)	/* Sanity check */
			{	printf ("\nError : gray value is %d\n\n", c) ;
				exit (1) ;
				} ;

			} ;
		colour [0] = colour [1] = colour [2] = c ;
		return ;
		} ;

	if (value >= 0.0)
	{	colour [0] = colour [1] = colour [2] = 255 ;
		return ;
		} ;

	value = fabs (value * (-180.0 / spec_floor_db) * 0.1) ;

	indx = lrintf (floor (value)) ;

	if (indx < 0)
	{	printf ("\nError : colour map array index is %d\n\n", indx) ;
		exit (1) ;
		} ;

	if (indx >= ARRAY_LEN (map))
	{	colour [0] = colour [1] = colour [2] = 0 ;
		return ;
		} ;

	rem = fmod (value, 1.0) ;

	colour [0] = lrintf ((1.0 - rem) * map [indx][0] + rem * map [indx + 1][0]) ;
	colour [1] = lrintf ((1.0 - rem) * map [indx][1] + rem * map [indx + 1][1]) ;
	colour [2] = lrintf ((1.0 - rem) * map [indx][2] + rem * map [indx + 1][2]) ;

	return ;
} /* get_colour_map_value */


static void
read_mono_audio (SNDFILE * file, sf_count_t filelen, double * data, int datalen, int indx, int total)
{
	sf_count_t start ;
	double error ;	 /* How much "start" is wrong by
			  * (start == exact_value_of_start + error) */
	sf_count_t k;

	memset (data, 0, datalen * sizeof (data [0])) ;

	start = (indx * filelen) / total - datalen / 2 ;

	error = start - (((double)indx * filelen) / total - datalen / 2.0) ;
	assert (-1.0 < error && error <= 0.0) ;

	if (start >= 0)
		sf_seek (file, start, SEEK_SET) ;
	else
	{	start = -start ;
		sf_seek (file, 0, SEEK_SET) ;
		data += start ;
		datalen -= start ;
		} ;

	sfx_mix_mono_read_double (file, data, datalen+(precise_timing?1:0)) ;

	/* The starting time is truncated to a sample boundary and is slighty
	 * earlier that the exact starting time.
	 * We compensate by reading in one more sample than this and
	 * interpolating between each sample and the following one.
	 */
	if ( precise_timing && error != 0.0 )
	{
		/* Proportion to take of this sample and of the following one.
		 * this + next == 1.0
		 */
		double this = 1.0 + error, next = -error ;

		for (k=0; k < datalen; k++)
			data[k] = data[k]*this + data[k+1]*next;
	}

	return ;
} /* read_mono_audio */

static void
apply_window (double * data, int datalen, enum WINDOW_FUNCTION window_function)
{	static double *window = NULL ;
	static int window_len = 0 ;
	int k ;

	if (window_len != datalen)
	{	window = realloc (window, datalen * sizeof (double)) ;
		if (window == NULL)
		{	printf ("%s : Not enough memory.\n", __func__) ;
			exit (1) ;
			} ;
		window_len = datalen ;

		switch (window_function)
		{	case KAISER :
				calc_kaiser_window (window, datalen, 20.0) ;
				break ;
			case NUTTALL:
				calc_nuttall_window (window, datalen) ;
				break ;
			case HANN :
				calc_hann_window (window, datalen) ;
				break ;
			default :
				printf ("Internal error: Unknown window_function.\n") ;
				exit (1) ;
			} ;
		} ;

	for (k = 0 ; k < datalen ; k++)
		data [k] *= window [k] ;

	return ;
} /* apply_window */

static void
render_spectrogram (cairo_surface_t * surface, double spec_floor_db, float **mag2d, double maxval, double left, double top, double width, double height, bool gray_scale)
{
	unsigned char colour [3] = { 0, 0, 0 } ;
	unsigned char *data ;
	double linear_spec_floor ;
	int w, h, stride ;

	stride = cairo_image_surface_get_stride (surface) ;

	data = cairo_image_surface_get_data (surface) ;
	memset (data, 0, stride * cairo_image_surface_get_height (surface)) ;

	linear_spec_floor = pow (10.0, spec_floor_db / 20.0) ;

	for (w = 0 ; w < width ; w ++)
		for (h = 0 ; h < height ; h++)
		{	int x, y ;

			mag2d [w][h] = mag2d [w][h] / maxval ;
			mag2d [w][h] = (mag2d [w][h] < linear_spec_floor) ? spec_floor_db : 20.0 * log10 (mag2d [w][h]) ;

			get_colour_map_value (mag2d [w][h], spec_floor_db, colour, gray_scale) ;

			y = height + top - 1 - h ;
			x = (w + left) * 4 ;
			data [y * stride + x + 0] = colour [2] ;
			data [y * stride + x + 1] = colour [1] ;
			data [y * stride + x + 2] = colour [0] ;
			data [y * stride + x + 3] = 0 ;
			} ;

	cairo_surface_mark_dirty (surface) ;
} /* render_spectrogram */

static void
render_heat_map (cairo_surface_t * surface, double magfloor, const RECT *r, bool gray_scale)
{
	unsigned char colour [3], *data ;
	int w, h, stride ;

	stride = cairo_image_surface_get_stride (surface) ;
	data = cairo_image_surface_get_data (surface) ;

	for (h = 0 ; h < r->height ; h++)
	{	get_colour_map_value (magfloor * (r->height - h) / (r->height + 1), magfloor, colour, gray_scale) ;

		for (w = 0 ; w < r->width ; w ++)
		{	int x, y ;

			x = (w + r->left) * 4 ;
			y = r->height + r->top - 1 - h ;

			data [y * stride + x + 0] = colour [2] ;
			data [y * stride + x + 1] = colour [1] ;
			data [y * stride + x + 2] = colour [0] ;
			data [y * stride + x + 3] = 0 ;
			} ;
		} ;

	cairo_surface_mark_dirty (surface) ;
} /* render_heat_map */

static inline void
x_line (cairo_t * cr, double x, double y, double len)
{	cairo_move_to (cr, x, y) ;
	cairo_rel_line_to (cr, len, 0.0) ;
	cairo_stroke (cr) ;
} /* x_line */

static inline void
y_line (cairo_t * cr, double x, double y, double len)
{	cairo_move_to (cr, x, y) ;
	cairo_rel_line_to (cr, 0.0, len) ;
	cairo_stroke (cr) ;
} /* y_line */

typedef struct
{	double value [15] ;
	double distance [15] ;
} TICKS ;

static inline int
calculate_ticks (double max, double distance, TICKS * ticks)
{	const int div_array [] =
	{	10, 10, 8, 6, 8, 10, 6, 7, 8, 9, 10, 11, 12, 12, 7, 14, 8, 8, 9
		} ;

	double scale = 1.0, scale_max ;
	int k, leading, divisions ;

	if (max < 0)
	{	printf ("\nError in %s : max < 0\n\n", __func__) ;
		exit (1) ;
		} ;

	while (scale * max >= ARRAY_LEN (div_array))
		scale *= 0.1 ;

	while (scale * max < 1.0)
		scale *= 10.0 ;

	leading = lround (scale * max) ;
	divisions = div_array [leading % ARRAY_LEN (div_array)] ;

	/* Scale max down. */
	scale_max = leading / scale ;
	scale = scale_max / divisions ;

	if (divisions > ARRAY_LEN (ticks->value) - 1)
	{	printf ("Error : divisions (%d) > ARRAY_LEN (ticks->value) (%d)\n", divisions, ARRAY_LEN (ticks->value)) ;
		exit (1) ;
		} ;

	for (k = 0 ; k <= divisions ; k++)
	{	ticks->value [k] = k * scale ;
		ticks->distance [k] = distance * ticks->value [k] / max ;
		} ;

	return divisions + 1 ;
} /* calculate_ticks */

static void
str_print_value (char * text, int text_len, double value)
{
	if (fabs (value) < 1e-10)
		snprintf (text, text_len, "0") ;
	else if (fabs (value) >= 10.0)
		snprintf (text, text_len, "%1.0f", value) ;
	else if (fabs (value) >= 1.0)
		snprintf (text, text_len, "%3.1f", value) ;
	else
		snprintf (text, text_len, "%4.2f", value) ;

	return ;
} /* str_print_value */

static void
render_spect_border (cairo_surface_t * surface, const char * filename, double left, double width, double seconds, double top, double height, double max_freq)
{
	char text [512] ;
	cairo_t * cr ;
	cairo_text_extents_t extents ;
	cairo_matrix_t matrix ;

	TICKS ticks ;
	int k, tick_count ;

	cr = cairo_create (surface) ;

	cairo_set_source_rgb (cr, 1.0, 1.0, 1.0) ;
	cairo_set_line_width (cr, BORDER_LINE_WIDTH) ;

	/* Print title. */
	cairo_select_font_face (cr, font_family, CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_NORMAL) ;
	cairo_set_font_size (cr, 1.0 * TITLE_FONT_SIZE) ;

	snprintf (text, sizeof (text), "Spectrogram: %s", filename) ;
	cairo_text_extents (cr, text, &extents) ;
	cairo_move_to (cr, left + 2, top - extents.height / 2) ;
	cairo_show_text (cr, text) ;

	/* Print labels. */
	cairo_select_font_face (cr, font_family, CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_NORMAL) ;
	cairo_set_font_size (cr, 1.0 * NORMAL_FONT_SIZE) ;

	/* Border around actual spectrogram. */
	cairo_rectangle (cr, left, top, width, height) ;

	tick_count = calculate_ticks (seconds, width, &ticks) ;
	for (k = 0 ; k < tick_count ; k++)
	{	/* Don't draw the tick if its further left than the right border. */
		if (left + ticks.distance [k] > width)
			continue ;

		y_line (cr, left + ticks.distance [k], top + height, TICK_LEN) ;
		if (k % 2 == 1)
			continue ;
		str_print_value (text, sizeof (text), ticks.value [k]) ;
		cairo_text_extents (cr, text, &extents) ;
		cairo_move_to (cr, left + ticks.distance [k] - extents.width / 2, top + height + 8 + extents.height) ;
		cairo_show_text (cr, text) ;
		} ;

	tick_count = calculate_ticks (max_freq, height, &ticks) ;
	for (k = 0 ; k < tick_count ; k++)
	{	x_line (cr, left + width, top + height - ticks.distance [k], TICK_LEN) ;
		if (k % 2 == 1)
			continue ;
		str_print_value (text, sizeof (text), ticks.value [k]) ;
		cairo_text_extents (cr, text, &extents) ;
		cairo_move_to (cr, left + width + 12, top + height - ticks.distance [k] + extents.height / 4.5) ;
		cairo_show_text (cr, text) ;
		} ;

	cairo_set_font_size (cr, 1.0 * NORMAL_FONT_SIZE) ;

	/* Label X axis. */
	snprintf (text, sizeof (text), "Time (secs)") ;
	cairo_text_extents (cr, text, &extents) ;
	cairo_move_to (cr, left + (width - extents.width) / 2, cairo_image_surface_get_height (surface) - 8) ;
	cairo_show_text (cr, text) ;

	/* Label Y axis (rotated). */
	snprintf (text, sizeof (text), "Frequency (Hz)") ;
	cairo_text_extents (cr, text, &extents) ;

	cairo_get_font_matrix (cr, &matrix) ;
	cairo_matrix_rotate (&matrix, -0.5 * M_PI) ;
	cairo_set_font_matrix (cr, &matrix) ;

	cairo_move_to (cr, cairo_image_surface_get_width (surface) - 12, top + (height + extents.width) / 2) ;
	cairo_show_text (cr, text) ;

	cairo_destroy (cr) ;
} /* render_spect_border */

static void
render_heat_border (cairo_surface_t * surface, double magfloor, const RECT *r)
{
	const char *decibels = "dB" ;
	char text [512] ;
	cairo_t * cr ;
	cairo_text_extents_t extents ;
	TICKS ticks ;
	int k, tick_count ;

	cr = cairo_create (surface) ;

	cairo_set_source_rgb (cr, 1.0, 1.0, 1.0) ;
	cairo_set_line_width (cr, BORDER_LINE_WIDTH) ;

	/* Border around actual spectrogram. */
	cairo_rectangle (cr, r->left, r->top, r->width, r->height) ;
	cairo_stroke (cr) ;

	cairo_select_font_face (cr, font_family, CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_NORMAL) ;
	cairo_set_font_size (cr, 1.0 * NORMAL_FONT_SIZE) ;

	cairo_text_extents (cr, decibels, &extents) ;
	cairo_move_to (cr, r->left + (r->width - extents.width) / 2, r->top - 5) ;
	cairo_show_text (cr, decibels) ;

	tick_count = calculate_ticks (fabs (magfloor), r->height, &ticks) ;
	for (k = 0 ; k < tick_count ; k++)
	{	x_line (cr, r->left + r->width, r->top + ticks.distance [k], TICK_LEN) ;
		if (k % 2 == 1)
			continue ;

		str_print_value (text, sizeof (text), -1.0 * ticks.value [k]) ;
		cairo_text_extents (cr, text, &extents) ;
		cairo_move_to (cr, r->left + r->width + 2 * TICK_LEN, r->top + ticks.distance [k] + extents.height / 4.5) ;
		cairo_show_text (cr, text) ;
		} ;

	cairo_destroy (cr) ;
} /* render_heat_border */

static void
interp_spec (float * mag, int maglen, const double *spec, int speclen)
{
	int k, lastspec = 0 ;

	mag [0] = spec [0] ;

	for (k = 1 ; k < maglen ; k++)
	{	double sum = 0.0 ;
		int count = 0 ;

		do
		{	sum += spec [lastspec] ;
			lastspec ++ ;
			count ++ ;
			}
		while (lastspec <= ceil ((k * speclen) / maglen)) ;

		mag [k] = sum / count ;
		} ;

	return ;
} /* interp_spec */

/* Pick the best FFT length good for FFTW?
**
** We use fftw_plan_r2r_1d() for which the documantation
** http://fftw.org/fftw3_doc/Real_002dto_002dReal-Transforms.html says:
**
** "FFTW is generally best at handling sizes of the form
** 2^a 3^b 5^c 7^d 11^e 13^f
** where e+f is either 0 or 1, and the other exponents are arbitrary."
*/

/* Helper function: does N have only 2, 3, 5 and 7 as its factors? */
static bool
is_2357 (int n)
{
	/* Just eliminate all factors os 2, 3, 5 and 7 and see if 1 remains */
	while (n % 2 == 0) n /= 2 ;
	while (n % 3 == 0) n /= 3 ;
	while (n % 5 == 0) n /= 5 ;
	while (n % 7 == 0) n /= 7 ;
	return (n == 1) ;
}

/* Helper function: is N a "fast" value for the FFT size? */
static bool
is_good_speclen (int n)
{
	/* It wants n, 11*n, 13*n but not (11*13*n)
	** where n only has as factors 2, 3, 5 and 7
	*/
	if (n % (11 * 13) == 0) return 0 ; /* No good */

	return is_2357 (n)	|| ((n % 11 == 0) && is_2357 (n / 11))
						|| ((n % 13 == 0) && is_2357 (n / 13)) ;
}

static void
render_to_surface (const RENDER * render, SNDFILE *infile, int samplerate, sf_count_t filelen, cairo_surface_t * surface)
{
	double * time_domain = NULL ;
	double * freq_domain = NULL ;
	double * single_mag_spec = NULL ;
	float ** mag_spec = NULL ; // Indexed by [w][h]

	fftw_plan plan ;
	double max_mag = 0.0 ;
	int width, height, w, speclen ;

	if (render->border)
	{	width = lrint (cairo_image_surface_get_width (surface) - LEFT_BORDER - RIGHT_BORDER) ;
		height = lrint (cairo_image_surface_get_height (surface) - TOP_BORDER - BOTTOM_BORDER) ;
		}
	else
	{	width = render->width ;
		height = render->height ;
		}

	if (width < 1)
	{	printf ("Error : 'width' parameter must be >= %d\n",
			render->border ? (int) (LEFT_BORDER + RIGHT_BORDER) + 1 : 1) ;
		exit (1) ;
		} ;

	if (height < 1)
	{	printf ("Error : 'height' parameter must be >= %d\n",
			render->border ? (int) (TOP_BORDER + BOTTOM_BORDER) + 1 : 1) ;
		exit (1) ;
		} ;

	/*
	** Choose a speclen value that is long enough to represent frequencies
	** down to 20Hz.
	*/
	speclen = height * (samplerate / 20 / height + 1) ;

	/* Find the nearest fast value for the FFT size. */
	{	int d ;	/* difference */

		for (d = 0 ; /* Will terminate */ ; d++)
		{	/* Logarithmically, the integer above is closer than
			** the integer below, so prefer it to the one below.
			*/
			if (is_good_speclen (speclen + d))
			{	speclen += d ;
				break ;
				}
			/* FFT length must also be >= the output height,
			** otherwise repeated pixel rows occur in the output.
			*/
			if (speclen - d >= height && is_good_speclen (speclen - d))
			{	speclen -= d ;
				break ;
				}
			}
		}

	/* + 1 because with --precise, read_mono_audio reads the following
	 * sample and interpolates to get the starting moment exact. */
	time_domain		= calloc ((2 * speclen + 1), sizeof (double)) ;
	freq_domain		= calloc (2 * speclen, sizeof (double)) ;
	single_mag_spec = calloc (speclen, sizeof (double)) ;
	mag_spec		= calloc (width, sizeof (float *)) ;
	if (time_domain == NULL || freq_domain == NULL || single_mag_spec == NULL || mag_spec == NULL)
	{	printf ("%s : Not enough memory.\n", __func__) ;
		exit (1) ;
		} ;
	for (w = 0 ; w < width ; w++)
	{	if ((mag_spec [w] = calloc (height, sizeof (float))) == NULL)
		{	printf ("%s : Not enough memory.\n", __func__) ;
			exit (1) ;
			} ;
		} ;

	plan = fftw_plan_r2r_1d (2 * speclen, time_domain, freq_domain, FFTW_R2HC, FFTW_MEASURE | FFTW_PRESERVE_INPUT) ;
	if (plan == NULL)
	{	printf ("%s : line %d : create plan failed.\n", __FILE__, __LINE__) ;
		exit (1) ;
		} ;

	for (w = 0 ; w < width ; w++)
	{	double single_max ;

		read_mono_audio (infile, filelen, time_domain, 2 * speclen, w, width) ;

		apply_window (time_domain, 2 * speclen, render->window_function) ;

		fftw_execute (plan) ;

		single_max = calc_magnitude (freq_domain, 2 * speclen, single_mag_spec) ;
		max_mag = MAX (max_mag, single_max) ;

		interp_spec (mag_spec [w], height, single_mag_spec, speclen) ;
		} ;

	fftw_destroy_plan (plan) ;
	free (time_domain) ;
	free (freq_domain) ;
	free (single_mag_spec) ;

	if (render->border)
	{	RECT heat_rect ;

		heat_rect.left = 12 ;
		heat_rect.top = TOP_BORDER + TOP_BORDER / 2 ;
		heat_rect.width = 12 ;
		heat_rect.height = height - TOP_BORDER / 2 ;

		render_spectrogram (surface, render->spec_floor_db, mag_spec, max_mag, LEFT_BORDER, TOP_BORDER, width, height, render->gray_scale) ;

		render_heat_map (surface, render->spec_floor_db, &heat_rect, render->gray_scale) ;

		render_spect_border (surface, render->filename, LEFT_BORDER, width, filelen / (1.0 * samplerate), TOP_BORDER, height, 0.5 * samplerate) ;
		render_heat_border (surface, render->spec_floor_db, &heat_rect) ;
		}
	else
		render_spectrogram (surface, render->spec_floor_db, mag_spec, max_mag, 0, 0, width, height, render->gray_scale) ;

	for (w = 0 ; w < width ; w++)
		free (mag_spec [w]) ;
	free (mag_spec) ;

	return ;
} /* render_to_surface */

static void
render_cairo_surface (const RENDER * render, SNDFILE *infile, int samplerate, sf_count_t filelen)
{
	cairo_surface_t * surface = NULL ;
	cairo_status_t status ;

	/*
	**	CAIRO_FORMAT_RGB24 	 each pixel is a 32-bit quantity, with the upper 8 bits
	**	unused. Red, Green, and Blue are stored in the remaining 24 bits in that order.
	*/

	surface = cairo_image_surface_create (CAIRO_FORMAT_RGB24, render->width, render->height) ;
	if (surface == NULL || cairo_surface_status (surface) != CAIRO_STATUS_SUCCESS)
	{	status = cairo_surface_status (surface) ;
		printf ("Error while creating surface : %s\n", cairo_status_to_string (status)) ;
		exit (1) ;
		} ;

	cairo_surface_flush (surface) ;

	render_to_surface (render, infile, samplerate, filelen, surface) ;

	status = cairo_surface_write_to_png (surface, render->pngfilepath) ;
	if (status != CAIRO_STATUS_SUCCESS)
	{	printf ("Error while creating PNG file : %s\n", cairo_status_to_string (status)) ;
		exit (1) ;
		} ;

	cairo_surface_destroy (surface) ;

	return ;
} /* render_cairo_surface */

static void
render_sndfile (const RENDER * render)
{
	SNDFILE *infile ;
	SF_INFO info ;

	if (render->log_freq)
	{	printf ("Error : --log-freq option not working yet.\n\n") ;
		exit (1) ;
		} ;

	memset (&info, 0, sizeof (info)) ;

	infile = sf_open (render->sndfilepath, SFM_READ, &info) ;
	if (infile == NULL)
	{	printf ("Error : failed to open file '%s' : \n%s\n", render->sndfilepath, sf_strerror (NULL)) ;
		exit (1) ;
		} ;

	render_cairo_surface (render, infile, info.samplerate, info.frames) ;

	sf_close (infile) ;

	return ;
} /* render_sndfile */

static void
usage_exit (const char * argv0, int error)
{
	const char * progname ;

	progname = strrchr (argv0, '/') ;
	progname = (progname == NULL) ? argv0 : progname + 1 ;

	printf ("\nUsage :\n\n    %s [options] <sound file> <img width> <img height> <png name>\n\n", progname) ;

	puts (
		"    Create a spectrogram as a PNG file from a given sound file. The\n"
		"    spectrogram image will be of the given width and height.\n"
		) ;

	puts (
		"    Options:\n"
		"        --dyn-range=<number>   : Dynamic range (ie 100 for 100dB range)\n"
		"        --no-border            : Drop the border, scales, heat map and title\n"
		/*-"        --log-freq             : Use a logarithmic frquency scale\n" -*/
		"        --gray-scale           : Output gray pixels instead of a heat map\n"
		"        --kaiser               : Use a Kaiser window function (the default)\n"
		"        --nuttall              : Use a Nuttall window function\n"
		"        --hann                 : Use a Hann window function\n"
		"        --precise              : Interpolate between samples for precise timing\n"
		) ;

	exit (error) ;
} /* usage_exit */

int
main (int argc, char * argv [])
{	RENDER render =
	{	NULL, NULL, NULL,
		0, 0,
		true, false, false,
		KAISER,
		SPEC_FLOOR_DB
		} ;
	int k ;

	if (argc < 5)
		usage_exit (argv [0], 0) ;

	for (k = 1 ; k < argc - 4 ; k++)
	{	double fval ;

		if (sscanf (argv [k], "--dyn-range=%lf", &fval) == 1)
		{	render.spec_floor_db = -1.0 * fabs (fval) ;
			continue ;
			}

		if (strcmp (argv [k], "--no-border") == 0)
		{	render.border = false ;
			continue ;
			} ;

		if (strcmp (argv [k], "--log-freq") == 0)
		{	render.log_freq = true ;
			continue ;
			} ;

		if (strcmp (argv [k], "--gray-scale") == 0)
		{	render.gray_scale = true ;
			continue ;
			} ;

		if (strcmp (argv [k], "--kaiser") == 0)
		{	render.window_function = KAISER ;
			continue ;
			} ;

		if (strcmp (argv [k], "--nuttall") == 0)
		{	render.window_function = NUTTALL ;
			continue ;
			} ;

		if (strcmp (argv [k], "--hann") == 0)
		{	render.window_function = HANN ;
			continue ;
			} ;

		if (strcmp (argv [k], "--precise") == 0)
		{	precise_timing = 1;
			continue ;
			} ;

		printf ("\nError : Bad command line argument '%s'\n", argv [k]) ;
		usage_exit (argv [0], 1) ;
		} ;

	render.sndfilepath = argv [k] ;
	render.width = parse_int_or_die (argv [k + 1], "width") ;
	render.height = parse_int_or_die (argv [k + 2], "height") ;
	render.pngfilepath = argv [k + 3] ;

	render.filename = strrchr (render.sndfilepath, '/') ;
	render.filename = (render.filename != NULL) ? render.filename + 1 : render.sndfilepath ;

	render_sndfile (&render) ;

	return 0 ;
} /* main */