/** * Computes a complex-to-complex dimension-2 fast Fourier transform. * Transforms a 3-D input array cx[n3][nfft][2*n1] of n3*nfft*n1 complex * numbers to a 3-D output array cy[n3][nfft][2*n1] of n3*nfft*n1 complex * numbers. * @param sign the sign (1 or -1) of the exponent used in the FFT. * @param n1 the 1st dimension of arrays. * @param n3 the 3rd dimension of arrays. * @param cx the input array. * @param cy the output array. */ public void complexToComplex2( int sign, int n1, int n3, float[][][] cx, float[][][] cy) { checkSign(sign); checkArray(2*n1,_nfft,n3,cx,"cx"); checkArray(2*n1,_nfft,n3,cy,"cy"); for (int i3=0; i3<n3; ++i3) complexToComplex2(sign,n1,cx[i3],cy[i3]); }
/** * Applies a forward space-to-frequency transform of a 2D array. * @param f the array to be transformed, a sampled function of space. * @return the transformed array, a sampled function of frequency. */ public float[][] applyForward(float[][] f) { ensureSamplingX2(f); float[][] fpad = pad(f); int nx2 = _sx2.getCount(); if (_complex) { _fft1c.complexToComplex1(_sign1,nx2,fpad,fpad); _fft2.complexToComplex2(_sign2,_nfft1,fpad,fpad); } else { _fft1r.realToComplex1(_sign1,nx2,fpad,fpad); _fft2.complexToComplex2(_sign2,_nfft1/2+1,fpad,fpad); } phase(fpad); center(fpad); return fpad; }
/** * Computes a complex-to-complex dimension-3 fast Fourier transform. * Transforms a 3-D input array cx[nfft][n2][2*n1] of nfft*n2*n1 complex * numbers to a 3-D output array cy[nfft][n2][2*n1] of nfft*n2*n1 complex * numbers. * @param sign the sign (1 or -1) of the exponent used in the FFT. * @param n1 the 1st dimension of arrays. * @param n2 the 2nd dimension of arrays. * @param cx the input array. * @param cy the output array. */ public void complexToComplex3( int sign, int n1, int n2, float[][][] cx, float[][][] cy) { checkSign(sign); checkArray(2*n1,n2,_nfft,cx,"cx"); checkArray(2*n1,n2,_nfft,cy,"cy"); float[][] cxi2 = new float[_nfft][]; float[][] cyi2 = new float[_nfft][]; for (int i2=0; i2<n2; ++i2) { for (int i3=0; i3<_nfft; ++i3) { cxi2[i3] = cx[i3][i2]; cyi2[i3] = cy[i3][i2]; } complexToComplex2(sign,n1,cxi2,cyi2); } }
extrapolate(xfft); _fft1.realToComplex1(-1,_nfft2,xfft,xfft); _fft2.complexToComplex2(-1,_nfft1/2+1,xfft,xfft); int nk1 = _nfft1/2+1; int nk2 = _nfft2; _fft2.complexToComplex2(1,_nfft1/2+1,xfft,xfft); _fft1.complexToReal1(1,_nfft2,xfft,xfft); copy(nx1,nx2,xfft,y);
/** * Applies a forward space-to-frequency transform of a 3D array. * @param f the array to be transformed, a sampled function of space. * @return the transformed array, a sampled function of frequency. */ public float[][][] applyForward(float[][][] f) { ensureSamplingX3(f); float[][][] fpad = pad(f); int nx2 = _sx2.getCount(); int nx3 = _sx3.getCount(); if (_complex) { _fft1c.complexToComplex1(_sign1,nx2,nx3,fpad,fpad); _fft2.complexToComplex2(_sign2,_nfft1,nx3,fpad,fpad); _fft3.complexToComplex3(_sign3,_nfft1,_nfft2,fpad,fpad); } else { _fft1r.realToComplex1(_sign1,nx2,nx3,fpad,fpad); _fft2.complexToComplex2(_sign2,_nfft1/2+1,nx3,fpad,fpad); _fft3.complexToComplex3(_sign3,_nfft1/2+1,_nfft2,fpad,fpad); } phase(fpad); center(fpad); return fpad; }
private void updateFfts(int nx1, int nx2) { if (_fft2==null || _h2fft==null || _nx1!=nx1 || _nx2!=nx2) { _nx1 = nx1; _nx2 = nx2; _nx3 = 0; _nfft1 = FftReal.nfftFast(_nx1+_nh1); _nfft2 = FftComplex.nfftFast(_nx2+_nh2); _fft1 = new FftReal(_nfft1); _fft2 = new FftComplex(_nfft2); _fft3 = null; _h1fft = null; _h2fft = new float[_nfft2][_nfft1+2]; _h3fft = null; float scale = 1.0f/(float)_nfft1/(float)_nfft2; for (int ih2=0; ih2<_nh2; ++ih2) { int jh2 = ih2-_kh2; if (jh2<0) jh2 += _nfft2; for (int ih1=0; ih1<_nh1; ++ih1) { int jh1 = ih1-_kh1; if (jh1<0) jh1 += _nfft1; _h2fft[jh2][jh1] = scale*_h2[ih2][ih1]; } } _fft1.realToComplex1(-1,_nfft2,_h2fft,_h2fft); _fft2.complexToComplex2(-1,_nfft1/2+1,_h2fft,_h2fft); } }
extrapolate(xfft); _fft1.realToComplex1(-1,_nfft2,_nfft3,xfft,xfft); _fft2.complexToComplex2(-1,_nfft1/2+1,_nfft3,xfft,xfft); _fft3.complexToComplex3(-1,_nfft1/2+1,_nfft2,xfft,xfft); int nk1 = _nfft1/2+1; _fft2.complexToComplex2(1,_nfft1/2+1,_nfft3,xfft,xfft); _fft1.complexToReal1(1,_nfft2,_nfft3,xfft,xfft); copy(nx1,nx2,nx3,xfft,y);
/** * Applies an inverse frequency-to-space transform of a 2D array. * @param g the array to be transformed, a sampled function of frequency. * @return the transformed array, a sampled function of space. */ public float[][] applyInverse(float[][] g) { ensureSamplingK2(g); float[][] gpad = (_overwrite)?g:copy(g); int nx1 = _sx1.getCount(); int nx2 = _sx2.getCount(); uncenter(gpad); unphase(gpad); if (_complex) { _fft2.complexToComplex2(-_sign2,_nfft1,gpad,gpad); _fft2.scale(_nfft1,nx2,gpad); _fft1c.complexToComplex1(-_sign1,nx2,gpad,gpad); _fft1c.scale(nx1,nx2,gpad); return ccopy(nx1,nx2,gpad); } else { _fft2.complexToComplex2(-_sign2,_nfft1/2+1,gpad,gpad); _fft2.scale(_nfft1/2+1,nx2,gpad); _fft1r.complexToReal1(-_sign1,nx2,gpad,gpad); _fft1r.scale(nx1,nx2,gpad); return copy(nx1,nx2,gpad); } }
_fft2.complexToComplex2(-1,_nfft1/2+1,_nfft3,_h3fft,_h3fft); _fft3.complexToComplex3(-1,_nfft1/2+1,_nfft2,_h3fft,_h3fft);
/** * Applies inverse 2D Fourier transform to an input wavenumber-domain image. * The output space-domain image is written to the appropriate part of the * steerable pyramid array. * @param lev level number. * @param dir basis filter orientation index for the input image. * @param cf input image in wavenumber domain (complex array). * @param spyr input/output 2D steerable pyramid. */ private void ftInverse(int lev,int dir, float[][] cf,float spyr[][][][]) { FftReal fft1; FftComplex fft2; int nf2 = cf.length; int nf1c = cf[0].length/2; int nf1 = (nf1c-1)*2; int mpad = round(20.0f/(1.0f+(float)lev)); int lfactor = (int)pow(2.0,(double)lev); int nl2 = (n2-1)/lfactor+1; int nl1 = (n1-1)/lfactor+1; fft1 = new FftReal(nf1); fft2 = new FftComplex(nf2); fft2.complexToComplex2(-1,nf1c,cf,cf); flipSign(2,cf); fft2.scale(nf1c,nf2,cf); fft1.complexToReal1(-1,nf2,cf,cf); fft1.scale(nf1,nf2,cf); copy(nl1,nl2,mpad,mpad,1,1,cf,0,0,1,1,spyr[lev][dir]); }
_fft3.complexToComplex3(-_sign3,_nfft1,_nfft2,gpad,gpad); _fft3.scale(_nfft1,_nfft2,nx3,gpad); _fft2.complexToComplex2(-_sign2,_nfft1,nx3,gpad,gpad); _fft2.scale(_nfft1,nx2,nx3,gpad); _fft1c.complexToComplex1(-_sign1,nx2,nx3,gpad,gpad); _fft3.complexToComplex3(-_sign3,_nfft1/2+1,_nfft2,gpad,gpad); _fft3.scale(_nfft1/2+1,_nfft2,nx3,gpad); _fft2.complexToComplex2(-_sign2,_nfft1/2+1,nx3,gpad,gpad); _fft2.scale(_nfft1/2+1,nx2,nx3,gpad); _fft1r.complexToReal1(-_sign1,nx2,nx3,gpad,gpad);
/** * Applies forward 2D Fourier transform. * @param level level number. * @param x input 2D image. * @return a 2D complex array, which is the forward 2D Fourier transform of * the input image. */ private float[][] ftForward(int level, float[][] x) { FftReal fft1; FftComplex fft2; int ny2 = x.length; int ny1 = x[0].length; int mpad = round(20.0f/(1.0f+(float)level)); int lfactor = (int)pow(2.0,(double)level); int nl2 = (n2-1)/lfactor+1; int nl1 = (n1-1)/lfactor+1; int nf1 = FftReal.nfftSmall(nl1+mpad*2); int nf1c = nf1/2+1; int nf2 = FftComplex.nfftSmall(nl2+mpad*2); float[][] xr = zerofloat(nf1,nf2); copy(ny1,ny2,0,0,x,mpad,mpad,xr); float[][] cx = czerofloat(nf1c,nf2); fft1 = new FftReal(nf1); fft2 = new FftComplex(nf2); fft1.realToComplex1(1,nf2,xr,cx); flipSign(2,cx); fft2.complexToComplex2(1,nf1c,cx,cx); return cx; }
public void test12Random() { int n1max = 26; int n2max = 26; for (int n2=2; n2<n2max; ++n2) { int n2fft = FftComplex.nfftSmall(n2); FftComplex fft2 = new FftComplex(n2fft); for (int n1=2; n1<n1max; ++n1) { int n1fft = FftReal.nfftSmall(n1); FftReal fft1 = new FftReal(n1fft); int nw = n1fft/2+1; float[][] rr = randfloat(n1fft,n2); float[][] rx = copy(rr); float[][] cy = czerofloat(nw,n2fft); fft1.realToComplex1( 1,n2,rx,cy); fft2.complexToComplex2( 1,nw,cy,cy); fft2.complexToComplex2(-1,nw,cy,cy); fft2.scale(nw,n2,cy); fft1.complexToReal1(-1,n2,cy,rx); fft1.scale(n1,n2,rx); assertRealEqual(n1,n2,rr,rx); } } }
public void test2Random() { int n1max = 26; int n2max = 26; for (int n2=2; n2<n2max; ++n2) { int n2fft = FftComplex.nfftSmall(n2); FftComplex fft2 = new FftComplex(n2fft); for (int n1=2; n1<n1max; ++n1) { int n1fft = FftComplex.nfftSmall(n1); FftComplex fft1 = new FftComplex(n1fft); float[][] cr = crandfloat(n1fft,n2fft); float[][] cx = ccopy(cr); float[][] cy = czerofloat(n1fft,n2fft); fft1.complexToComplex1( 1,n2fft,cx,cy); fft2.complexToComplex2( 1,n1fft,cy,cy); fft1.complexToComplex1(-1,n2fft,cy,cx); fft2.complexToComplex2(-1,n1fft,cx,cx); fft1.scale(n1fft,n2fft,cx); fft2.scale(n1fft,n2fft,cx); assertEqual(cr,cx); } } }
flipSign(3, cf); fft3.scale(nf1c,nf2,nf3,cf); fft2.complexToComplex2(-1,nf1c,nf3,cf,cf); flipSign(2, cf); fft2.scale(nf1c,nf2,nf3,cf);
public void test3Random() { int n1 = 11; int n2 = 12; int n3 = 13; int n1fft = FftComplex.nfftSmall(n1); int n2fft = FftComplex.nfftSmall(n2); int n3fft = FftComplex.nfftSmall(n3); FftComplex fft1 = new FftComplex(n1fft); FftComplex fft2 = new FftComplex(n2fft); FftComplex fft3 = new FftComplex(n3fft); float[][][] cr = crandfloat(n1fft,n2fft,n3fft); float[][][] cx = ccopy(cr); fft1.complexToComplex1( 1,n2fft,n3fft,cx,cx); fft2.complexToComplex2( 1,n1fft,n3fft,cx,cx); fft3.complexToComplex3( 1,n1fft,n2fft,cx,cx); fft1.complexToComplex1(-1,n2fft,n3fft,cx,cx); fft2.complexToComplex2(-1,n1fft,n3fft,cx,cx); fft3.complexToComplex3(-1,n1fft,n2fft,cx,cx); fft1.scale(n1fft,n2fft,n3fft,cx); fft2.scale(n1fft,n2fft,n3fft,cx); fft3.scale(n1fft,n2fft,n3fft,cx); assertEqual(cr,cx); }
public void test2() { int n1max = 26; int n2max = 26; for (int n2=2; n2<n2max; ++n2) { int n2fft = FftComplex.nfftSmall(n2); FftComplex fft2 = new FftComplex(n2fft); for (int n1=2; n1<n1max; ++n1) { int n1fft = FftComplex.nfftSmall(n1); FftComplex fft1 = new FftComplex(n1fft); float[][] c1 = czerofloat(n1fft,n2fft); c1[1][2] = 1.0f; float[][] cx = ccopy(c1); fft1.complexToComplex1(1,n2fft,cx,cx); fft2.complexToComplex2(1,n1fft,cx,cx); float ra = 0.0f; float rb1 = 2.0f*FLT_PI/(float)n1fft; float rb2 = 2.0f*FLT_PI/(float)n2fft; float[][] amp = fillfloat(1.0f,n1fft,n2fft); float[][] phs = rampfloat(ra,rb1,rb2,n1fft,n2fft); float[][] cc = polar(amp,phs); assertEqual(cc,cx); fft1.complexToComplex1(-1,n2fft,cx,cx); fft2.complexToComplex2(-1,n1fft,cx,cx); fft1.scale(n1fft,n2fft,cx); fft2.scale(n1fft,n2fft,cx); assertEqual(c1,cx); } } }
fft1.realToComplex1(1,nf2,nf3,xr,cx); flipSign(2, cx); fft2.complexToComplex2(1,nf1c,nf3,cx,cx); flipSign(3, cx); fft3.complexToComplex3(1,nf1c,nf2,cx,cx);
float[][] cx = rx; fft1.realToComplex1(1,n2,rx,cx); fft2.complexToComplex2(1,nw,cx,cx); float ra = 0.0f; float rb1 = 2.0f*FLT_PI/(float)n1fft; float[][] cc = polar(amp,phs); assertComplexEqual(nw,n2fft,cc,cx); fft2.complexToComplex2(-1,nw,cx,cx); fft2.scale(nw,n2,cx); fft1.complexToReal1(-1,n2,cx,rx);