integralapi.cpp

本程序利用VC实现均匀量化和非均匀量化信噪比的仿真。并画出曲线

//integralapi.cpp

#include "stdafx.h"
#include "integralapi.h"
#include <io.h>
#include <errno.h>

#include <math.h>
#include <direct.h>

/////////////////////////////////////////////////////////////////
  double   WINAPI  f1(double x,double q,double a)
  {
    double f;
	f=(x-q)*(x-q)/a;   
	//确定要积分的函数，
	//即量化误差平方与密度函数的积


	return (f);
  }


/////////对f1函数进行在a1到a2的定积分/////////////////////////

double   WINAPI  integral1(double a1,double a2,double q,double a) 


  {
    
	
	double s,h,y;
	int n,i;
	s=(f1(a1,q,a)+f1(a2,q,a))/2.0;
	n=100;
	h=(a2-a1)/n;
	for(i=1;i<n;i++)
		s=s+f1((a1+i*h),q,a);
	y=s*h;
	return (y);
 
  }


////////////////////////////////////////////////////
///////////////////////////计算均匀量化的量化信噪比

double  WINAPI  snruniform(double a,int m)

{
    int d=4096/m;    //均匀量化间隔
	int a3;    //
	a3=int(a/d);   //确定a所在的位置哪个段
    int mi,mj,q,i;
    double nq,s,jf,snr1,snr;
	nq=0;
	mi=0;

   //实现g1（x）
	for(i=0;i<a3;i++)
	{
      mj=mi+d;
      q=(mi+mj)/2;
      jf=integral1(mi,mj,q,a) ;
	  nq=nq+jf;
	  mi=mi+d;

	}
    //实现g2（x）
      q=mi+d/2;
	  jf=integral1(mi,a,q,a) ;
	  nq=nq+jf;    //量化噪声功率


	  s=integral1(0,a,0,a) ;  //信号功率
	  
      
	  snr1=s/nq;    //信号量化噪声功率比
	  snr=10*log10(snr1);  //化为对数形式，log10函数是10为底的对数

    return(snr);

}


//////////////////////////////////////////////////////////
/////////对输入的信号进行13折线的变换，从x变化为y///////////
double  WINAPI  alaw(double x)
{
     double y;
     int a1=0;        //大段数  从0至7   共8级
	 int m;  // m大段起点  
	 
	 ///////计算a所在的大段的位置  
     m=16;
	 while(m<=x)
	 {
      m=2*m;
	  a1=a1+1;
	 }
    if(a1==0)
		y=16*x;
	else
	
        y=16*pow(2,(1-a1))*x+256*(a1-1);
	return(y);

}


//////////////////////////////////////////////////////////
////////////////////////计算A律十三折线量化的量化信噪比//////
///////////a为输入的最大幅度

   double  WINAPI  snralaw(double a)
   {
      
	 double jf,y;       //
	 int q,i;        //
     int mi=0;        //
	 int mj;          //
	 double nq;          //  噪声功率
     double s,snr1,snr;  //s信号功率，snr1 信号噪声功率比，snr以dB为单位
   
     y=alaw(a);   //对输入的信号进行13折线的变换，从a变化为y
     int d=4906/256;    //均匀量化间隔
	 int a3;    //
	 a3=int(y/d);   //确定a所在的位置哪个段   

	 nq=0;
   	mi=0;


	for(i=0;i<a3;i++)
	{
      mj=mi+d;
      q=(mi+mj)/2;
      jf=integral1(mi,mj,q,a) ;
	  nq=nq+jf;
	  mi=mi+d;

	}

      q=mi+d/2;
	  jf=integral1(mi,y,q,a) ;
	  nq=nq+jf;    //量化噪声功率


	  s=integral1(0,y,0,a) ;  //信号功率
	  
      
	  snr1=s/nq;    //信号量化噪声功率比
	  snr=10*log10(snr1);  //化为对数形式，log10函数是10为底的对数

    return(snr);

     
   }