C++ cabs函数代码示例

double ov_matrix12x12_rowsumnorm(matrix12x12 A) {
	
  double norm, nrm;
  int i, j;

  norm = 0.0;
  for (i=0; i<12; i++) {

    nrm = 0.0;
    for (j=0; j<12; j++)
      nrm += cabs(A[i][j]);

    if (nrm > norm)
      norm = nrm;
  }

  return norm;
}

// Compute Real Cepstrum Of Signal
void RealCepstrum(int n, float *signal, float *realCepstrum)
{
  float *realTime, *imagTime, *realFreq, *imagFreq;
  int i;

  realTime = new float[n];
  imagTime = new float[n];
  realFreq = new float[n];
  imagFreq = new float[n];

  // Compose Complex FFT Input

  for(i = 0; i < n; i++)
  {
    realTime[i] = signal[i];
    imagTime[i] = 0.0f;
  }

  // Perform DFT

  DFT(n, realTime, imagTime, realFreq, imagFreq);

  // Calculate Log Of Absolute Value

  for(i = 0; i < n; i++)
  {
    realFreq[i] = logf(cabs(realFreq[i], imagFreq[i]));
    imagFreq[i] = 0.0f;
  }

  // Perform Inverse FFT

  InverseDFT(n, realTime, imagTime, realFreq, imagFreq);

  // Output Real Part Of FFT
  for(i = 0; i < n; i++)
    realCepstrum[i] = realTime[i];

  delete realTime;
  delete imagTime;
  delete realFreq;
  delete imagFreq;
}

PetscErrorCode DiffIE(Vec originalVec, Vec reducedVec, Vec Distance, Vec nParticlePerCube)
{
    // Find the distance between 2 solutions originalVec and reducedVec
    PetscInt    s,t;
    PetscInt    n,xstart,xend;
    PetscScalar tmp,v,diff;
    Vec         vec;
    VecScatter  scat;
    IS          is;

    PetscFunctionBegin;
    VecGetOwnershipRange(originalVec,&xstart,&xend);
    VecGetSize(reducedVec,&n);

    VecCreate(PETSC_COMM_SELF,&vec);
    VecSetSizes(vec,PETSC_DECIDE,n);
    VecSetFromOptions(vec);

    ISCreateStride(PETSC_COMM_WORLD,n,0,1,&is);
    VecScatterCreate(reducedVec,PETSC_NULL,vec,is,&scat);
    //VecScatterCreateToAll(reducedVec,&scat,&vec);
    VecScatterBegin(scat,reducedVec,vec,INSERT_VALUES,SCATTER_FORWARD);
    VecScatterEnd(scat,reducedVec,vec,INSERT_VALUES,SCATTER_FORWARD);

    for(s=xstart;s<xend;s++)
    {
        t = ScatIE.FindCube(s);  
        VecGetValues(vec,1,&t,&tmp);
        VecGetValues(originalVec,1,&s,&v);

        diff = cabs(v - tmp);
        VecSetValues(Distance,1,&t,&diff,ADD_VALUES);
        tmp = 1;
        VecSetValues(nParticlePerCube,1,&t,&tmp,ADD_VALUES);        
    }   
  
    VecDestroy(&vec);
    VecScatterDestroy(&scat);
    ISDestroy(&is);

    PetscFunctionReturn(0);
}

/* writes a time-domain waveform to stdout as tab-separated values */
int output_td_waveform(REAL8TimeSeries * h_plus, REAL8TimeSeries * h_cross, struct params p)
{
    double t0;
    size_t j;
    t0 = XLALGPSGetREAL8(&h_plus->epoch);
    if (p.amp_phase) {
        REAL8Sequence *amp = XLALCreateREAL8Sequence(h_plus->data->length);
        REAL8Sequence *phi = XLALCreateREAL8Sequence(h_plus->data->length);
        double phi0;

        /* compute the amplitude and phase of h+ - i hx */
        for (j = 0; j < h_plus->data->length; ++j) {
            double complex z = h_plus->data->data[j] - I * h_cross->data->data[j];
            amp->data[j] = cabs(z);
            phi->data[j] = carg(z);
        }

        /* unwrap the phase */
        XLALREAL8VectorUnwrapAngle(phi, phi);

        /* make phase in range -pi to +pi at end of waveform */
        /* extrapolate the end of the waveform using last and second last points */
        phi0 = 2 * phi->data[phi->length - 1] - phi->data[phi->length - 2];
        // phi0 = phi->data[phi->length - 1];
        phi0 -= fmod(phi0 + copysign(LAL_PI, phi0), 2.0 * LAL_PI) - copysign(LAL_PI, phi0);
        for (j = 0; j < phi->length; ++j)
            phi->data[j] -= phi0;

        fprintf(stdout, "# time (s)\th_abs (strain)\t h_arg (rad)\n");
        for (j = 0; j < h_plus->data->length; ++j)
            fprintf(stdout, "%.9f\t%e\t%e\n", t0 + j * h_plus->deltaT, amp->data[j], phi->data[j]);

        XLALDestroyREAL8Sequence(phi);
        XLALDestroyREAL8Sequence(amp);
    } else {
        fprintf(stdout, "# time (s)\th_+ (strain)\th_x (strain)\n");
        for (j = 0; j < h_plus->data->length; ++j)
            fprintf(stdout, "%.9f\t%e\t%e\n", t0 + j * h_plus->deltaT, h_plus->data->data[j], h_cross->data->data[j]);
    }
    return 0;
}

// calculates the inverse of the matrix M outputs to M_1 //
int 
inverse( GLU_complex M_1[ NCNC ] , 
	 const GLU_complex M[ NCNC ] )
{
#if (defined HAVE_LAPACKE_H)
  const int n = NC , lda = NC ;
  int ipiv[ NC + 1 ] ;
  memcpy( M_1 , M , NCNC * sizeof( GLU_complex ) ) ;
  int info = LAPACKE_zgetrf( LAPACK_ROW_MAJOR , n , n , M_1 , lda , ipiv ) ;
  info = LAPACKE_zgetri( LAPACK_ROW_MAJOR , n , M_1 , lda, ipiv ) ;
  return info ;
#elif (defined CLASSICAL_ADJOINT_INV)
  // define the adjunct //
  GLU_complex adjunct[ NCNC ] GLUalign ; 
  register GLU_complex deter = cofactor_transpose( adjunct , M ) ;
  // here we worry about numerical stability //
  if( cabs( deter ) < NC * PREC_TOL ) {
    fprintf( stderr , "[INVERSE] Matrix is singular !!! "
	     "deter=%1.14e %1.14e \n" , creal( deter ) , cimag( deter ) ) ; 
    write_matrix( M ) ; 
    return GLU_FAILURE ; 
  }
  // obtain inverse of M from 1/( detM )*adj( M ) //
  size_t i ;
  deter = 1.0 / deter ;
  for( i = 0 ; i < NCNC ; i++ ) { M_1[i] = adjunct[i] * deter ; }
  return GLU_SUCCESS ;
#else
  #if NC == 2 
  // use the identity, should warn for singular matrices
  const double complex INV_detM = 1.0 / ( det( M ) ) ;
  M_1[ 0 ] =  M[ 3 ] * INV_detM ;
  M_1[ 1 ] = -M[ 1 ] * INV_detM ;
  M_1[ 2 ] = -M[ 2 ] * INV_detM ;
  M_1[ 3 ] =  M[ 0 ] * INV_detM ;
  return GLU_SUCCESS ;
  #else 
  return gauss_jordan( M_1 , M ) ;
  #endif
#endif 
}

static double_complex v_matrix_element(unsigned *total_num_eval,
                                       const struct integrand_params *params,
                                       double r_max, double abs_err,
                                       double rel_err, unsigned limit)
{
    double est_err;
    unsigned num_eval;
    const double_complex integral = gk_cquad(integrand, params, 0, r_max,
                                             abs_err, rel_err, limit,
                                             &est_err, &num_eval);
    const double_complex element = 2 / M_PI * integral;
    if (!(est_err < dmax(abs_err, rel_err * cabs(integral)))) {
        fprintf(stderr, "WARNING: integral did not converge: "
                "%i evals; error is %g, required %g\n",
                num_eval, est_err, abs_err);
        fflush(stderr);
    }
    if (total_num_eval)
        *total_num_eval += num_eval;
    return element;
}

double compute_mandelbrot(double _Complex c, int maxiter){

    double _Complex z = 0;
    int value = 0;
    int  diverge = 0;
    int n;
    for (n=0; n < maxiter; n++){
        z = cpow (z, 2) + c;
        if (cabs(z) > 2){
            value =  n;
            diverge = 1;
            break;
        }
    }
    if (diverge == 1){
        return value;
    }
    else{
        return 2;
    }
}

void CT_Freq_Res_mag(double *num,double *den,int N,int points,double *response,double Wmax)
{
	double complex response_num[points];
	double complex response_den[points];
	double w;
	int i,j;
	double complex S;
	for(w=0,j=0;w< Wmax; w= w+(Wmax/points)){
		S = 1;
		response_num[j] = 0;
		response_den[j] = 0;
		for(i=N-1;i>=0;i--){
			response_num[j] += (num[i]*S);
			response_den[j] += (den[i]*S);
			S *= I*w;
		}
		response[j] = cabs(response_num[j]/response_den[j]);
		j++;

	}
}

Signal::Signal() {
    value.resize(N);
    int n;
    for (n = 0; n < N; n++) {
        value[n] =  A0 * cos(phase(n)) + A0 * sin(phase(n)) * _Complex_I;
    }

    ofstream infile;
    infile.open("fq.dat");
    for (int n = 1; n < N; n++) {
        infile << (phase(n) - phase(n - 1)) / (Ts * 2 * M_PI) << " " << (n * Ts) << endl;
        frequency.push_back((phase(n) - phase(n - 1)) / (Ts * 2 * M_PI));
    }
    infile.close();

    infile.open("s.dat");
    for (int n = 0; n < N; n++) {
        infile << cabs(value[n]) * cos(carg(value[n])) << " " << (n * Ts) << endl;
    }
    infile.close();
}

int main()
{
  FILE *fp_r = fopen("Eth_r.txt", "r");
  FILE *fp_i = fopen("Eth_i.txt", "r");

  Complex datas[n];

  memset(datas, 0, sizeof(Complex)*n);
  for(int j=0; j<1500; j++)
  {
    double rm, im;
    fscanf(fp_r, "%lf", &rm);
    fscanf(fp_i, "%lf", &im);
    datas[j]   = rm + I*im;
  }
  
  printf("read file\n");

  FILE *outR = fopen("re.txt", "w");
  FILE *outI = fopen("im.txt", "w");
  for(int i=0; i<1500; i++){
    fprintf(outR, "%lf\n", creal(datas[i]));
    fprintf(outI, "%lf\n", cimag(datas[i]));
  }
  fclose(outR);
  fclose(outI);
  
  Complex *res = cdft(datas, n);
  
  FILE *fp = fopen("rest2.txt", "w");
  for(int i=0; i<1500; i++)
  {
    fprintf(fp, "%lf\n", cabs(res[i]) / sqrt(n));
  }

free(res);
  fclose(fp);  
  printf("finish\n");
  return 0;
}

int laguer(double complex a[], const int m, double complex *x, int *its, const int maxit) {
  int iter, i, j;
  double abx, abp, abm, err;
  double complex dx,x1,b,d,f,g,h,sq,gp,gm,g2;
  static double frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0};
  for (iter = 1; iter <= maxit; iter++) { 
    *its = iter;
    b = a[m];
    err = cabs(b);
    d = 0.;
    f = 0.;
    abx = cabs(*x);
    for (j = m-1; j >= 0; j--) {
      f = (*x) * f + d;
      d = (*x) * d + b;
      b = (*x) * b + a[j];
      err = cabs(b) + abx * err;
    }
    err *= epss;
    if (cabs(b) <= err) return(0);
    g = d / b;
    g2 = g * g;
    h = g2 - 2. * f / b;
    sq = csqrt((double)(m-1) * ((double)(m)*h - g2));
    gp = g + sq;
    gm = g - sq;
    abp = cabs(gp);
    abm = cabs(gm);
    if (abp < abm) gp = gm;
    dx=((dmax(abp,abm) > 0. ?
	 ((double complex)(m))/gp :
	 (1. + abx)*(cos((double)iter) + _Complex_I*sin((double)iter))));
    x1 = (*x) - dx;
    
    if (creal(*x) == creal(x1) && cimag(*x) == cimag(x1)) {
      return(0);
    }
    if (iter % MT) {
      *x=x1;
    }
    else {
      *x = (*x) - frac[iter/MT]*dx;
    }
  }
  fprintf(stderr, "Too many iterations in laguer\n");
  return(-1);
}

//Calculate noise buffer
void noise_estimate(){
	int i;
	float mag,s_factor;
	for (i=0;i<FFTLEN;i++)
	{    
		s_factor=1;      
		mag=cabs(buffer[i]);               
		N_i[i] = M1[i];
		if(N_i[i]>M2[i]){N_i[i]=M2[i];}
		if(N_i[i]>M3[i]){N_i[i]=M3[i];}
		if(N_i[i]>M4[i]){N_i[i]=M4[i];}
		SNR20=2-mag/N_i[i];//first order maclaurin
		SNR=SNR20*20;
		if(SNR>20){SNR=20;}
		if(SNR<0){SNR=0;}
		SNR_GAIN=-0.2*SNR+6;
		if(add1==1){
		if(i<=2||(i>=30&&i<=223)||i>=253){s_factor=50;}
		}
		N[i]=SNR_GAIN*s_factor*alpha*(1-lpf_weight_noise)*N_i[i]+lpf_weight_noise*N[i];
	}	
}

double complex
cpow(double complex a, double complex z)
{
	double complex w;
	double x, y, r, theta, absa, arga;

	x = creal(z);
	y = cimag(z);
	absa = cabs(a);
	if (absa == 0.0) {
		return (0.0 + 0.0 * I);
	}
	arga = carg(a);
	r = pow(absa, x);
	theta = x * arga;
	if (y != 0.0) {
		r = r * exp(-y * arga);
		theta = theta + y * log(absa);
	}
	w = r * cos(theta) + (r * sin(theta)) * I;
	return w;
}

//Calculate noise buffer
void noise_estimate(){
	int i;
	float mag,s_factor;
	for (i=0;i<FFTLEN;i++)
	{    
		s_factor=1;      
		mag=cabs(buffer[i]);               
		N_i[i] = M1[i];
		if(N_i[i]>M2[i]){N_i[i]=M2[i];}
		if(N_i[i]>M3[i]){N_i[i]=M3[i];}
		if(N_i[i]>M4[i]){N_i[i]=M4[i];}
		SNR20=2-mag/N_i[i];//first order maclaurin
		SNR=SNR20*20;
		if(SNR>10){SNR=10;}
		if(SNR<-5){SNR=-5;}
		alpha=-0.2*SNR+5;
		if(add1==1){
		if(i<=5||i>=80){s_factor=2;}
		}
		N[i]=s_factor*alpha*(1-lpf_weight_noise)*N_i[i]+lpf_weight_noise*N[i];
	}	
}

// NOTE: S is both an input AND output.
// Combines eigenspectra, fills S, but also computes the normalized absolute difference from what
// was previously in S.  This difference = || S_new - S_old ||/N is returned in diff.
void combine_eig_spec( const fftw_complex *Sk, uint_t K, uint_t N, const double *wk, uint_t nwk, fftw_complex *S, double &diff){
    
    uint_t iwt; //index of weight to use
    
    diff=0;
    double absdiff;
    fftw_complex Sl;

    for (uint_t ii=0; ii<N; ii++){
        Sl=S[ii];                               //stores previous value
        S[ii]=0;
        for ( uint_t jj=0; jj< K; jj++){
            iwt=jj*nwk+(ii%nwk);
            S[ii]=S[ii]+wk[iwt]*Sk[jj*N+ii];   //weighted average
        }
        absdiff=cabs(S[ii]-Sl);          //absolute difference
        diff=diff+absdiff;
    }
    
    diff=sqrt(diff)/N;
    
}