//cheapchop.c
//Measures DC voltages using the CheapChop circuit and a soundcard

//	(C) Marko Cebokli S57UUU	aug 2007
//	GNU/GPL licence:  www.gnu.org/licences/licences.html#GPL

//compile:	gcc -lm cheapchop.c -o cheapchop

#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <linux/soundcard.h>
#include <stdlib.h>

double PI=3.14159265358979;


//--------------------------------------------------------------
//open sound card
//bits = 8 or 16
//chans = 1 or 2
//srate = what is possible depends on the card...
void sc_open(int *fd, int bits, int chans, int srate)
{
int arg,status;

*fd = open("/dev/dsp", O_RDWR);
if (fd < 0)
	{
	perror("open of /dev/dsp failed");
	exit(EXIT_FAILURE);
	}

/* set sampling parameters */
arg = bits;      /* sample size */

status = ioctl(*fd, SOUND_PCM_WRITE_BITS, &arg);
if (status == -1)
	{
	perror("SOUND_PCM_WRITE_BITS ioctl failed");
	exit(EXIT_FAILURE);
	}
if (arg != bits)
	{
	perror("unable to set sample size");
	printf("\nDesired size = %d bits,   actual %d",bits,arg);
	exit(EXIT_FAILURE);
	}

arg = chans;  /* mono or stereo */
status = ioctl(*fd, SOUND_PCM_WRITE_CHANNELS, &arg);
if (status == -1)
	{
	perror("SOUND_PCM_WRITE_CHANNELS ioctl failed");
	exit(EXIT_FAILURE);
	}
if (arg != chans)
	{
	perror("unable to set number of channels");
	printf("\nDesired %d chans,   actual %d",chans,arg);
	exit(EXIT_FAILURE);
	}

arg = srate;      /* sampling rate */
status = ioctl(*fd, SOUND_PCM_WRITE_RATE, &arg);
if (status == -1)
	{
	perror("SOUND_PCM_WRITE_RATE ioctl failed");
	exit(EXIT_FAILURE);
	}
if (arg != srate)
	{
	perror("unable to set sample rate");
	printf("\nDesired srate = %d Hz,   actual %d",srate,arg);
	exit(EXIT_FAILURE);
	}

}

//--------------------------------------------------------------
//make a simple FIR bandpass filter
//nf	length of filter
//fc	center freq of filter
//fs	sampling frequency
//fr[], fi[]	real and imag coefficients
//g	returns gain
void makefilter(int nf, float fc, float fs, float fr[], float fi[], float *g)
{
int i;
float w,n1;

*g=0.0;
n1=(float)nf/2.0-0.5;	//simetric carrier
for (i=0;i<nf;i++)
	{
//	w=0.5-0.5*cos((float)i/(float)nf*2.0*PI);   //hanning window
	w=0.5-0.5*cos(((float)i+0.5)/(float)nf*2.0*PI);  //symmetric
	fr[i]=w*cosf(((float)i-n1)*PI*fc/fs*2.0);
	fi[i]=w*sinf(((float)i-n1)*PI*fc/fs*2.0);
	*g=*g+w;
	}
*g=*g/(float)nf;
}

//---------------------------------------------------------------
//bandpass detector
//nf = length of the filter (see makefilter function above)
//n = where in the array data[] the samples are taken
//g = gain of the filter
float detektor(int nf, int n, float data[], float fr[], float fi[], float g)
{
int i;
float dr,di,dc,dcr,dci;

dcr=0.0;dci=0.0;
for (i=0;i<nf;i++)
	{
	dcr=dcr+fr[i];
	dci=dci+fi[i];
	}
dr=0.0;di=0.0;dc=0.0;
for (i=0;i<nf;i++)
	{
	dr=dr+data[i+n]*fr[i];
	di=di+data[i+n]*fi[i];
	dc=dc+data[i+n];
	}

dc=dc/(float)nf;
dr=dr-dcr*dc;
di=di-dci*dc;
return sqrtf(dr*dr+di*di)/(float)nf/g*2.0;
}

//**********************************************************
int main()
{
short buf[200000];
float fr[64],fi[64];
float data[100000];
float d1,sat1,sat2,fd,g;
float avg1,avg2;
int i,j,k,nf,navg,length,status;
int scfd;

printf("\n");

fd=5000.0;	//frequency, to which the detector is tuned
nf=64;		//filter length
navg=64;	//amount of averaging
length=nf*navg;

makefilter(nf,fd,48000.0,fr,fi,&g);

sc_open(&scfd,16,2,48000);

while (1) 
	{ /* loop until Control-C */

	status = read(scfd, buf, 4*length);	//2chans * 2bytes/samp
	if (status != (4*length))
		{
		perror("read wrong number of bytes");
		}

	sat1=0;sat2=0;
	for (i=0;i<length;i++)
		{
		if ((buf[2*i]>32760)||(buf[2*i]<-32760))
			sat1++;
		if ((buf[2*i+1]>32760)||(buf[2*i+1]<-32760))
			sat2++;
		data[i]=(float)buf[2*i];
		}
	avg1=0.0;j=0;
	for (i=0;i<navg;i++)
		{
		d1=detektor(nf,j,data,fr,fi,g);
		avg1=avg1+d1;
		j=j+nf;
		}
	avg1=avg1/(float)navg;
	for (i=0;i<length;i++)
		{
		data[i]=(float)buf[2*i+1];
		}
	avg2=0.0;j=0;
	for (i=0;i<navg;i++)
		{
		d1=detektor(nf,j,data,fr,fi,g);
		avg2=avg2+d1;
		j=j+nf;
		}
	avg2=avg2/(float)navg;
//	avg1=avg1/5.062;	//line inp min gain = 5.062/mV
//	avg2=avg2/5.062;
	if ((sat1==0)&&(sat2==0))
		printf("\nCH1 = %7.1f     CH2 = %7.1f",avg1,avg2);
	if ((sat1!=0)&&(sat2==0))
		printf("\nCH1 = %7.1f**   CH2 = %7.1f",avg1,avg2);
	if ((sat1==0)&&(sat2!=0))
	printf("\nCH1 = %7.1f   CH2 = %7.1f**",avg1,avg2);
	if ((sat1!=0)&&(sat2!=0))
	printf("\nCH1 = %7.1f**   CH2 = %7.1f**",avg1,avg2);
}

printf("\n\n");		//never reaches this really :-)
return 0;
}