//SIDI1FFT.C   single bit sampling on LPT port, then display FFT
//INTEGER version

//             Marko Cebokli  feb 2004   GNU/GPL licence

//             www.gnu.org/licences/licences.html#GPL

//      version 1 - software sample timing


//  LPT port connection:

//				   pin 11 channel 1 I             bit 7
//                 pin 10 channel 1 Q             bit 6
//                 pin 12 channel 2 I             bit 5
//                 pin 13 channel 2 Q             bit 4
//                 pin 18..25  ground


#include <stdio.h>
#include <math.h>
#include <dos.h>     //inportb
#include <time.h>
#include <graphics.h>
#include <conio.h>   //getch() kbhit()
#include <stdlib.h>  //rand()

double PI=3.14159265358979;

unsigned int sig[32767];
int wr[32767],wi[32767];
int bi[32767];
int xr[32767],xi[32767];
float graf[16383];

//--------------------------------------------------------------

void fft(int n, int xr[], int xi[], int wr[], int wi[], int bi[])
{
int b,dp,a,t,p;
long xwr,xwi;

for (b=0;b<n;b++)
  if (bi[b]>b)
	{xwr=xr[b];xwi=xi[b];
	xr[b]=xr[bi[b]];xi[b]=xi[bi[b]];
	xr[bi[b]]=xwr;xi[bi[b]]=xwi;}

b=1;dp=n>>1;
do
  {a=0;
  do
	{t=a+b;p=0;
	do
	  {xwr=(long)xr[a+b]*wr[p]/65536-(long)xi[a+b]*wi[p]/65536;
	  xwi=(long)xr[a+b]*wi[p]/65536+(long)xi[a+b]*wr[p]/65536;
	  xr[a+b]=xr[a]/2-xwr;xi[a+b]=xi[a]/2-xwi;
	  xr[a]=xr[a]/2+xwr;xi[a]=xi[a]/2+xwi;
	  p=p+dp;a=a++;}
	while (a<t);
	a=a+b;}
  while (a<n);
  dp=dp>>1;b=b<<1;}
while (b<n);
}


//--------------------------------------------------------------------

void fftmojprip(int n,int wr[], int wi[], int bi[])
{
int i,j,k,l,m;
m=(int)(log((float)n)/log(1.9999999999999));
for (l=0;l<n;l++)    //bit reverse adrese
  {i=0;  j=l<<1;
  for (k=0;k<m;k++)
	i=(j=j>>1)%2+(i=i<<1);
  bi[l]=i;}

for (l=0;l<n;l++)    //sin/cos tabela
  {wr[l]=32767*cos((float)l/(float)n*2*PI);
  wi[l]=-32767*sin((float)l/(float)n*2*PI);}
}

//---------------------------------------------------------------------

void pisikompl(int n, int xr[], int xi[])
{
int i;
for (i=0;i<n;i++)
  printf("%8d %8d |",xr[i],xi[i]);
}

//-------------------------------------------------------------------

void siggen(int n, float h, int* xr, int* xi)
{
int i;
for (i=0;i<n;i++)
  {
//  xr[i]=100;xi[i]=0;    //DC
//  xr[i]=32767*cos((float)i/(float)n*h*2*PI);xi[i]=0;      //sine
  xr[i]=32767*((float)rand()/RAND_MAX-.5);xi[i]=0.;       //random
//  xr[i]=65535*((i/32)%2-.5);xi[i]=0.;                 //square
  }
}

//-----------------------------------------------------------

void vzorcenje(long n)     //this one does the sampling
{                    //array is written each time to equalize timing
int i;
long j;
unsigned char a;

for (j=0;j<n;j++)  //131072
  {
  i=(int)(j>>2);
  a=inportb(0x379)>>4;    //D25:   *pin11  pin10  pin12  pin13  (pin15)
  sig[i]=(sig[i]<<4)+a;
  }
}

//----------------------------------------------------------

void hitrovz(long n)		//fast sampling without bit concatenation
{
long j;

for (j=0;j<n;j++)
  {
  sig[j]=inportb(0x379);
  }
}

//----------------------------------------------------------

void preflintaj1(long n, int k) //podatki iz "hitrovz" v real/imag vhodne podatke FFT
{
long j;
int w;
for (j=0;j<n;j++)
  {
  w=32767;      //brez windowing, ce nimas koprocesorja
//  w=(int)(32767.0*sin((float)j/n*3.141592654));  //window (brez veze?)
  if (k==1)
	{
	if ((sig[j]&0x40)==0x40) xi[j]=w; else xi[j]=-w;   //Q1
	if ((sig[j]&0x80)==0x80) xr[j]=w; else xr[j]=-w;   //I1
	//xr[j]=0;
	}
  else
	{
	if ((sig[j]&0x10)==0x10) xi[j]=w; else xi[j]=-w;   //Q2
	if ((sig[j]&0x20)==0x20) xr[j]=-w; else xr[j]=w;   //I2
	//xr[j]=0;
	}
  }
}

//----------------------------------------------------------

float checkspeed()       //measure speed of computer (sample rate etc)
{
clock_t sta,mid,sto;
float srate,brate;

//printf("\n clktck=%8.2f \n",CLK_TCK);
sta=clock();
vzorcenje(131072);  //do the sampling
mid=clock();
hitrovz(32767);hitrovz(32767);hitrovz(32767);hitrovz(32767);
hitrovz(32767);hitrovz(32767);hitrovz(32767);hitrovz(32767);
hitrovz(32767);hitrovz(32767);hitrovz(32767);hitrovz(32767);
hitrovz(32767);hitrovz(32767);hitrovz(32767);hitrovz(32767);
hitrovz(32767);hitrovz(32767);hitrovz(32767);hitrovz(32767);
hitrovz(32767);hitrovz(32767);hitrovz(32767);hitrovz(32767);
hitrovz(32767);hitrovz(32767);hitrovz(32767);hitrovz(32767);
hitrovz(32767);hitrovz(32767);hitrovz(32767);hitrovz(32767);
sto=clock();
//printf("\n cas vzorcenja= %d",mid-sta);
printf("\n sampling time = %d",sto-mid);
srate=131072.0/(mid-sta)*CLK_TCK;
//printf("\n sampling rate= %8.0f /sec",srate);
srate=8*131072.0/(sto-mid)*CLK_TCK;
printf("\n sampling rate= %8.0f /sec",srate);
return srate;
}

//-----------------------------------------------------------------

void risigraf1(float graf[], long tock, float miny, float maxy)
{
int x,y;
long j;
//line(0,0,getmaxx(),getmaxy());


setviewport(50,11,562,227,1);
clearviewport();

setcolor(GREEN);
moveto(0,0);lineto(512,0);lineto(512,216);lineto(0,216);lineto(0,0);
for (y=36;y<216;y=y+36) line(0,y,512,y);
for (x=52;x<500;x=x+51) line(x,0,x,216);
setcolor(WHITE);

x=0;y=216-(graf[(int)(tock/2)]-miny)/(maxy-miny)*216;
moveto(x,y);
for (j=1;j<(int)(tock/2);j++)
  {
  x=(float)j/(float)tock*512;
  y=216-(graf[j+(int)(tock/2)]-miny)/(maxy-miny)*216;
  lineto(x,y);
  }

x=256;y=216-(graf[0]-miny)/(maxy-miny)*216;
moveto(x,y);
for (j=1;j<(int)(tock/2);j++)
  {
  x=256+(float)j/(float)tock*512;
  y=216-(graf[j]-miny)/(maxy-miny)*216;
  lineto(x,y);
  }
}
//-----------------------------------------------------------------

void risigraf2(float graf[], long tock, float miny, float maxy)
{
int x,y;
long j;
//line(0,0,getmaxx(),getmaxy());


setviewport(50,240,562,456,1);
clearviewport();

setcolor(GREEN);
moveto(0,0);lineto(512,0);lineto(512,216);lineto(0,216);lineto(0,0);
for (y=36;y<216;y=y+36) line(0,y,512,y);
for (x=52;x<500;x=x+51) line(x,0,x,216);
setcolor(WHITE);

x=0;y=216-(graf[(int)(tock/2)]-miny)/(maxy-miny)*216;
moveto(x,y);
for (j=1;j<(int)(tock/2);j++)
  {
  x=(float)j/(float)tock*512;
  y=216-(graf[j+(int)(tock/2)]-miny)/(maxy-miny)*216;
  lineto(x,y);
  }

x=256;y=216-(graf[0]-miny)/(maxy-miny)*216;
moveto(x,y);
for (j=1;j<(int)(tock/2);j++)
  {
  x=256+(float)j/(float)tock*512;
  y=216-(graf[j]-miny)/(maxy-miny)*216;
  lineto(x,y);
  }

}

//***********************************************************************

int main()
{
int gdriver = DETECT, gmode, errorcode;
int i,n,size;
char irqmask;
char st[30];
float srate;

size=16384;
n=size;
printf("\n\n--------------------SIDI1FFT---------------------\n");
srate=checkspeed();
fftmojprip(n,wr,wi,bi);
//printf("\n W=");pisikompl(5,wr,wi);
//for (i=0;i<n;i++) printf(" %d",bi[i]);

hitrovz(size);
preflintaj1(size,1);

//siggen(n,3,xr,xi);

//printf("\n x=");pisikompl(8,xr,xi);
fft(n,xr,xi,wr,wi,bi);
//printf("\n X=");pisikompl(8,xr,xi);
printf("\n\nPress ENTER to continue");
getch();

irqmask=inportb(0x21);
outportb(0x21,irqmask+1);      //stop timer interrupt


initgraph(&gdriver, &gmode, "c:\\bc\\bgi");

do
  {
  hitrovz(size);
  preflintaj1(size,1);
//  siggen(n,20,xr,xi);
  fft(n,xr,xi,wr,wi,bi);
  for (i=0;i<size;i++)
	{
	//graf[i]=sqrt((float)xr[i]*xr[i]+(float)xi[i]*xi[i]);
	graf[i]=20*log10(sqrt((float)xr[i]*xr[i]+(float)xi[i]*xi[i])+0.000000000001)-24.5;
	}
  risigraf1(graf,size,0,60);

  hitrovz(size);
  preflintaj1(size,2);
//  siggen(n,20,xr,xi);
  fft(n,xr,xi,wr,wi,bi);
  for (i=0;i<size;i++)
	{
	//graf[i]=sqrt((float)xr[i]*xr[i]+(float)xi[i]*xi[i]);
	graf[i]=20*log10(sqrt((float)xr[i]*xr[i]+(float)xi[i]*xi[i])+0.000000000001)-24.5;
	}
  risigraf2(graf,size,0,60);

  setviewport(0,0,639,479,1);

  outtextxy(275,1,"SIDI1FFT");

  outtextxy(10,20,"CH1");
  outtextxy(10,250,"CH2");

  sprintf(st,"-%2.3f MHz",srate/2000000);
  outtextxy(18,470,st);
  sprintf(st,"+%2.3f MHz",srate/2000000);
  outtextxy(532,470,st);

  }
while(kbhit()==0);

closegraph();

outportb(0x21,irqmask);

return 0;
}