C++ cexpf函数代码示例

/* Performs 1st order linear interpolation with out-of-bound interpolation */
void srslte_interp_linear_offset_cabs(cf_t *input, cf_t *output, 
                               uint32_t M, uint32_t len, 
                               uint32_t off_st, uint32_t off_end) 
{
  uint32_t i, j;
  float mag0=0, mag1=0, arg0=0, arg1=0, mag=0, arg=0;

  for (i=0;i<len-1;i++) {
    mag0 = cabsf(input[i]);
    mag1 = cabsf(input[i+1]);
    arg0 = cargf(input[i]);
    arg1 = cargf(input[i+1]);
    if (i==0) {
      for (j=0;j<off_st;j++) {
        mag = mag0 - (j+1)*(mag1-mag0)/M;
        arg = arg0 - (j+1)*(arg1-arg0)/M;
        output[j] = mag * cexpf(I * arg);
      }
    }
    for (j=0;j<M;j++) {
      mag = mag0 + j*(mag1-mag0)/M;
      arg = arg0 + j*(arg1-arg0)/M;
      output[i*M+j+off_st] = mag * cexpf(I * arg);
    }
  }
  if (len > 1) {
    for (j=0;j<off_end;j++) {
      mag = mag1 + j*(mag1-mag0)/M;
      arg = arg1 + j*(arg1-arg0)/M;
      output[i*M+j+off_st] = mag * cexpf(I * arg);
    }
  }
}

void autotest_nco_crcf_mix_block_up()
{
    // options
    unsigned int buf_len = 4096;
    float        phase   = 0.7123f;
    float        freq    = 0; //0.1324f;
    float        tol     = 1e-2f;

    // create object
    nco_crcf nco = nco_crcf_create(LIQUID_NCO);
    nco_crcf_set_phase    (nco, phase);
    nco_crcf_set_frequency(nco, freq);

    // generate signal
    float complex buf_0[buf_len];
    float complex buf_1[buf_len];
    unsigned int i;
    for (i=0; i<buf_len; i++)
        buf_0[i] = cexpf(_Complex_I*2*M_PI*randf());

    // mix signal
    nco_crcf_mix_block_up(nco, buf_0, buf_1, buf_len);

    // compare result to expected
    float phi = phase;
    for (i=0; i<buf_len; i++) {
        float complex v = buf_0[i] * cexpf(_Complex_I*phi);
        CONTEND_DELTA( crealf(buf_1[i]), crealf(v), tol);
        CONTEND_DELTA( cimagf(buf_1[i]), cimagf(v), tol);
        phi += freq;
    }

    // destroy object
    nco_crcf_destroy(nco);
}

/* Performs 1st order linear interpolation with out-of-bound interpolation */
void interp_linear_offset(cf_t *input, cf_t *output, int M, int len, int off_st, int off_end) {
	int i, j;
	float mag0, mag1, arg0, arg1, mag, arg;

	for (i=0;i<len-1;i++) {
		mag0 = cabsf(input[i]);
		mag1 = cabsf(input[i+1]);
		arg0 = cargf(input[i]);
		arg1 = cargf(input[i+1]);
		if (i==0) {
			for (j=0;j<off_st;j++) {
				mag = mag0 - (j+1)*(mag1-mag0)/M;
				arg = arg0 - (j+1)*(arg1-arg0)/M;
				output[j] = mag * cexpf(I * arg);
			}
		}
		for (j=0;j<M;j++) {
			mag = mag0 + j*(mag1-mag0)/M;
			arg = arg0 + j*(arg1-arg0)/M;
			output[i*M+j+off_st] = mag * cexpf(I * arg);
		}
	}
	if (len > 1) {
		for (j=0;j<off_end;j++) {
			mag = mag1 + j*(mag1-mag0)/M;
			arg = arg1 + j*(arg1-arg0)/M;
			output[i*M+j+off_st] = mag * cexpf(I * arg);
		}
	}
}

void myradon2_lop(bool adj, bool add, int nm, int nd, sf_complex *mm, sf_complex *dd)
/*< radon linear operator >*/
{
	int ix, ip;
	sf_complex sumc;

    	if (nm != np || nd != nx) sf_error("%s: mismatched data sizes",__FILE__);
	
	sf_cadjnull(adj, add, nm, nd, mm, dd);

	if(adj){// mm(p,w)=sum_{ix=0}^{nx} dd(xx[ix],w)*exp(i*w*p*xx[ix])
		for(ip=0; ip<np; ip++) // loop over slopes
		{
			sumc=sf_cmplx(0,0);
			for(ix=0; ix<nx; ix++) 
				sumc+=cexpf(I*w*p[ip]*xx[ix])*dd[ix];
			mm[ip]=sumc;
		}
	}else{// dd(xx,w)=sum_{ip=0}^{np} mm(p[ip],w)*exp(-i*w*p[ip]*xx)
		for(ix=0; ix<nx; ix++) 
		{
			sumc=sf_cmplx(0,0);
			for(ip=0; ip<np; ip++)
				sumc+=cexpf(-I*w*p[ip]*xx[ix])*mm[ip];
			dd[ix]=sumc;
		}
	}
}

/*  _m          :   number of correlators                   */
PRESYNC() PRESYNC(_create)(TC *         _v,
                           unsigned int _n,
                           float        _dphi_max,
                           unsigned int _m)
{
    // validate input
    if (_n < 1) {
        fprintf(stderr, "error: bpresync_%s_create(), invalid input length\n", EXTENSION_FULL);
        exit(1);
    } else if (_m == 0) {
        fprintf(stderr, "error: bpresync_%s_create(), number of correlators must be at least 1\n", EXTENSION_FULL);
        exit(1);
    }

    // allocate main object memory and initialize
    PRESYNC() _q = (PRESYNC()) malloc(sizeof(struct PRESYNC(_s)));
    _q->n = _n;
    _q->m = _m;

    _q->n_inv = 1.0f / (float)(_q->n);

    unsigned int i;

    // create internal receive buffers
    _q->rx_i = WINDOW(_create)(_q->n);
    _q->rx_q = WINDOW(_create)(_q->n);

    // create internal array of frequency offsets
    _q->dphi = (float*) malloc( _q->m*sizeof(float) );

    // create internal synchronizers
    _q->sync_i = (DOTPROD()*) malloc( _q->m*sizeof(DOTPROD()) );
    _q->sync_q = (DOTPROD()*) malloc( _q->m*sizeof(DOTPROD()) );

    // buffer
    T vi_prime[_n];
    T vq_prime[_n];
    for (i=0; i<_q->m; i++) {

        // generate signal with frequency offset
        _q->dphi[i] = (float)i / (float)(_q->m-1)*_dphi_max;
        unsigned int k;
        for (k=0; k<_q->n; k++) {
            vi_prime[k] = REAL( _v[k] * cexpf(-_Complex_I*k*_q->dphi[i]) );
            vq_prime[k] = IMAG( _v[k] * cexpf(-_Complex_I*k*_q->dphi[i]) );
        }

        _q->sync_i[i] = DOTPROD(_create)(vi_prime, _q->n);
        _q->sync_q[i] = DOTPROD(_create)(vq_prime, _q->n);
    }

    // allocate memory for cross-correlation
    _q->rxy = (float*) malloc( _q->m*sizeof(float) );

    // reset object
    PRESYNC(_reset)(_q);

    return _q;
}

static void init_wave(int init, 
		      int nx, float dx,
		      int nz, float dz,
		      sf_complex *pp /* [nx] */,
		      float wov, int nw, int iw)
{
    int ix;
    float x,x0,z0,phase,amp;
    
    x0 = nx*dx/3;
    z0 = nz*dz/3;

    switch(init) {
	case 1: /*  planar wave @ 15deg */
	    for (ix=0; ix < nx; ix++) {
		x = (ix+1)*dx - x0;
		phase  = wov*x*sinf(15*SF_PI/180.);
		pp[ix] = cexpf(sf_cmplx(0.,phase));
	    }
	    break;
	case 2: /* expanding spherical wave */
	    for (ix=0; ix < nx; ix++) {
		x = (ix+1)*dx - x0;
		phase  = wov*hypotf(z0,x);
		pp[ix] = cexpf(sf_cmplx(0.,phase));
	    }
	    break;
	case 3: /* point source */
	    for (ix=0; ix < nx; ix++) {
		pp[ix]=sf_cmplx(0.,0.);
	    }
	    pp[nx/3-1] = sf_cmplx(1.,0.);
	    break;
	case 4: /* collapsing spherical wave */
	    for (ix=0; ix < nx; ix++) {
		x = (ix+1)*dx - x0;
		phase  = -wov*hypotf(z0,x);
		pp[ix] = cexpf(sf_cmplx(0.,phase));
	    }
	    break;
	default:
	    sf_error("Unknown init=%d",init);
    }

    amp = (nw-iw+1.0)/nw;
    amp = cosf((1-amp)*(0.5*SF_PI));
    amp *= amp;

    for (ix=0; ix < nx; ix++) {
#ifdef SF_HAS_COMPLEX_H
	pp[ix] *= amp;
#else
	pp[ix] = sf_crmul(pp[ix],amp);
#endif
    }
}

// [AC-Adaptive]
static void radial_self_delays(unsigned int N, float shifts[N], const float phi[N], const long dims[DIMS], const complex float* in)
{
	unsigned int d = 2;
	unsigned int flags = (1 << d);

	assert(N == dims[d]);

	long dims1[DIMS];
	md_select_dims(DIMS, ~flags, dims1, dims);

	complex float* tmp1 = md_alloc(DIMS, dims1, CFL_SIZE);
	complex float* tmp2 = md_alloc(DIMS, dims1, CFL_SIZE);

	long pos[DIMS] = { 0 };

	for (unsigned int i = 0; i < dims[d]; i++) {

		pos[d] = i;
		md_copy_block(DIMS, pos, dims1, tmp1, dims, in, CFL_SIZE);

		// find opposing spoke

		float mdelta = 0.;
		int mindex = 0;

		for (unsigned int j = 0; j < dims[d]; j++) {

			float delta = cabsf(cexpf(1.i * phi[j]) - cexpf(1.i * phi[i]));

			if (mdelta <= delta) {

				mdelta = delta;
				mindex = j;
			}
		}

		pos[d] = mindex;
		md_copy_block(DIMS, pos, dims1, tmp2, dims, in, CFL_SIZE);


		unsigned int d2 = 1;
		float rshifts[DIMS];
		md_flip(DIMS, dims1, MD_BIT(d2), tmp2, tmp2, CFL_SIZE); // could be done by iFFT in est_subpixel_shift
		est_subpixel_shift(DIMS, rshifts, dims1, MD_BIT(d2), tmp2, tmp1);

		float mshift = rshifts[d2] / 2.; // mdelta

		shifts[i] = mshift;
	}

	md_free(tmp1);
	md_free(tmp2);
}

void xkolmog(sf_complex *trace1, sf_complex *trace2)
/*< convert Fourier-domain cross-correlation to minimum-phase >*/ 
{
    int i1;
    const double eps=1.e-32;

    for (i1=0; i1 < nk; i1++) {
#ifdef SF_HAS_COMPLEX_H
	fft1[i1] = clogf(trace1[i1]+eps)/nk;
#else
	fft1[i1] = sf_crmul(clogf(sf_cadd(trace1[i1],sf_cmplx(eps,0.))),
			    1.0/nk);
#endif
    }

    /* Inverse transform */
    kiss_fft(invs,(const kiss_fft_cpx *) fft1, (kiss_fft_cpx *) trace1);

#ifdef SF_HAS_COMPLEX_H
    trace1[0]    *= 0.5; trace2[0]    = trace1[0];
    trace1[nk/2] *= 0.5; trace2[nk/2] = trace1[nk/2];
#else
    trace1[0]    = sf_crmul(trace1[0],   0.5); trace2[0]    = trace1[0];
    trace1[nk/2] = sf_crmul(trace1[nk/2],0.5); trace2[nk/2] = trace1[nk/2];
#endif
    for (i1=1+nk/2; i1 < nk; i1++) {
	trace2[nk-i1] = trace1[i1];
	trace1[i1] = sf_cmplx(0.,0.);
	trace2[i1] = sf_cmplx(0.,0.);
    }

    /* Fourier transform */
    kiss_fft(forw,(const kiss_fft_cpx *) trace1, (kiss_fft_cpx *) fft1);
    kiss_fft(forw,(const kiss_fft_cpx *) trace2, (kiss_fft_cpx *) fft2);
    
    for (i1=0; i1 < nk; i1++) {
#ifdef SF_HAS_COMPLEX_H
	fft1[i1] = cexpf(fft1[i1])/nk;
	fft2[i1] = cexpf(fft2[i1])/nk;
#else
	fft1[i1] = sf_crmul(cexpf(fft1[i1]),1./nk);
	fft2[i1] = sf_crmul(cexpf(fft2[i1]),1./nk);
#endif
    }

    /* Inverse transform */
    kiss_fft(invs,(const kiss_fft_cpx *) fft1, (kiss_fft_cpx *) trace1);
    kiss_fft(invs,(const kiss_fft_cpx *) fft2, (kiss_fft_cpx *) trace2);

    for (i1=0; i1 < nk; i1++) {
	trace2[i1] = conjf(trace2[i1]);
    }
}

sf_complex Green(float r1,float r2,float r3,float s1,float s2,float s3,float omega) 
/*< Green's function >*/
{
    double tt,amp;
    sf_complex val;

    GreenTtAmp(r1,r2,r3,s1,s2,s3,&tt,&amp);
#ifdef SF_HAS_COMPLEX_H
    val=amp*cexpf(sf_cmplx(0.,omega*tt)); 
#else
    val=sf_crmul(cexpf(sf_cmplx(0.,omega*tt)),amp); 
#endif
    return (val);
}

void BornScatteredField(float xx,float xy,float xz,float *output)
/*< Born scattering field >*/
{
    int iw;
    double omega,scale;

    sf_complex  U,dU;
    sf_complex val,fkern;

    ZeroArray(u,  nt); 
    ZeroArray(du, nt); 
    *output=0.;

    for (iw=0; iw<nw; iw++) {
	omega=ow+dw*iw;

	scale =cos((SF_PI/2)*(((double) iw+1)/((double) nw+1)));
	scale*=scale*omega;

#ifdef SF_HAS_COMPLEX_H
	/* background field */
	U=Green(xx,xy,xz,sx,sy,sz,omega)*scale;
	    
	/* scattered field */
	val=Green(px,py,pz,sx,sy,sz,omega)*scale;
	val *= Green(xx,xy,xz,px,py,pz,omega);
	dU = val*(-omega*omega*dv);

	U += dU;
	fkern=cexpf(sf_cmplx(0.,-omega*0.6));
	val=fkern*U;
#else
	/* background field */
	U=sf_crmul(Green(xx,xy,xz,sx,sy,sz,omega),scale);
	    
	/* scattered field */
	val=sf_crmul(Green(px,py,pz,sx,sy,sz,omega),scale);
	val=sf_cmul(val,Green(xx,xy,xz,px,py,pz,omega));
	dU=sf_crmul(val,-omega*omega*dv);

	U=sf_cadd(U,dU);
	fkern=cexpf(sf_cmplx(0.,-omega*0.6));
	val=sf_cmul(fkern,U);
#endif

	*output+=crealf(val);
    }
    return;
}    

// 
// Test AC gain control
//
void autotest_agc_crcf_ac_gain_control()
{
    // set paramaters
    float gamma = 0.1f;             // nominal signal level
    float bt    = 0.1f;             // bandwidth-time product
    float tol   = 0.001f;           // error tolerance
    float dphi  = 0.1f;             // NCO frequency

    // create AGC object and initialize
    agc_crcf q = agc_crcf_create();
    agc_crcf_set_bandwidth(q, bt);

    unsigned int i;
    float complex x;
    float complex y;
    for (i=0; i<256; i++) {
        x = gamma * cexpf(_Complex_I*i*dphi);
        agc_crcf_execute(q, x, &y);
    }

    if (liquid_autotest_verbose)
        printf("gamma : %12.8f, rssi : %12.8f\n", gamma, agc_crcf_get_signal_level(q));

    // Check results
    CONTEND_DELTA( agc_crcf_get_gain(q), 1.0f/gamma, tol);

    // destroy AGC object
    agc_crcf_destroy(q);
}

void fbdip_init(float rad, int m1, int m2, 
	int *rect, int niter, float dip0, bool vb)
/*< initialize >*/
{
	int n, nn[2];
	n1 = m1;
	n2 = m2;
	verb = vb;

	u1 = sf_floatalloc2(n1, n2);
	u2 = sf_floatalloc2(n1, n2);
	u3 = sf_floatalloc2(n1, n2);
	u4 = sf_floatalloc2(n1, n2);
	u5 = sf_floatalloc2(n1, n2);
	p = sf_complexalloc2(n1, n2);

	r=rad;
	p0 = rad*cexpf(sf_cmplx(0, dip0));
	
	if(rect[0]>0 && rect[1]>0)
	{
		n = n1*n2;
		nn[0] = n1;
		nn[1] = n2;
		sf_divn_init (2, n, nn, rect, niter, false);
		use_divn=true;
	}else 	use_divn=false;
}

// create freqmod object
//  _kf     :   modulation factor
FREQMOD() FREQMOD(_create)(float _kf)
{
    // validate input
    if (_kf <= 0.0f || _kf > 1.0) {
        fprintf(stderr,"error: freqmod_create(), modulation factor %12.4e out of range [0,1]\n", _kf);
        exit(1);
    }

    // create main object memory
    FREQMOD() q = (freqmod) malloc(sizeof(struct FREQMOD(_s)));

    // set modulation factor
    q->kf  = _kf;
    q->ref = q->kf * (1<<16);

    // initialize look-up table
    q->sincos_table_len = 1024;
    q->sincos_table     = (TC*) malloc( q->sincos_table_len*sizeof(TC) );
    unsigned int i;
    for (i=0; i<q->sincos_table_len; i++) {
        q->sincos_table[i] = cexpf(_Complex_I*2*M_PI*(float)i / (float)(q->sincos_table_len) );
    }

    // reset modem object
    FREQMOD(_reset)(q);

    // return object
    return q;
}

void benchmark_gmskmodem_demodulate(struct rusage *_start,
                                    struct rusage *_finish,
                                    unsigned long int *_num_iterations)
{
    // options
    unsigned int k=2;   // filter samples/symbol
    unsigned int m=3;   // filter delay (symbols)
    float BT=0.3f;      // bandwidth-time product

    // create modem object
    gmskdem demod = gmskdem_create(k, m, BT);

    float complex x[k];
    unsigned int symbol_out = 0;
    
    unsigned long int i;
    for (i=0; i<k; i++)
        x[i] = randnf()*cexpf(_Complex_I*2*M_PI*randf());

    // start trials
    getrusage(RUSAGE_SELF, _start);
    for (i=0; i<(*_num_iterations); i++) {
        gmskdem_demodulate(demod, x, &symbol_out);
        gmskdem_demodulate(demod, x, &symbol_out);
        gmskdem_demodulate(demod, x, &symbol_out);
        gmskdem_demodulate(demod, x, &symbol_out);
    }
    getrusage(RUSAGE_SELF, _finish);
    *_num_iterations *= 4;

    // destroy modem objects
    gmskdem_destroy(demod);
}

/*------------------------------------------------------------*/
void cam3_ssf(
    sf_complex       w /* frequency */,
    sf_complex   ***wx /* wavefield */,
    cub3d          cub,
    cam3d          cam,
    tap3d          tap,
    slo3d          slo,
    int            imz,
    int         ompith
    )
/*< Wavefield extrapolation by SSF >*/
{
    sf_complex co,cs,cr,w2,khy,kss,krr;
    float s, kmy, d,dsc,drc;
    int imy,imx,ihx,jmy, js,jr;

    co = sf_cmplx(0,0);
#ifdef SF_HAS_COMPLEX_H
    w2 = w*w;
#else
    w2 = sf_cmul(w,w);
#endif
    
#ifdef SF_HAS_COMPLEX_H
    LOOP( s = slo->so[ompith][ cam->jy[imy] ][ cam->is[ihx][imx] ] 
	  +   slo->so[ompith][ cam->jy[imy] ][ cam->ir[ihx][imx] ];
	  wx[ihx][imy][imx] *= cexpf(-w*s* cub->amz.d/2); );