C++ crealf函数代码示例

void wlanframesync_debug_print(wlanframesync _q,
                               const char * _filename)
{
#if DEBUG_WLANFRAMESYNC
    if (_q->agc_rx          == NULL ||
        _q->debug_x         == NULL ||
        _q->debug_rssi      == NULL ||
        _q->debug_framesyms == NULL)
    {
        fprintf(stderr,"error: wlanframe_debug_print(), debugging objects don't exist; enable debugging first\n");
        return;
    }

    FILE * fid = fopen(_filename,"w");
    if (!fid) {
        fprintf(stderr,"error: wlanframe_debug_print(), could not open '%s' for writing\n", _filename);
        return;
    }
    fprintf(fid,"%% %s : auto-generated file\n", _filename);

    fprintf(fid,"close all;\n");
    fprintf(fid,"clear all;\n");
    fprintf(fid,"n = %u;\n", DEBUG_WLANFRAMESYNC_BUFFER_LEN);
    unsigned int i;
    float complex * rc;
    float * r;

    fprintf(fid,"x = zeros(1,n);\n");
    windowcf_read(_q->debug_x, &rc);
    for (i=0; i<DEBUG_WLANFRAMESYNC_BUFFER_LEN; i++)
        fprintf(fid,"x(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(rc[i]), cimagf(rc[i]));
    fprintf(fid,"figure;\n");
    fprintf(fid,"plot(0:(n-1),real(x),0:(n-1),imag(x));\n");
    fprintf(fid,"xlabel('sample index');\n");
    fprintf(fid,"ylabel('received signal, x');\n");

    // write agc_rssi
    fprintf(fid,"\n\n");
    fprintf(fid,"agc_rssi = zeros(1,%u);\n", DEBUG_WLANFRAMESYNC_BUFFER_LEN);
    windowf_read(_q->debug_rssi, &r);
    for (i=0; i<DEBUG_WLANFRAMESYNC_BUFFER_LEN; i++)
        fprintf(fid,"agc_rssi(%4u) = %12.4e;\n", i+1, r[i]);
    fprintf(fid,"figure;\n");
    fprintf(fid,"plot(agc_rssi)\n");
    fprintf(fid,"ylabel('RSSI [dB]');\n");
    
    // write frame symbols
    fprintf(fid,"framesyms = zeros(1,n);\n");
    windowcf_read(_q->debug_framesyms, &rc);
    for (i=0; i<DEBUG_WLANFRAMESYNC_BUFFER_LEN; i++)
        fprintf(fid,"framesyms(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(rc[i]), cimagf(rc[i]));
    fprintf(fid,"figure;\n");
    fprintf(fid,"plot(real(framesyms),imag(framesyms),'x','MarkerSize',2);\n");
    fprintf(fid,"axis([-1 1 -1 1]*1.5);\n");
    fprintf(fid,"axis square;\n");
    fprintf(fid,"grid on;\n");
    fprintf(fid,"xlabel('real');\n");
    fprintf(fid,"ylabel('imag');\n");

    // write gain arrays
    fprintf(fid,"\n\n");
    fprintf(fid,"G0a = zeros(1,64);\n");
    fprintf(fid,"G0b = zeros(1,64);\n");
    fprintf(fid,"G1a = zeros(1,64);\n");
    fprintf(fid,"G1b = zeros(1,64);\n");
    fprintf(fid,"G   = zeros(1,64);\n");
    for (i=0; i<64; i++) {
        unsigned int k = (i + 32) % 64;
        fprintf(fid,"G0a(%3u) = %12.8f + j*%12.8f;\n", k+1, crealf(_q->G0a[i]), cimagf(_q->G0a[i]));
        fprintf(fid,"G0b(%3u) = %12.8f + j*%12.8f;\n", k+1, crealf(_q->G0b[i]), cimagf(_q->G0b[i]));
        fprintf(fid,"G1a(%3u) = %12.8f + j*%12.8f;\n", k+1, crealf(_q->G1a[i]), cimagf(_q->G1a[i]));
        fprintf(fid,"G1b(%3u) = %12.8f + j*%12.8f;\n", k+1, crealf(_q->G1b[i]), cimagf(_q->G1b[i]));
        fprintf(fid,"G(%3u)   = %12.8f + j*%12.8f;\n", k+1, crealf(_q->G[i]),   cimagf(_q->G[i]));
    }
    fprintf(fid,"%% apply timing offset (backoff) phase shift\n");
    fprintf(fid,"f = -32:31;\n");
    fprintf(fid,"b = 2;\n");
    fprintf(fid,"G0a = G0a.*exp(j*b*2*pi*f/64);\n");
    fprintf(fid,"G0b = G0b.*exp(j*b*2*pi*f/64);\n");
    fprintf(fid,"G1a = G1a.*exp(j*b*2*pi*f/64);\n");
    fprintf(fid,"G1b = G1b.*exp(j*b*2*pi*f/64);\n");
    fprintf(fid,"G   = G.*exp(j*b*2*pi*f/64);\n");

    fprintf(fid,"figure;\n");
    fprintf(fid,"subplot(2,1,1);\n");
    fprintf(fid,"  plot(f,abs(G1a),'x', f,abs(G1b),'x', f,abs(G),'-k','LineWidth',2);\n");
    fprintf(fid,"  ylabel('G (mag)');\n");
    fprintf(fid,"subplot(2,1,2);\n");
    fprintf(fid,"  plot(f,arg(G1a),'x', f,arg(G1b),'x', f,arg(G),'-k','LineWidth',2);\n");
    fprintf(fid,"  ylabel('G (phase)');\n");

    fclose(fid);
    printf("wlanframesync/debug: results written to '%s'\n", _filename);
#else
    fprintf(stderr,"wlanframesync_debug_print(): compile-time debugging disabled\n");
#endif
}

//.........这里部分代码省略.........
    if(!complx) rcurr  = sf_floatalloc(nzx2);

    cwave  = sf_complexalloc(nk);
    cwavem  = sf_complexalloc(nk);
    wave   = sf_complexalloc2(nk,m2);
    wave2  = sf_complexalloc2(nzx2,m2);

    for (iz=0; iz < nzx2; iz++) {
	curr[iz] = sf_cmplx(0.,0.);
	if(!complx) rcurr[iz]= 0.;
    }

    /* MAIN LOOP */
    for (it=0; it<nt; it++) {
	if(verb) sf_warning("it=%d;",it);

	/* matrix multiplication */
	for (im = 0; im < m2; im++) {
	    for (ix = 0; ix < nx; ix++) {
		for (iz=0; iz < nz; iz++) {
		    i = iz+ix*nz;  /* original grid */
		    j = iz+ix*nz2; /* padded grid */
#ifdef SF_HAS_COMPLEX_H
		    currm[j] = lt[im][i]*curr[j];
#else
		    currm[j] = sf_cmul(lt[im][i], curr[j]);
#endif
		}
	    }
	    cfft2(currm,wave[im]);
	}

/* approach no.2 */
	for (ik = 0; ik < nk; ik++) {
	    c = sf_cmplx(0.,0.);
	    for (im = 0; im < m2; im++) {
#ifdef SF_HAS_COMPLEX_H
		c += wave[im][ik]*rt[ik][im];
#else
		c = sf_cadd(c,sf_cmul(wave[im][ik],rt[ik][im])); /* complex multiplies complex */
#endif
	    }
	    cwave[ik] = c;
	}

	for (im = 0; im < m2; im++) {
	    for (ik = 0; ik < nk; ik++) {
#ifdef SF_HAS_COMPLEX_H
		cwavem[ik] = cwave[ik]*rt[ik][im];
#else
		cwavem[ik] = sf_cmul(cwave[ik],rt[ik][im]); /* complex multiplies complex */
#endif
	    }
	    icfft2(wave2[im],cwavem);
	}

/*
	for (im = 0; im < m2; im++) {
	    for (ik = 0; ik < nk; ik++) {
#ifdef SF_HAS_COMPLEX_H
		cwave[ik] = wave[im][ik]*rt[ik][im]*rt[ik][im];
#else
		cwave[ik] = sf_cmul(sf_cmul(wave[im][ik],rt[ik][im]),rt[ik][im]);
#endif
	    }
	    icfft2(wave2[im],cwave);
	}
*/
	for (ix = 0; ix < nx; ix++) {
	    for (iz=0; iz < nz; iz++) {
		i = iz+ix*nz;  /* original grid */
		j = iz+ix*nz2; /* padded grid */
#ifdef SF_HAS_COMPLEX_H
		c = ww[it] * rr[i]; // source term
#else
		c = sf_crmul(ww[it], rr[i]); // source term
#endif
		for (im = 0; im < m2; im++) {
#ifdef SF_HAS_COMPLEX_H
		    c += lt[im][i]*wave2[im][j];
#else
		    c = sf_cadd(c,sf_cmul(lt[im][i], wave2[im][j]));
#endif
		}

		curr[j] = c;
		if (!complx) rcurr[j] = crealf(c);

	    }

	    /* write wavefield to output */
	    if (complx)
		sf_complexwrite(curr+ix*nz2,nz,Fo);
	    else
		sf_floatwrite(rcurr+ix*nz2,nz,Fo);
	}
    }
    if(verb) sf_warning("."); 
    exit (0);
}

//.........这里部分代码省略.........
    for (i=0; i<num_samples_resampled; i++)
        y[i] += nstd*(randnf() + _Complex_I*randnf());


    // 
    // create and run symbol synchronizer
    //

    symsync_crcf d = symsync_crcf_create_rnyquist(ftype_rx, k, m, beta, num_filters);
    symsync_crcf_set_lf_bw(d,bt);
    symsync_crcf_set_output_rate(d,k_out);

    unsigned int num_samples_sync=0;
    unsigned int nn;
    unsigned int num_symbols_sync = 0;
    float tau_hat[num_samples];
    for (i=ds; i<num_samples_resampled; i++) {
        tau_hat[num_samples_sync] = symsync_crcf_get_tau(d);
        symsync_crcf_execute(d, &y[i], 1, &z[num_samples_sync], &nn);

        // decimate
        unsigned int j;
        for (j=0; j<nn; j++) {
            if ( (num_samples_sync%k_out)==0 )
                sym_out[num_symbols_sync++] = z[num_samples_sync];
            num_samples_sync++;
        }
    }
    symsync_crcf_destroy(d);

    // print last several symbols to screen
    printf("output symbols:\n");
    for (i=num_symbols_sync-10; i<num_symbols_sync; i++)
        printf("  sym_out(%2u) = %8.4f + j*%8.4f;\n", i+1, crealf(sym_out[i]), cimagf(sym_out[i]));

    //
    // export output file
    //

    FILE* fid = fopen(OUTPUT_FILENAME,"w");
    fprintf(fid,"%% %s, auto-generated file\n\n", OUTPUT_FILENAME);
    fprintf(fid,"close all;\nclear all;\n\n");

    fprintf(fid,"k=%u;\n",k);
    fprintf(fid,"m=%u;\n",m);
    fprintf(fid,"beta=%12.8f;\n",beta);
    fprintf(fid,"k_out=%u;\n",k_out);
    fprintf(fid,"num_filters=%u;\n",num_filters);
    fprintf(fid,"num_symbols=%u;\n",num_symbols);

    for (i=0; i<h_len; i++)
        fprintf(fid,"h(%3u) = %12.5f;\n", i+1, h[i]);

    for (i=0; i<num_symbols; i++)
        fprintf(fid,"s(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(s[i]), cimagf(s[i]));

    for (i=0; i<num_samples; i++)
        fprintf(fid,"x(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(x[i]), cimagf(x[i]));
        
    for (i=0; i<num_samples_resampled; i++)
        fprintf(fid,"y(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(y[i]), cimagf(y[i]));
        
    for (i=0; i<num_samples_sync; i++)
        fprintf(fid,"z(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(z[i]), cimagf(z[i]));
        
    for (i=0; i<num_symbols_sync; i++)

int main(int argc, char*argv[]) {
    srand( time(NULL) );
    // parameters
    float phase_offset = M_PI/10;
    float frequency_offset = 0.001f;
    float SNRdB = 30.0f;
    float pll_bandwidth = 0.02f;
    modulation_scheme ms = LIQUID_MODEM_QPSK;
    unsigned int n=256;     // number of iterations

    int dopt;
    while ((dopt = getopt(argc,argv,"uhs:b:n:P:F:m:")) != EOF) {
        switch (dopt) {
        case 'u':
        case 'h':   usage();    return 0;
        case 's':   SNRdB = atof(optarg);           break;
        case 'b':   pll_bandwidth = atof(optarg);   break;
        case 'n':   n = atoi(optarg);               break;
        case 'P':   phase_offset = atof(optarg);    break;
        case 'F':   frequency_offset= atof(optarg); break;
        case 'm':
            ms = liquid_getopt_str2mod(optarg);
            if (ms == LIQUID_MODEM_UNKNOWN) {
                fprintf(stderr,"error: %s, unknown/unsupported modulation scheme \"%s\"\n", argv[0], optarg);
                return 1;
            }
            break;
        default:
            exit(1);
        }
    }
    unsigned int d=n/32;      // print every "d" lines

    FILE * fid = fopen(OUTPUT_FILENAME,"w");
    fprintf(fid, "%% %s : auto-generated file\n", OUTPUT_FILENAME);
    fprintf(fid, "clear all;\n");
    fprintf(fid, "phi=zeros(1,%u);\n",n);
    fprintf(fid, "r=zeros(1,%u);\n",n);

    // objects
    nco_crcf nco_tx = nco_crcf_create(LIQUID_VCO);
    nco_crcf nco_rx = nco_crcf_create(LIQUID_VCO);

    modem mod   = modem_create(ms);
    modem demod = modem_create(ms);

    unsigned int bps = modem_get_bps(mod);

    // initialize objects
    nco_crcf_set_phase(nco_tx, phase_offset);
    nco_crcf_set_frequency(nco_tx, frequency_offset);
    nco_crcf_pll_set_bandwidth(nco_rx, pll_bandwidth);

    float noise_power = powf(10.0f, -SNRdB/20.0f);

    // print parameters
    printf("PLL example :\n");
    printf("modem : %u-%s\n", 1<<bps, modulation_types[ms].name);
    printf("frequency offset: %6.3f, phase offset: %6.3f, SNR: %6.2fdB, pll b/w: %6.3f\n",
            frequency_offset, phase_offset, SNRdB, pll_bandwidth);

    // run loop
    unsigned int i, M=1<<bps, sym_in, sym_out, num_errors=0;
    float phase_error;
    float complex x, r, v, noise;
    for (i=0; i<n; i++) {
        // generate random symbol
        sym_in = rand() % M;

        // modulate
        modem_modulate(mod, sym_in, &x);

        // channel
        //r = nco_crcf_cexpf(nco_tx);
        nco_crcf_mix_up(nco_tx, x, &r);

        // add complex white noise
        crandnf(&noise);
        r += noise * noise_power;

        // 
        //v = nco_crcf_cexpf(nco_rx);
        nco_crcf_mix_down(nco_rx, r, &v);

        // demodulate
        modem_demodulate(demod, v, &sym_out);
        num_errors += count_bit_errors(sym_in, sym_out);

        // error estimation
        //phase_error = cargf(r*conjf(v));
        phase_error = modem_get_demodulator_phase_error(demod);

        // perfect error estimation
        //phase_error = nco_tx->theta - nco_rx->theta;

        // print every line in a format that octave can read
        fprintf(fid, "phi(%u) = %10.6E;\n", i+1, phase_error);
        fprintf(fid, "r(%u) = %10.6E + j*%10.6E;\n",
                i+1, crealf(v), cimagf(v));

//.........这里部分代码省略.........

// ellip_azpkf()
//
// Compute analog zeros, poles, gain of low-pass elliptic
// filter, grouping complex conjugates together. If
// the filter order is odd, the single real pole is at the
// end of the array.
//  _n      :   filter order
//  _ep     :   epsilon_p, related to pass-band ripple
//  _es     :   epsilon_s, related to stop-band ripple
//  _za     :   output analog zeros [length: floor(_n/2)]
//  _pa     :   output analog poles [length: _n]
//  _ka     :   output analog gain
void ellip_azpkf(unsigned int _n,
                 float _ep,
                 float _es,
                 float complex * _za,
                 float complex * _pa,
                 float complex * _ka)
{
    // filter specifications
    float fp = 1.0f / (2.0f * M_PI);    // pass-band cutoff
    float fs = 1.1*fp;                  // stop-band cutoff
    //float Gp = 1/sqrtf(1 + _ep*_ep);    // pass-band gain
    //float Gs = 1/sqrtf(1 + _es*_es);    // stop-band gain

    // number of iterations for elliptic integral
    // approximations
    unsigned int n=7;

    float Wp = 2*M_PI*fp;
    float Ws = 2*M_PI*fs;

    // ripples passband, stopband
    float ep = _ep;
    float es = _es;
#if LIQUID_DEBUG_ELLIP_PRINT
    printf("ep, es      : %12.8f, %12.8f\n", _ep, _es);
#endif

    float k  = Wp/Ws;           // 0.8889f;
    float k1 = ep/es;           // 0.0165f;
#if LIQUID_DEBUG_ELLIP_PRINT
    printf("k           : %12.4e\n", k);
    printf("k1          : %12.4e\n", k1);
#endif

    float K,  Kp;
    float K1, K1p;
    ellipkf(k, n, &K,  &Kp);    // K  = 2.23533416, Kp  = 1.66463780
    ellipkf(k1,n, &K1, &K1p);   // K1 = 1.57090271, K1p = 5.49361753
#if LIQUID_DEBUG_ELLIP_PRINT
    printf("K,  Kp      : %12.8f, %12.8f\n", K,  Kp);
    printf("K1, K1p     : %12.8f, %12.8f\n", K1, K1p);
#endif

    float N      = (float)_n;       // ceilf(Nexact) = 5
#if LIQUID_DEBUG_ELLIP_PRINT
    float Nexact = (K1p/K1)/(Kp/K); // 4.69604063
    printf("N (exact)   : %12.8f\n", Nexact);
    printf("N           : %12.8f\n", N);
#endif

    k = ellipdegf(N,k1,n);      // 0.91427171
#if LIQUID_DEBUG_ELLIP_PRINT
    printf("k           : %12.4e\n", k);
#endif

#if LIQUID_DEBUG_ELLIP_PRINT
    float fs_new = fp/k;        // 4.37506723
    printf("fs_new      : %12.8f\n", fs_new);
#endif

    unsigned int L = (unsigned int)(floorf(N/2.0f)); // 2
    unsigned int r = ((unsigned int)N) % 2;
    float u[L];
    unsigned int i;
    for (i=0; i<L; i++) {
        float t = (float)i + 1.0f;
        u[i] = (2.0f*t - 1.0f)/N;
#if LIQUID_DEBUG_ELLIP_PRINT
        printf("u[%3u]      : %12.8f\n", i, u[i]);
#endif
    }
    float complex zeta[L];
    for (i=0; i<L; i++) {
        zeta[i] = ellip_cdf(u[i],k,n);
#if LIQUID_DEBUG_ELLIP_PRINT
        printf("zeta[%3u]   : %12.8f + j*%12.8f\n", i, crealf(zeta[i]), cimagf(zeta[i]));
#endif
    }

    // compute filter zeros
    float complex za[L];
    for (i=0; i<L; i++) {
        za[i] = _Complex_I * Wp / (k*zeta[i]);
#if LIQUID_DEBUG_ELLIP_PRINT
        printf("za[%3u]     : %12.8f + j*%12.8f\n", i, crealf(za[i]), cimagf(za[i]));
#endif
    }

//.........这里部分代码省略.........

//.........这里部分代码省略.........

   /* check output units */
   expectedUnit = lalDimensionlessUnit;
   expectedUnit.unitNumerator[LALUnitIndexADCCount] = 1;
   expectedUnit.unitNumerator[LALUnitIndexSecond] = 1;
   result = XLALUnitCompare(&expectedUnit, &(goodOutput.sampleUnits));
   if (optVerbose)
   {
     if ( XLALUnitAsString( unitString, LALUnitTextSize, &(goodOutput.sampleUnits)) == NULL ) {
       return SZEROPADANDFFTTESTC_EFLS;
     }
     printf( "Units are \"%s\", ", unitString );
     if ( XLALUnitAsString( unitString, LALUnitTextSize, &expectedUnit) == NULL ) {
       return SZEROPADANDFFTTESTC_EFLS;
     }
     printf( "should be \"%s\"\n", unitString );
   }

   if (result != 0)
   {
     printf("  FAIL: Valid data test #1\n");
     if (optVerbose)
     {
       printf("Exiting with error: %s\n",
              SZEROPADANDFFTTESTC_MSGEFLS);
     }
     return SZEROPADANDFFTTESTC_EFLS;
   }

   /* check output values */
   if (optVerbose)
   {
     printf("hBarTilde(0)=%g + %g i, should be %g\n",
            crealf(goodOutput.data->data[0]), cimagf(goodOutput.data->data[0]),
            crealf(expectedOutputDataData[0]));
   }
   if ( fabsf(crealf(goodOutput.data->data[0]) - crealf(expectedOutputDataData[0]))
        /* / expectedOutputDataData[0].re */> SZEROPADANDFFTTESTC_TOL
        || fabsf(cimagf(goodOutput.data->data[0])) > SZEROPADANDFFTTESTC_TOL )
   {
     printf("  FAIL: Valid data test\n");
     if (optVerbose)
       {
         printf("Exiting with error: %s\n", SZEROPADANDFFTTESTC_MSGEFLS);
       }
     return SZEROPADANDFFTTESTC_EFLS;
   }

   for (i=1; i<SZEROPADANDFFTTESTC_LENGTH; ++i)
   {
     f = i * SZEROPADANDFFTTESTC_DELTAF;
     if (optVerbose)
     {
       printf("hBarTilde(%f Hz)=%g + %g i, should be %g + %g i\n",
              f, crealf(goodOutput.data->data[i]), cimagf(goodOutput.data->data[i]),
              crealf(expectedOutputDataData[i]), cimagf(expectedOutputDataData[i]));
     }
     if (fabsf(crealf(goodOutput.data->data[i]) - crealf(expectedOutputDataData[i]))
         /* / expectedOutputDataData[0].re */> SZEROPADANDFFTTESTC_TOL
         || fabsf(cimagf(goodOutput.data->data[i]) - cimagf(expectedOutputDataData[i]))
         /* / expectedOutputDataData[0].re */> SZEROPADANDFFTTESTC_TOL)
     {
       printf("  FAIL: Valid data test\n");
       if (optVerbose)
       {
         printf("Exiting with error: %s\n", SZEROPADANDFFTTESTC_MSGEFLS);

float (crealf)(float complex z)
{
	return crealf(z);
}

float
my_float_real(float complex c)
{
  return crealf(c);
}

//.........这里部分代码省略.........
						//                    c += curr[j];

						for(im=0; im<m; im++){
#ifdef SF_HAS_COMPLEX_H
							c += lt[im][i]*wave[im][j];
#else
							c += sf_cmul(lt[im][i], wave[im][j]);
#endif
						}
						curr[j]=c;
					}
				}

				if (taper!=0) {
					if (it%taper == 0) {
						cfft2(curr,cwave);
						for (ik = 0; ik < nk; ik++) {
#ifdef SF_HAS_COMPLEX_H
							cwavem[ik] = cwave[ik]*ktp[ik];
#else
							cwavem[ik] = sf_crmul(cwave[ik],ktp[ik]);
#endif
						}
						icfft2(curr,cwavem);
					}
				}

				if(it%scalet==0){
#ifdef _OPENMP
#pragma omp parallel for private(ix, iz)
#endif
					for(ix=0; ix<rnx; ix++){
						for(iz=0; iz<rnz; iz++){
							fwf[wfit][ix][iz]=crealf(curr[(ix+nb)*fnz+(iz+nb)]);
						}
					}
					wfit++;
				}
			} //end of it

			/* check wfnt */
			if(wfit != wfnt) sf_error("At this point, wfit should be equal to wfnt");

			/* backward propagation starts from here... */
#ifdef _OPENMP
#pragma omp parallel for private(iz)
#endif
			for(iz=0; iz<fnzx; iz++){
				curr[iz]=sf_cmplx(0.,0.);
			}

			wfit=wfnt-1;
			for(it=nt-1; it>=0; it--){
				if(verb) sf_warning("Backward propagation it=%d/%d",it+1, nt);
#ifdef _OPENMP
#pragma omp parallel for private(ix)
#endif
				for(ix=0; ix<nr; ix++){
					curr[(nr0+ix*ndr)*fnz+gpz]+=dd[ix][it];
				}

				cfft2(curr, cwave);

				for(im=0; im<m; im++){
#ifdef _OPENMP
#pragma omp parallel for private(ik)

void LALTfrWv (LALStatus *stat, REAL4Vector* sig, TimeFreqRep *tfr, TimeFreqParam *param)
{

  INT4    nf;
  INT4    time;
  INT4    column, row;
  INT4    taumax, tau;
  REAL4Vector  *lacf = NULL;      /* local autocorrelation function */
  COMPLEX8Vector  *vtmp = NULL;
  RealFFTPlan  *plan = NULL;

  INITSTATUS(stat);
  ATTATCHSTATUSPTR (stat);

  /* Make sure the arguments are not NULL: */
  ASSERT (sig, stat, TFR_ENULL, TFR_MSGENULL);
  ASSERT (tfr, stat, TFR_ENULL, TFR_MSGENULL);
  ASSERT (param, stat, TFR_ENULL, TFR_MSGENULL);

  /* Make sure the data pointers are not NULL: */
  ASSERT (sig->data, stat, TFR_ENULL, TFR_MSGENULL);
  ASSERT (tfr->timeInstant, stat, TFR_ENULL, TFR_MSGENULL);
  ASSERT (tfr->freqBin, stat, TFR_ENULL, TFR_MSGENULL);
  ASSERT (tfr->map, stat, TFR_ENULL, TFR_MSGENULL);

  /* Make sure the requested TFR type corresponds to what will be done */
  ASSERT (tfr->type == WignerVille , stat, TFR_ENAME, TFR_MSGENAME);
  ASSERT (param->type == WignerVille , stat, TFR_ENAME, TFR_MSGENAME);

  /* Make sure the number of freq bins is a positive number: */
  nf = tfr->fRow;
  ASSERT (nf > 0 , stat, TFR_EFROW, TFR_MSGEFROW);

  /* Make sure the number of freq bins is a power of 2: */
  while(!(nf & 1))
      nf = nf>>1;

  ASSERT (nf == 1, stat, TFR_EFROW, TFR_MSGEFROW);

  /* Make sure the timeInstant indicates existing time instants */
  for (column=0 ; column<tfr->tCol ; column++)
    {
      if ((tfr->timeInstant[column] < 0) || (tfr->timeInstant[column] > (INT4)(sig->length-1)))
	{
	  ASSERT (tfr->timeInstant[column] > 0, stat, TFR_EBADT, TFR_MSGEBADT);
	  ASSERT (tfr->timeInstant[column] < (INT4)sig->length, stat, TFR_EBADT, TFR_MSGEBADT);
	}
    }
    TRY(LALCreateForwardRealFFTPlan(stat->statusPtr, &plan,(UINT4)tfr->fRow,0),stat);

  TRY(LALSCreateVector(stat->statusPtr, &lacf, tfr->fRow), stat);
  TRY(LALCCreateVector(stat->statusPtr, &vtmp, tfr->fRow/2+1), stat);

  for (column = 0; column < tfr->tCol; column++)
    {

      for (row = 0; row < tfr->fRow; row++)
	lacf->data[row] = 0.0;

      time = tfr->timeInstant[column];
      taumax = MIN (time, (INT4)(sig->length -1 - time));
      taumax = MIN (taumax, (tfr->fRow / 2 - 1));

      for (tau = -taumax; tau <= taumax; tau++)
	{
	  row = (tfr->fRow+tau)%tfr->fRow;
	  lacf->data[row] =   sig->data[time + tau]*sig->data[time - tau];
        }

      tau=tfr->fRow/2;
      if ((time<=(INT4)sig->length-tau-1)&(time>=tau))
	lacf->data[tau] =  sig->data[time+tau]*sig->data[time-tau];

      LALForwardRealFFT(stat->statusPtr, vtmp, lacf, plan);

      for (row = 0; row < (tfr->fRow/2+1); row++)
	tfr->map[column][row] = crealf(vtmp->data[row]);

    }
  /* Reflecting frequency halfing in WV distrob so multiply by 1/2 */
  for (row = 0; row < (tfr->fRow/2+1) ; row++)
    tfr->freqBin[row] = (REAL4) row / (2 *tfr->fRow);

  TRY(LALSDestroyVector(stat->statusPtr, &lacf), stat);
  TRY(LALCDestroyVector(stat->statusPtr, &vtmp), stat);

  TRY(LALDestroyRealFFTPlan(stat->statusPtr, &plan), stat);

  DETATCHSTATUSPTR (stat);

  (void)param;

  /* normal exit */
  RETURN (stat);
}

//.........这里部分代码省略.........
    us = sf_input("us");

    if (!sf_histint(us,"n3",&ns)) sf_error("No ns=.");
    if (!sf_histint(us,"n4",&nw)) sf_error("No nw=.");
    if (!sf_histfloat(us,"d4",&dw)) sf_error("No dw=.");
    if (!sf_histfloat(us,"o4",&ow)) sf_error("No ow=.");

    /* read receiver wavefield */	
    if (NULL == sf_getstring("ur"))
	sf_error("Need receiver wavefield ur=");
    ur = sf_input("ur");    
    
    /* read wavelet */
    if (NULL == sf_getstring("wvlt"))
	sf_error("Need wvlt=");
    wvlt = sf_input("wvlt");
    
    f = sf_complexalloc2(nw,ns);
    sf_complexread(f[0],nw*ns,wvlt);
    sf_fileclose(wvlt);

    /* read list */
    if (NULL == sf_getstring("list"))
	sf_error("Need list=");
    list = sf_input("list");

    pp = sf_intalloc2(2,ns);
    sf_intread(pp[0],2*ns,list);
    sf_fileclose(list);

    /* allocate memory */
    swave = sf_complexalloc3(n1,n2,ns);
    rwave = sf_complexalloc3(n1,n2,ns);
    stemp = sf_complexalloc2(ns,ns);
    rtemp = sf_complexalloc2(ns,ns);

    ahesss = sf_floatalloc2(n1*n2,ns);
    ahessr = sf_floatalloc2(n1*n2,ns);

    ahess = sf_floatalloc(n1*n2);

    /* loop over frequency */
    for (iw=0; iw < nw; iw++) {
	omega = (double) 2.*SF_PI*(ow+iw*dw);	

	/* read wavefields */
	sf_complexread(swave[0][0],n1*n2*ns,us);
	sf_complexread(rwave[0][0],n1*n2*ns,ur);

#ifdef _OPENMP
#pragma omp parallel num_threads(uts) private(is,ip,i)
#endif
	{
#ifdef _OPENMP
#pragma omp for
#endif
	    for (is=0; is < ns; is++) {
		for (ip=0; ip < ns; ip++) {
		    /* temps */
		    stemp[is][ip] = -omega*omega/conjf(f[ip][iw])
			*rwave[is][pp[ip][1]][pp[ip][0]];
		    
		    /* tempr */
		    rtemp[is][ip] = -omega*omega/conjf(f[ip][iw])
			*conjf(swave[is][pp[ip][1]][pp[ip][0]]);
		}
	    }

	    /* loop over model */
#ifdef _OPENMP
#pragma omp for
#endif
	    for (i=0; i < n1*n2; i++) {
		for (is=0; is < ns; is++) {
		    for (ip=0; ip < ns; ip++) {
			ahesss[ip][i] += crealf(
			    conjf(swave[ip][0][i]*swave[is][0][i])*stemp[is][ip]);
			ahessr[ip][i] += crealf(
			    conjf(swave[ip][0][i])*rwave[is][0][i]*rtemp[is][ip]);
		    }
		}
	    }
	}	
    }

    /* assemble */
#ifdef _OPENMP
#pragma omp parallel for num_threads(uts) private(i,ip)
#endif
    for (i=0; i < n1*n2; i++) {
	for (ip=0; ip < ns; ip++) {
	    ahess[i] += powf(ahesss[ip][i]+ahessr[ip][i],2.);
	}
    }

    /* output hessian */
    sf_floatwrite(ahess,n1*n2,out);

    exit(0);
}

int main() {
    int             ftype   = LIQUID_IIRDES_BUTTER; // filter prototype
    unsigned int    n       =     7;                // Hilbert filter order
    float           Ap      =  0.1f;                // pass-band ripple [dB]
    float           As      = 60.0f;                // stop-band attenuation [dB]
    float           fc      = 0.123456;             // signal center frequency
    unsigned int num_input_samples=180;             // number of samples

    // derived values
    unsigned int num_total_samples = num_input_samples + 50; // allow for filter settling

    // create Hilbert transform objects
    iirhilbf qi = iirhilbf_create(ftype,n,Ap,As);    // interpolator
    iirhilbf qd = iirhilbf_create(ftype,n,Ap,As);    // decimator
    iirhilbf_print(qi);

    // data arrays
    float complex x[  num_total_samples];   // complex input
    float         y[2*num_total_samples];   // real output
    float complex z[  num_total_samples];   // complex output

    // initialize input array
    unsigned int i;
    for (i=0; i<num_total_samples; i++) {
        x[i]  = cexpf(_Complex_I*2*M_PI*fc*i) +
                cexpf(_Complex_I*2*M_PI*fc*i*1.3f)*0.1f;
        x[i] *= (i < num_input_samples) ? 1.855f*hamming(i,num_input_samples) : 0.0f;
    }

    // execute interpolator (complex to real conversion)
    iirhilbf_interp_execute_block(qi, x, num_total_samples, y);
    
    // execute decimator (real to complex conversion)
    iirhilbf_decim_execute_block(qd, y, num_total_samples, z);

    // destroy Hilbert transform object
    iirhilbf_destroy(qi);
    iirhilbf_destroy(qd);

    // 
    // export results to file
    //
    FILE*fid = fopen(OUTPUT_FILENAME,"w");
    fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
    fprintf(fid,"clear all;\n");
    fprintf(fid,"close all;\n");
    fprintf(fid,"num_input_samples=%u;\n", num_input_samples);
    fprintf(fid,"num_total_samples=%u;\n", num_total_samples);
    fprintf(fid,"tx = 0:(num_total_samples-1);\n");
    fprintf(fid,"ty = [0:(2*num_total_samples-1)]/2;\n");
    fprintf(fid,"tz = tx;\n");

    for (i=0; i<num_total_samples; i++) {
        // print results
        fprintf(fid,"x(%3u) = %12.4e + %12.4ej;\n",   i+1, crealf(x[i]), cimagf(x[i]));
        fprintf(fid,"y(%3u) = %12.4e;\n",           2*i+1, y[2*i+0]);
        fprintf(fid,"y(%3u) = %12.4e;\n",           2*i+2, y[2*i+1]);
        fprintf(fid,"z(%3u) = %12.4e + %12.4ej;\n",   i+1, crealf(z[i]), cimagf(z[i]));
    }

    fprintf(fid,"figure;\n");
    fprintf(fid,"subplot(3,1,1);\n");
    fprintf(fid,"  plot(tx,real(x),'Color',[0.00 0.25 0.50],'LineWidth',1.3,...\n");
    fprintf(fid,"       tx,imag(x),'Color',[0.00 0.50 0.25],'LineWidth',1.3);\n");
    fprintf(fid,"  legend('real','imag','location','northeast');\n");
    fprintf(fid,"  ylabel('transformed/complex');\n");
    fprintf(fid,"  axis([0 num_total_samples -2 2]);\n");
    fprintf(fid,"  grid on;\n");
    fprintf(fid,"subplot(3,1,2);\n");
    fprintf(fid,"  plot(ty,y,'Color',[0.00 0.25 0.50],'LineWidth',1.3);\n");
    fprintf(fid,"  ylabel('original/real');\n");
    fprintf(fid,"  axis([0 num_total_samples -2 2]);\n");
    fprintf(fid,"  grid on;\n");
    fprintf(fid,"subplot(3,1,3);\n");
    fprintf(fid,"  plot(tz,real(z),'Color',[0.00 0.25 0.50],'LineWidth',1.3,...\n");
    fprintf(fid,"       tz,imag(z),'Color',[0.00 0.50 0.25],'LineWidth',1.3);\n");
    fprintf(fid,"  legend('real','imag','location','northeast');\n");
    fprintf(fid,"  ylabel('transformed/complex');\n");
    fprintf(fid,"  axis([0 num_total_samples -2 2]);\n");
    fprintf(fid,"  grid on;\n");

    // plot results
    fprintf(fid,"nfft=4096;\n");
    fprintf(fid,"%% compute normalized windowing functions\n");
    fprintf(fid,"X=20*log10(abs(fftshift(fft(x/num_input_samples,nfft))));\n");
    fprintf(fid,"Y=20*log10(abs(fftshift(fft(y/num_input_samples,nfft))));\n");
    fprintf(fid,"Z=20*log10(abs(fftshift(fft(z/num_input_samples,nfft))));\n");
    fprintf(fid,"f =[0:(nfft-1)]/nfft-0.5;\n");
    fprintf(fid,"figure; plot(f+0.5,X,'LineWidth',1,'Color',[0.50 0.50 0.50],...\n");
    fprintf(fid,"             f*2,  Y,'LineWidth',2,'Color',[0.00 0.50 0.25],...\n");
    fprintf(fid,"             f+0.5,Z,'LineWidth',1,'Color',[0.00 0.25 0.50]);\n");
    fprintf(fid,"grid on;\n");
    fprintf(fid,"axis([-1.0 1.0 -80 20]);\n");
    fprintf(fid,"xlabel('normalized frequency');\n");
    fprintf(fid,"ylabel('PSD [dB]');\n");
    fprintf(fid,"legend('original/cplx','transformed/real','regenerated/cplx','location','northeast');");

    fclose(fid);
    printf("results written to %s\n", OUTPUT_FILENAME);

//.........这里部分代码省略.........

// receive the 'SIGNAL' field
void wlanframesync_execute_rxsignal(wlanframesync _q)
{
    _q->timer++;
    if (_q->timer < 80)
        return;

    // reset timer
    _q->timer = 0;

    // run fft
    float complex * rc;
    windowcf_read(_q->input_buffer, &rc);
    memmove(_q->x, &rc[16-2], 64*sizeof(float complex));

    // compute fft, storing result into _q->X
    FFT_EXECUTE(_q->fft);
  
    // recover symbol, correcting for gain, pilot phase, etc.
    wlanframesync_rxsymbol(_q);
    
    // demodulate, decode, ...
    memset(_q->signal_int, 0x00, 6*sizeof(unsigned char));

    _q->signal_int[0] |= crealf(_q->X[38]) > 0.0f ? 0x80 : 0x00;
    _q->signal_int[0] |= crealf(_q->X[39]) > 0.0f ? 0x40 : 0x00;
    _q->signal_int[0] |= crealf(_q->X[40]) > 0.0f ? 0x20 : 0x00;
    _q->signal_int[0] |= crealf(_q->X[41]) > 0.0f ? 0x10 : 0x00;
    _q->signal_int[0] |= crealf(_q->X[42]) > 0.0f ? 0x08 : 0x00;
    //  43 : pilot
    _q->signal_int[0] |= crealf(_q->X[44]) > 0.0f ? 0x04 : 0x00;
    _q->signal_int[0] |= crealf(_q->X[45]) > 0.0f ? 0x02 : 0x00;
    _q->signal_int[0] |= crealf(_q->X[46]) > 0.0f ? 0x01 : 0x00;
    _q->signal_int[1] |= crealf(_q->X[47]) > 0.0f ? 0x80 : 0x00;
    _q->signal_int[1] |= crealf(_q->X[48]) > 0.0f ? 0x40 : 0x00;
    _q->signal_int[1] |= crealf(_q->X[49]) > 0.0f ? 0x20 : 0x00;
    _q->signal_int[1] |= crealf(_q->X[50]) > 0.0f ? 0x10 : 0x00;
    _q->signal_int[1] |= crealf(_q->X[51]) > 0.0f ? 0x08 : 0x00;
    _q->signal_int[1] |= crealf(_q->X[52]) > 0.0f ? 0x04 : 0x00;
    _q->signal_int[1] |= crealf(_q->X[53]) > 0.0f ? 0x02 : 0x00;
    _q->signal_int[1] |= crealf(_q->X[54]) > 0.0f ? 0x01 : 0x00;
    _q->signal_int[2] |= crealf(_q->X[55]) > 0.0f ? 0x80 : 0x00;
    _q->signal_int[2] |= crealf(_q->X[56]) > 0.0f ? 0x40 : 0x00;
    //  57 : pilot
    _q->signal_int[2] |= crealf(_q->X[58]) > 0.0f ? 0x20 : 0x00;
    _q->signal_int[2] |= crealf(_q->X[59]) > 0.0f ? 0x10 : 0x00;
    _q->signal_int[2] |= crealf(_q->X[60]) > 0.0f ? 0x08 : 0x00;
    _q->signal_int[2] |= crealf(_q->X[61]) > 0.0f ? 0x04 : 0x00;
    _q->signal_int[2] |= crealf(_q->X[62]) > 0.0f ? 0x02 : 0x00;
    _q->signal_int[2] |= crealf(_q->X[63]) > 0.0f ? 0x01 : 0x00;
    //   0 : NULL
    _q->signal_int[3] |= crealf(_q->X[ 1]) > 0.0f ? 0x80 : 0x00;
    _q->signal_int[3] |= crealf(_q->X[ 2]) > 0.0f ? 0x40 : 0x00;
    _q->signal_int[3] |= crealf(_q->X[ 3]) > 0.0f ? 0x20 : 0x00;
    _q->signal_int[3] |= crealf(_q->X[ 4]) > 0.0f ? 0x10 : 0x00;
    _q->signal_int[3] |= crealf(_q->X[ 5]) > 0.0f ? 0x08 : 0x00;
    _q->signal_int[3] |= crealf(_q->X[ 6]) > 0.0f ? 0x04 : 0x00;
    //   7 : pilot
    _q->signal_int[3] |= crealf(_q->X[ 8]) > 0.0f ? 0x02 : 0x00;
    _q->signal_int[3] |= crealf(_q->X[ 9]) > 0.0f ? 0x01 : 0x00;
    _q->signal_int[4] |= crealf(_q->X[10]) > 0.0f ? 0x80 : 0x00;
    _q->signal_int[4] |= crealf(_q->X[11]) > 0.0f ? 0x40 : 0x00;
    _q->signal_int[4] |= crealf(_q->X[12]) > 0.0f ? 0x20 : 0x00;
    _q->signal_int[4] |= crealf(_q->X[13]) > 0.0f ? 0x10 : 0x00;
    _q->signal_int[4] |= crealf(_q->X[14]) > 0.0f ? 0x08 : 0x00;
    _q->signal_int[4] |= crealf(_q->X[15]) > 0.0f ? 0x04 : 0x00;
    _q->signal_int[4] |= crealf(_q->X[16]) > 0.0f ? 0x02 : 0x00;
    _q->signal_int[4] |= crealf(_q->X[17]) > 0.0f ? 0x01 : 0x00;
    _q->signal_int[5] |= crealf(_q->X[18]) > 0.0f ? 0x80 : 0x00;
    _q->signal_int[5] |= crealf(_q->X[19]) > 0.0f ? 0x40 : 0x00;
    _q->signal_int[5] |= crealf(_q->X[20]) > 0.0f ? 0x20 : 0x00;
    //  21 : pilot
    _q->signal_int[5] |= crealf(_q->X[22]) > 0.0f ? 0x10 : 0x00;
    _q->signal_int[5] |= crealf(_q->X[23]) > 0.0f ? 0x08 : 0x00;
    _q->signal_int[5] |= crealf(_q->X[24]) > 0.0f ? 0x04 : 0x00;
    _q->signal_int[5] |= crealf(_q->X[25]) > 0.0f ? 0x02 : 0x00;
    _q->signal_int[5] |= crealf(_q->X[26]) > 0.0f ? 0x01 : 0x00;

    // decode SIGNAL field
    wlanframesync_decode_signal(_q);

    // validate proper decoding
    if (!_q->signal_valid) {
        // reset synchronizer and return
        wlanframesync_reset(_q);
        return;
    }

    // set state
    _q->state = WLANFRAMESYNC_STATE_RXDATA;
}

// frame detection
void wlanframesync_execute_seekplcp(wlanframesync _q)
{
    _q->timer++;

    // TODO : only check every 100 - 150 (decimates/reduced complexity)
    if (_q->timer < 64)
        return;

    // reset timer
    _q->timer = 0;

    // read contents of input buffer
    float complex * rc;
    windowcf_read(_q->input_buffer, &rc);
    
    // estimate gain
    // TODO : use gain from result of FFT
    unsigned int i;
    float g = 0.0f;
    for (i=16; i<80; i+=4) {
        // compute |rc[i]|^2 efficiently
        g += crealf(rc[i  ])*crealf(rc[i  ]) + cimagf(rc[i  ])*cimagf(rc[i  ]);
        g += crealf(rc[i+1])*crealf(rc[i+1]) + cimagf(rc[i+1])*cimagf(rc[i+1]);
        g += crealf(rc[i+2])*crealf(rc[i+2]) + cimagf(rc[i+2])*cimagf(rc[i+2]);
        g += crealf(rc[i+3])*crealf(rc[i+3]) + cimagf(rc[i+3])*cimagf(rc[i+3]);
    }
    g = 64.0f / (g + 1e-12f);
    
    // save gain (permits dynamic invocation of get_rssi() method)
    _q->g0 = g;

    // estimate S0 gain
    wlanframesync_estimate_gain_S0(_q, &rc[16], _q->G0a);
    
    // compute S0 metrics
    float complex s_hat;
    wlanframesync_S0_metrics(_q, _q->G0a, &s_hat);
    s_hat *= g;

    float tau_hat  = cargf(s_hat) * (float)(16.0f) / (2*M_PI);
#if DEBUG_WLANFRAMESYNC_PRINT
    printf(" - gain=%12.3f, rssi=%8.2f dB, s_hat=%12.4f <%12.8f>, tau_hat=%8.3f\n",
            sqrt(g),
            -10*log10(g),
            cabsf(s_hat), cargf(s_hat),
            tau_hat);
#endif

    // 
    if (cabsf(s_hat) > WLANFRAMESYNC_S0A_ABS_THRESH) {

        int dt = (int)roundf(tau_hat);
        // set timer appropriately...
        _q->timer = (16 + dt) % 16;
        //_q->timer += 32; // add delay to help ensure good S0 estimate (multiple of 16)
        _q->state = WLANFRAMESYNC_STATE_RXSHORT0;

#if DEBUG_WLANFRAMESYNC_PRINT
        printf("********** frame detected! ************\n");
        printf("    s_hat   :   %12.8f <%12.8f>\n", cabsf(s_hat), cargf(s_hat));
        printf("  tau_hat   :   %12.8f\n", tau_hat);
        printf("    dt      :   %12d\n", dt);
        printf("    timer   :   %12u\n", _q->timer);
#endif
    }

}

int main(int argc, char* argv[])
{
    sf_file Fw_s=NULL, Fw_r=NULL;
    sf_file Fi=NULL, Fm=NULL, Fr=NULL;
    sf_axis ag,at,aw,ar;
    int ig,it,iw,ir;
    int ng,nt,nw,nr;

    int method;
    bool  verb;
    bool   adj;

    sf_complex **dat_s=NULL, *wfl_s=NULL;
    sf_complex **dat_r=NULL, *wfl_r=NULL;
    float      **img=NULL;
    float      **aa=NULL,**bb=NULL,**mm=NULL;

    sf_complex **ab=NULL;
    float      **a0=NULL,**b0=NULL;

    float w,ws,wr,w0,dw;
    char *met="";

    sf_init(argc,argv);

    if(! sf_getbool("verb", &verb))     verb=false;
    if(! sf_getint("method",&method)) method=0;    /* extrapolation method */
    if(! sf_getbool("adj",  &adj))       adj=false;/* y=modeling; n=migration */
    Fm = sf_input("abm");
    Fr = sf_input("abr");

    at=sf_iaxa(Fm,2); sf_setlabel(at,"t"); /* 'extrapolation axis' */
    ar=sf_iaxa(Fr,1); sf_setlabel(ar,"r"); /* a,b reference */
    if(method==0) sf_setn(ar,1); /* pure F-D */
    nr=sf_n(ar); 

    Fw_s = sf_input ( "in");
    if (SF_COMPLEX !=sf_gettype(Fw_s)) sf_error("Need complex source");

    /* 'position axis' (could be angle) */
    ag = sf_iaxa(Fw_s,1); ng=sf_n(ag); sf_setlabel(ag,"g"); 

    /* 'extrapolation axis' (could be time) */
    at = sf_iaxa(Fw_s,2); nt=sf_n(at); sf_setlabel(at,"t");
    aw = sf_iaxa(Fw_s,3); sf_setlabel(aw,"w"); /* frequency */

    if(adj) { /* modeling */
	Fw_r = sf_output("out"); 
	sf_settype(Fw_r,SF_COMPLEX);

	Fi   = sf_input ("img");
	if (SF_FLOAT !=sf_gettype(Fi)) sf_error("Need float image");

	sf_oaxa(Fw_r,ag,1);
	sf_oaxa(Fw_r,at,2);
	sf_oaxa(Fw_r,aw,3);
    } else {  /* migration */
	Fw_r = sf_input ("rwf");
	if (SF_COMPLEX !=sf_gettype(Fw_r)) sf_error("Need complex data");

	Fi   = sf_output("out"); 
	sf_settype(Fi,SF_FLOAT);
		
	sf_oaxa(Fi,ag,1);
	sf_oaxa(Fi,at,2);
	sf_putint(Fi,"n3",1);
    }

    img   = sf_floatalloc2  (ng,nt);
    dat_s = sf_complexalloc2(ng,nt);
    dat_r = sf_complexalloc2(ng,nt);
    wfl_s = sf_complexalloc (ng);
    wfl_r = sf_complexalloc (ng);

    if(verb) {
	sf_raxa(ag);
	sf_raxa(at);
	sf_raxa(aw);
	sf_raxa(ar);
    }

    /* read ABM */
    aa = sf_floatalloc2  (ng,nt);
    bb = sf_floatalloc2  (ng,nt);
    mm = sf_floatalloc2  (ng,nt);

    sf_floatread(aa[0],ng*nt,Fm); /* a coef */
    sf_floatread(bb[0],ng*nt,Fm); /* b coef */
    sf_floatread(mm[0],ng*nt,Fm); /* mask */

    /* read ABr */
    ab = sf_complexalloc2(nr,nt);
    a0 = sf_floatalloc2  (nr,nt);
    b0 = sf_floatalloc2  (nr,nt);

    sf_complexread(ab[0],nr*nt,Fr);
    for(it=0;it<nt;it++) {
	for(ir=0;ir<nr;ir++) {
	    a0[it][ir] = crealf(ab[it][ir]);
	    b0[it][ir] = cimagf(ab[it][ir]);
//.........这里部分代码省略.........