demodulatormex.c

本程序利用改进的CHOW算法

#include "mex.h"
#include "mat.h"
#include <math.h>
#define  bits_out plhs[0]
#define eie0_out plhs[1]
#define eie1_out plhs[2]
#define	bits_in prhs[0]
#define baud_in prhs[1]
#define number_of_bits prhs[2]
#define	tag_in prhs[3]

//function bits_out=demodulatormex(samples,numberofbits,tag);

//Pass in a modulated waveform to be digitally demodulated.
//Also, specify the baudrate. Finally, specify the type of modulation to use.
//tag = 1 (FSK), 2 (BPSK), 3 (QPSK), 4 (16PSK), 5 (16QAM) - unimplemented, 6 (DPSK) unimplemented
//returns message (output_bits) if FSK,BPSK {0,1} if QPSK ......

// This function was written by me, Joe Williams, a Stanford EE Grad Student. It is used for my 
// EE359 project (Wireless Comm). You may use this code in whatever manner you want as long as
// you give me credit.

void mexFunction(int nlhs, 
                     mxArray *plhs[],
                     int nrhs, 
                     const mxArray *prhs[])
    {
	
	  double *output_bits,*eie0,*eie1;
	  double *input_bits; 
      MATFile *ph;
	  double baudrate,fa,fb,Fs,num_bits,totaldatasamples,v0a,v0b,v1a,v1b;
      int buflen,tag,samplesperbit,i,j,inputlength,bit;
	  double f0,f1,wave0[100000],wave1[100000],wave0c[100000],wave1c[100000],wave2[100000],wave3[100000],wave[16][100000],val,tempval,roundval;
      double *var;
    
	  baudrate=mxGetScalar(baud_in);
      fa=900000000;
      fb=899900000;
	  tag=mxGetScalar(tag_in);
	  num_bits=mxGetScalar(number_of_bits);
	  Fs=9000000*20; // * 100
      input_bits = mxGetPr(bits_in);
	  inputlength = mxGetN(bits_in);
	  samplesperbit=Fs/(baudrate/100);
	  //mexPrintf("%d",samplesperbit);
      //samplesperbit=9000; 
	  totaldatasamples = samplesperbit*num_bits;
      bits_out = mxCreateDoubleMatrix(1, num_bits, mxREAL);
      output_bits = mxGetPr(bits_out);
      eie0_out = mxCreateDoubleMatrix(1, num_bits, mxREAL);
      eie0 = mxGetPr(eie0_out);
	  eie1_out = mxCreateDoubleMatrix(1, num_bits, mxREAL);
      eie1 = mxGetPr(eie1_out);

	  f0=fa/Fs*.02*3.1415926539;
	  f1=fb/Fs*.02*3.1415926539;
	 
	  switch(tag){//construct 1 period of signal
	  case 0://FSK noncoh
		for(i=0;i<samplesperbit;i++){
	      wave0[i]=sin(f0*(float)i);
		  wave0c[i]=cos(f0*(float)i);
		}
		for(i=0;i<samplesperbit;i++){
	      wave1[i]=sin(f1*(float)i);
		  wave1c[i]=sin(f1*(float)i);
		}
	    break;
	  case 1://FSK coh
	    for(i=0;i<samplesperbit;i++)
	      wave0[i]=sin(f0*(float)i);
	    for(i=0;i<samplesperbit;i++)
	      wave1[i]=sin(f1*(float)i);
	    break;
	  case 2://BPSK
		for(i=0;i<samplesperbit;i++)
	      wave0[i]=sin(f0*(float)i);
	    for(i=0;i<samplesperbit;i++)
	      wave1[i]=sin(f0*(float)i+3.1415926539);
	    break;
	  case 3://QPSK
		for(i=0;i<samplesperbit;i++)
		  wave0[i]=sin(f0*(float)i);
		for(i=0;i<samplesperbit;i++)
		  wave1[i]=cos(f0*(float)i);
        break;
	  case 4://16PSK
		for(i=0;i<samplesperbit;i++)
		  wave0[i]=sin(f0*(float)i);
		for(i=0;i<samplesperbit;i++)
		  wave1[i]=cos(f0*(float)i);
        break;
      }
		

	  //deconstruct signal
	  if(tag==0){//FSK non-coherent
       /* for(i=0;i<num_bits;i++){
		  v0a=0.0;v1a=0.0;v0b=0.0;v1b=0.0;
	      for(j=0;j<samplesperbit/100;j++){
		    v0a=v0a+(input_bits[i*samplesperbit+j])*(wave0[j]);
		    v1a=v1a+(input_bits[i*samplesperbit+j])*(wave1[j]);
		    v0b=v0b+(input_bits[i*samplesperbit+j])*(wave0c[j]);
		    v1b=v1b+(input_bits[i*samplesperbit+j])*(wave1c[j]);
		  }
		  output_bits[i]=0;
		  if((pow(v0a,2)+pow(v0b,2))<(pow(v1a,2)+pow(v1b,2))){
		    output_bits[i]=1;
			mexPrintf("%f vs %f\n",pow(v1a,2)+pow(v1b,2),pow(v0a,2)+pow(v0b,2));
		  }
		}*/
        mexPrintf("Not working for FA and FB so close!");
	  }
      if(tag==1){//FSK coherent
        for(i=0;i<num_bits;i++){
		  v0a=0.0;v1a=0.0;
		  for(j=0;j<samplesperbit;j++){
			v0a=v0a+(input_bits[i*samplesperbit+j])*(wave0[j]);
		    v1a=v1a+(input_bits[i*samplesperbit+j])*(wave1[j]);
          }
		  output_bits[i]=0;eie0[i]=v0a;eie1[i]=v1a;
		  if(v1a>v0a) output_bits[i]=1;
		}
	  }
	  else if(tag==2){//BPSK
        for(i=0;i<num_bits;i++){
		  v0a=0.0;v1a=0.0;
		  for(j=0;j<samplesperbit;j++){
			v0a=v0a+(input_bits[i*samplesperbit+j])*(wave0[j]);
		    v1a=v1a+(input_bits[i*samplesperbit+j])*(wave1[j]);
          }
		  output_bits[i]=0;
		  if(v1a>v0a) output_bits[i]=1;
		  eie0[i]=v0a;eie1[i]=v1a;
		}
	  }
	  else if(tag==3){//QPSK
        for(i=0;i<num_bits;i++){
		  v0a=0.0;v1a=0.0;
		  for(j=0;j<samplesperbit;j++){
			v0a=v0a+(input_bits[i*samplesperbit+j])*(wave0[j]);
		    v1a=v1a+(input_bits[i*samplesperbit+j])*(wave1[j]);
          }
		  output_bits[i]=0;
          if((v0a<0.0)&&(v1a<0.0)) output_bits[i]=2;
          else if((v0a<0)&&(v1a>0)) output_bits[i]=1;
          else if((v0a>0)&&(v1a<0)) output_bits[i]=3;
		  eie0[i]=v0a;eie1[i]=v1a;
		}
	  }
	  else if(tag==4){//16PSK
        for(i=0;i<num_bits;i++){
		  v0a=0.0;v1a=0.0;
		  for(j=0;j<samplesperbit;j++){
			v0a=v0a+(input_bits[i*samplesperbit+j])*(wave0[j]);
		    v1a=v1a+(input_bits[i*samplesperbit+j])*(wave1[j]);
          }
		  eie0[i]=v0a;eie1[i]=v1a;
		  val=atan(v1a/v0a);
		  tempval=8.0*val/3.1415926539-(int)(8.0*val/3.1415926539);
		  if(tempval>0.0){
			if(tempval<0.5)
			  roundval=(int)(8.0*val/3.1415926539);
			else
			  roundval=(int)(8.0*val/3.1415926539)+1;
		  }
		  else{
		    if(tempval>-0.5)
			  roundval=(int)(8.0*val/3.1415926539);
			else
			  roundval=(int)(8.0*val/3.1415926539)-1;
		  }
          output_bits[i]=roundval;
          if((v0a<0)&&(v1a<0)) output_bits[i]+=8;
          else if((v0a>0)&&(v1a<0)) output_bits[i]+=16;
          else if((v0a<0)&&(v1a>0)) output_bits[i]+=8;
          if(output_bits[i]==16) output_bits[i]=0;
		}
	  }
	  //mexPrintf("Done demodulating signal.\n");
}