本文整理汇总了C++中sf_crmul函数的典型用法代码示例。如果您正苦于以下问题：C++ sf_crmul函数的具体用法？C++ sf_crmul怎么用？C++ sf_crmul使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了sf_crmul函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: cbanded_solve
void cbanded_solve (sf_complex *b)
/*< multiply by inverse (in place) >*/
{
int k, m;
sf_complex t;
for (k = 1; k < band; k++) {
t = b[k];
for (m = 1; m <= k; m++) {
#ifdef SF_HAS_COMPLEX_H
t -= o[m-1][k-m] * b[k-m];
#else
t = sf_csub(t,sf_cmul(o[m-1][k-m],b[k-m]));
#endif
}
b[k] = t;
}
for (k = band; k < n; k++) {
t = b[k];
for (m = 1; m <= band; m++) {
#ifdef SF_HAS_COMPLEX_H
t -= o[m-1][k-m] * b[k-m];
#else
t = sf_csub(t,sf_cmul(o[m-1][k-m],b[k-m]));
#endif
}
b[k] = t;
}
for (k = n-1; k >= n - band; k--) {
#ifdef SF_HAS_COMPLEX_H
t = b[k]/d[k];
#else
t = sf_crmul(b[k],1./d[k]);
#endif
for (m = 0; m < n - k - 1; m++) {
#ifdef SF_HAS_COMPLEX_H
t -= conjf(o[m][k]) * b[k+m+1];
#else
t = sf_csub(t,sf_cmul(conjf(o[m][k]),b[k+m+1]));
#endif
}
b[k] = t;
}
for (k = n - band - 1; k >= 0; k--) {
#ifdef SF_HAS_COMPLEX_H
t = b[k]/d[k];
#else
t = sf_crmul(b[k],1./d[k]);
#endif
for (m = 0; m < band; m++) {
#ifdef SF_HAS_COMPLEX_H
t -= conjf(o[m][k]) * b[k+m+1];
#else
t = sf_csub(t,sf_cmul(conjf(o[m][k]),b[k+m+1]));
#endif
}
b[k] = t;
}
}示例2: ctoeplitz_solve
void ctoeplitz_solve (const sf_complex *r /* top row of the matrix */,
sf_complex *f /* inverted in place */)
/*< apply the solver >*/
{
int i,j;
sf_complex e,c,w, bot;
float v;
v=crealf(r[0]);
#ifdef SF_HAS_COMPLEX_H
f[0] /= v;
#else
f[0] = sf_crmul(f[0],1./v);
#endif
for (j=1; j < n; j++) {
e = cdprod(j,a,r);
#ifdef SF_HAS_COMPLEX_H
c = -e/v;
#else
c = sf_crmul(e,-1./v);
#endif
v += crealf(c)*crealf(e) + cimagf(c)*cimagf(e);
for (i=1; i <= j/2; i++) {
#ifdef SF_HAS_COMPLEX_H
bot = a[j-i] + c*conjf(a[i]);
a[i] += c*conjf(a[j-i]);
#else
bot = sf_cadd(a[j-i],sf_cmul(c,conjf(a[i])));
a[i] = sf_cadd(a[i],sf_cmul(c,conjf(a[j-i])));
#endif
a[j-i] = bot;
}
a[j] = c;
w = cdprod(j,f,r);
#ifdef SF_HAS_COMPLEX_H
c = (f[j]-w)/v;
#else
c = sf_crmul(sf_csub(f[j],w),1./v);
#endif
for (i=0; i < j; i++) {
#ifdef SF_HAS_COMPLEX_H
f[i] += c*conjf(a[j-i]);
#else
f[i] = sf_cadd(f[i],sf_cmul(c,conjf(a[j-i])));
#endif
}
f[j] = c;
}
}示例3: xkolmog
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]);
}
}示例4: sf_ompfft3a1
/*------------------------------------------------------------*/
void sf_ompfft3a1(bool inv /* inverse/forward flag */,
kiss_fft_cpx ***pp /* [n1][n2][n3] */,
ompfft3d fft,
int ompith)
/*< apply FFT on axis 1 >*/
{
int i1, i2, i3;
if (inv) {
/* IFT 1 */
for (i3=0; i3 < fft->n3; i3++) {
for(i2=0; i2 < fft->n2; i2++) {
#ifdef _OPENMP
#pragma omp critical
#endif
kiss_fft(fft->invs[ompith],pp[i3][i2],pp[i3][i2]);
}
}
/* scaling */
for (i3=0; i3 < fft->n3; i3++) {
for (i2=0; i2 < fft->n2; i2++) {
for(i1=0; i1 < fft->n1; i1++) {
pp[i3][i2][i1] = sf_crmul(pp[i3][i2][i1],fft->scale);
}
}
}
} else {
/* scaling */
for (i3=0; i3 < fft->n3; i3++) {
for (i2=0; i2 < fft->n2; i2++) {
for(i1=0; i1 < fft->n1; i1++) {
pp[i3][i2][i1] = sf_crmul(pp[i3][i2][i1],fft->scale);
}
}
}
/* FFT 1 */
for (i3=0; i3 < fft->n3; i3++) {
for(i2=0; i2 < fft->n2; i2++) {
#ifdef _OPENMP
#pragma omp critical
#endif
kiss_fft(fft->forw[ompith],pp[i3][i2],pp[i3][i2]);
}
}
}
}示例5: fft2
/*------------------------------------------------------------*/
void fft2(bool inv /* inverse/forward flag */,
kiss_fft_cpx **pp /* [1...n2][1...n1] */)
/*< Apply 2-D FFT >*/
{
int i1,i2;
if (inv) {
/* IFT 1 */
for(i2=0; i2 < n2; i2++) {
kiss_fft(invs1,pp[i2],pp[i2]);
}
/* IFT 2 */
for(i1=0; i1 < n1; i1++) {
kiss_fft_stride(invs2,pp[0]+i1,trace2,n1);
for(i2=0; i2<n2; i2++) {
pp[i2][i1] = trace2[i2];
}
}
/* scaling */
for (i2=0; i2<n2; i2++) {
for(i1=0; i1 < n1; i1++) {
pp[i2][i1] = sf_crmul(pp[i2][i1],fftscale);
}
}
} else {
/* scaling */
for (i2=0; i2<n2; i2++) {
for(i1=0; i1 < n1; i1++) {
pp[i2][i1] = sf_crmul(pp[i2][i1],fftscale);
}
}
/* FFT 2 */
for(i1=0; i1 < n1; i1++) {
kiss_fft_stride(forw2,pp[0]+i1,trace2,n1);
for(i2=0; i2<n2; i2++) {
pp[i2][i1] = trace2[i2];
}
}
/* FFT 1 */
for(i2=0; i2 < n2; i2++) {
kiss_fft(forw1,pp[i2],pp[i2]);
}
}
}示例6: sf_fft3a3
/*------------------------------------------------------------*/
void sf_fft3a3(bool inv /* inverse/forward flag */,
kiss_fft_cpx ***pp /* [n1][n2][n3] */,
sf_fft3d fft)
/*< apply FFT on axis 3 >*/
{
int i1, i2, i3;
if (inv) {
/* IFT 3 */
for (i2=0; i2 < fft->n2; i2++) {
for(i1=0; i1 < fft->n1; i1++) {
kiss_fft_stride(fft->invs,pp[0][0]+i1+i2*fft->n1,fft->trace,fft->n1*fft->n2);
for(i3=0; i3 < fft->n3; i3++) {
pp[i3][i2][i1] = fft->trace[i3];
}
}
}
/* scaling */
for (i3=0; i3 < fft->n3; i3++) {
for (i2=0; i2 < fft->n2; i2++) {
for(i1=0; i1 < fft->n1; i1++) {
pp[i3][i2][i1] = sf_crmul(pp[i3][i2][i1],fft->scale);
}
}
}
} else {
/* scaling */
for (i3=0; i3 < fft->n3; i3++) {
for (i2=0; i2 < fft->n2; i2++) {
for(i1=0; i1 < fft->n1; i1++) {
pp[i3][i2][i1] = sf_crmul(pp[i3][i2][i1],fft->scale);
}
}
}
/* FFT 3 */
for (i2=0; i2 < fft->n2; i2++) {
for(i1=0; i1 < fft->n1; i1++) {
kiss_fft_stride(fft->forw,pp[0][0]+i1+i2*fft->n1,fft->trace,fft->n1*fft->n2);
for(i3=0; i3 < fft->n3; i3++) {
pp[i3][i2][i1] = fft->trace[i3];
}
}
}
}
}示例7: BornScatteredField
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;
} 示例8: icfft2
void icfft2(sf_complex *out /* [n1*n2] */,
sf_complex *inp /* [nk*n2] */)
/*< 2-D inverse FFT >*/
{
int i1, i2;
#pragma omp parallel for private(i2,i1) default(shared)
for (i2=0; i2<local_n0; i2++) {
for (i1=0; i1<nk; i1++) {
dd[i2*nk+i1]=inp[i2*nk+i1];
}
}
fftwf_execute(icfg);
/* FFT centering and normalization*/
#pragma omp parallel for private(i2,i1) default(shared)
for (i2=0; i2<local_n0; i2++) {
for (i1=0; i1<n1; i1++) {
#ifdef SF_HAS_COMPLEX_H
out[i2*n1+i1] = ((((i2+local_0_start)%2==0)==(i1%2==0))? wt:-wt) * cc[i2*n1+i1];
#else
out[i2*n1+i1] = sf_crmul(cc[i2*n1+i1],((((i2+local_0_start)%2==0)==(i1%2==0))? wt:-wt));
#endif
}
}
}示例9: explsourcet
void explsourcet(sf_complex *curr/*@[email protected]*/,
sf_complex *vwavlet,
int vit, int vsx, int vsz,
int nx2, int nz2,
srcpar vps/*decay parameters*/)
/*< explosive source >*/
{
float phi = 0.0;
int cent = (int)vps->range/2;
int ix, iz;
if (vps->decay ==1){
#ifdef _OPENMP
#pragma omp parallel for private(ix,iz,phi)
#endif
for (ix=0; ix<2*cent; ix++)
for (iz=0; iz<2*cent; iz++) {
phi = exp( -1*vps->alpha*vps->alpha*((ix-cent)*(ix-cent)+(iz-cent)*(iz-cent)) );
#ifdef SF_HAS_COMPLEX_H
curr[(vsx-cent+ix)*nz2+(vsz-cent+iz)] += vwavlet[vit]*phi;
#else
curr[(vsx-cent+ix)*nz2+(vsz-cent+iz)] += sf_crmul(vwavlet[vit],phi);
#endif
}
} else {
curr[vsx*nz2+vsz] += vwavlet[vit];
}
}示例10: icfft2
void icfft2(sf_complex *out /* [n1*n2] */,
sf_complex *inp /* [nkk*n2] */)
/*< 2-D inverse FFT >*/
{
int i1, i2;
#ifdef SF_HAS_FFTW
fftwf_execute(icfg);
#else
for (i1=0; i1 < nkk; i1++) {
kiss_fft_stride(icfg2,(kiss_fft_cpx *) (inp+i1),ctrace2,nkk);
for (i2=0; i2<n2; i2++) {
temp[i2][i1] = ctrace2[i2];
}
}
for (i2=0; i2 < n2; i2++) {
kiss_fft_stride(icfg1,temp[i2],(kiss_fft_cpx *) cc[i2],1);
}
#endif
/* FFT centering and normalization*/
for (i2=0; i2<n2; i2++) {
for (i1=0; i1<n1; i1++) {
#ifdef SF_HAS_COMPLEX_H
out[i2*n1+i1] = (((i2%2==0)==(i1%2==0))? wt:-wt) * cc[i2][i1];
#else
out[i2*n1+i1] = sf_crmul(cc[i2][i1],(((i2%2==0)==(i1%2==0))? wt:-wt));
#endif
}
}
}示例11: triple
static void triple (int o, int d, int nx, int nb, const sf_complex* x, sf_complex* tmp, bool box)
{
int i;
float wt;
sf_complex xi;
for (i=0; i < nx + 2*nb; i++) {
tmp[i] = sf_cmplx(0.,0.);
}
if (box) {
wt = 1.0/(2*nb-1);
for (i=0; i < nx; i++) {
#ifdef SF_HAS_COMPLEX_H
xi = wt*x[o+i*d];
tmp[i+1] += xi;
tmp[i+2*nb] -= xi;
#else
xi = sf_crmul(x[o+i*d],wt);
tmp[i+1] = sf_cadd(tmp[i+1],xi);
tmp[i+2*nb] = sf_cadd(tmp[i+2*nb],sf_cneg(xi));
#endif
}
} else {
wt = 1.0/(nb*nb);
for (i=0; i < nx; i++) {
#ifdef SF_HAS_COMPLEX_H
xi = wt*x[o+i*d];
tmp[i] -= xi;
tmp[i+nb] += 2*xi;
tmp[i+2*nb] -= xi;
#else
xi = sf_crmul(x[o+i*d],wt);
tmp[i] = sf_cadd(tmp[i],sf_cneg(xi));
tmp[i+nb] = sf_cadd(tmp[i+nb],sf_crmul(xi,2.));
tmp[i+2*nb] = sf_cadd(tmp[i+2*nb],sf_cneg(xi));
#endif
}
}
}示例12: sf_fft3a1
/*------------------------------------------------------------*/
void sf_fft3a1(bool inv /* inverse/forward flag */,
kiss_fft_cpx ***pp /* [n1][n2][n3] */,
sf_fft3d fft)
/*< apply FFT on axis 1 >*/
{
int i1, i2, i3;
if (inv) {
/* IFT 1 */
for (i3=0; i3 < fft->n3; i3++) {
for(i2=0; i2 < fft->n2; i2++) {
kiss_fft(fft->invs,pp[i3][i2],pp[i3][i2]);
}
}
/* scaling */
for (i3=0; i3 < fft->n3; i3++) {
for (i2=0; i2 < fft->n2; i2++) {
for(i1=0; i1 < fft->n1; i1++) {
pp[i3][i2][i1] = sf_crmul(pp[i3][i2][i1],fft->scale);
}
}
}
} else {
/* scaling */
for (i3=0; i3 < fft->n3; i3++) {
for (i2=0; i2 < fft->n2; i2++) {
for(i1=0; i1 < fft->n1; i1++) {
pp[i3][i2][i1] = sf_crmul(pp[i3][i2][i1],fft->scale);
}
}
}
/* FFT 1 */
for (i3=0; i3 < fft->n3; i3++) {
for(i2=0; i2 < fft->n2; i2++) {
kiss_fft(fft->forw,pp[i3][i2],pp[i3][i2]);
}
}
}
}示例13: init_wave
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
}
}示例14: rweone_fft
void rweone_fft(
bool inv,
kiss_fft_cpx *d)
/*< apply FFT >*/
{
int ig;
if(inv) {
kiss_fft(invs,d,d);
for(ig=0;ig<ag.n;ig++) { d[ig] = sf_crmul(d[ig],ffts); }
} else {
for(ig=0;ig<ag.n;ig++) { d[ig] = sf_crmul(d[ig],ffts); }
kiss_fft(forw,d,d);
}
}示例15: RytovSensitivity
void RytovSensitivity(float xx,float xy,float xz,float *output)
/*< Rytov sensitivity >*/
{
int iw;
double omega,scale;
sf_complex U,dU;
sf_complex val;
*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;
#ifdef SF_HAS_COMPLEX_H
/* background field */
U=Green(rx,ry,rz,sx,sy,sz,omega)*scale;
/* scattered field */
val=Green(xx,xy,xz,sx,sy,sz,omega)*scale;
val*=Green(rx,ry,rz,xx,xy,xz,omega);
dU=val*(-omega*omega*dv);
val=U*omega;
*output -= scale*cimagf(dU/U);
#else
/* background field */
U=sf_crmul(Green(rx,ry,rz,sx,sy,sz,omega),scale);
/* scattered field */
val=sf_crmul(Green(xx,xy,xz,sx,sy,sz,omega),scale);
val=sf_cmul(val,Green(rx,ry,rz,xx,xy,xz,omega));
dU=sf_crmul(val,-omega*omega*dv);
val=sf_crmul(U,omega);
*output -= scale*cimagf(sf_cdiv(dU,U));
#endif
}
return;
}