/*
 * sig_detect: how small a signal can be detected by simply sampling
 * and averaging?
 *
 * created, 2003-01-30, bsy@cs.ucsd.edu.
 *
 * -s specifies the signal strength relative to the noise, in dB.
 * -n specifies the number of samples to average.
 *
 * since this is a simulation and we know a priori whether the signal
 * is present or not, so we can perfectly assign the sample values to
 * buckets -- so we just do the calculations per bucket separately.
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>

int	debug = 0;

double	sample_mean(double	sig_dB,
		    int		nsamples)
{
	double	sum;
	double	noise, sig, sample;
	int	i;

	/*
	 * this code assumes no overflow, no loss of precision (!) in
	 * the double precision floating point arithmetic for
	 * accumulating the sum
	 */
	for (i = 0, sum = 0.0; i < nsamples; i++) {
		noise = drand48();
		if (sig_dB != -HUGE_VAL)
			sig = drand48() * pow(10.0,sig_dB/10.0);
		else
			sig = 0.0;
		sample = sig + noise;
		if (debug) {
			printf("noise:  %-15.13g\n",noise);
			printf("signal: %-15.13g\n",sig);
			printf("sample: %-15.13g\n",sample);
		}
		sum += sample;
	}
	return sum / nsamples;
}

void	usage(void)
{
	fprintf(stderr,"Usage: sig_detect [-dvV] [-i iterations] [-n nsamples] [-s sig_strengthFloat_dB]\n");
}

int	main(int	ac,
	     char	**av)
{
	double	sig_strength = -30.0;
	int	nsamples = 1000000;
	int	opt;
	int	verbose = 0;
	int	iterations = 1;
	double	mean_sig, mean_no_sig;
	int	detected, ndetected;

	while (EOF != (opt = getopt(ac,av,"di:n:s:vV"))) switch (opt) {
	case 'd':
		++debug;
		break;
	case 'i':
		iterations = atoi(optarg);
		break;
	case 'n':
		nsamples = atoi(optarg);
		break;
	case 's':
		sig_strength = atof(optarg);
		break;
	case 'v':
		++verbose;
		break;
	case 'V':
		--verbose;
		break;
	default:
		usage();
		exit(1);
	}
		
	if (sig_strength >= 0.0) {
		fprintf(stderr,"sig_detect: Signal strength must be less than 0 dB\n");
		usage();
		exit(1);
	}

	srand48(time((time_t *)0) ^ getpid());

	if (verbose) {
		printf("CSE127, bsy\n\n");
		printf("Signal Detection module\n");
		printf("------ --------- ------\n\n");
		printf("Signal strength:   %g dB\n",sig_strength);
		printf("Number of samples: %d\n\n",nsamples);
	}
	ndetected = 0;
	while (--iterations >= 0) {
		mean_no_sig = sample_mean(-HUGE_VAL,nsamples);
		printf("no signal average: %-15.013g\n",mean_no_sig);
		mean_sig = sample_mean(sig_strength,nsamples);
		detected = mean_sig > mean_no_sig;
		printf("   signal average: %-15.013g%s\n",mean_sig,detected?" *":"");
		ndetected += detected;
	}
	printf("\n%d/%d detection rate\n",ndetected,iterations);
	return 0;
}