xmippDoc/html/xmipp__fft_8h.html

#include <data/matrix1d.h>
#include <data/xmipp_fft.h>
double discreteTransform(double w, int N1) {
   if (w==0) return 2*N1+1;
   else return sin(w*(N1+0.5))/sin(0.5*w);
}
double continuousTransform(double W, double T1) {
   if (W==0) return 2*T1;
   else return 2*sin(W*T1)/W;
}
int main() {
    try {
      MultidimArray<double> x(65);
  x.setXmippOrigin();
  double T=0.5;
  double T1=6;
  int N1=(int)CEIL(T1/T);
  // Fill x with a pulse from -N1 to N1 (-T1 to T1 in continuous)
  FOR_ALL_ELEMENTS_IN_ARRAY1D(x)
     if (ABS(i)<=N1) x(i)=1;
  // Compute the Fourier transform
  MultidimArray< std::complex<double> > X;
  MultidimArray<double> Xmag;
  FourierTransform(x,X);
  FFT_magnitude(X,Xmag);
  // Compute the frequency axes
  MultidimArray<double> contfreq(XSIZE(X)), digfreq(XSIZE(X));
         FOR_ALL_ELEMENTS_IN_ARRAY1D(X)
             FFT_IDX2DIGFREQ(i,XSIZE(X),digfreq(i));
  digfreq*=2*PI;
  contfreq=digfreq/T;
  // Show all Fourier transforms
  FOR_ALL_ELEMENTS_IN_ARRAY1D(X) {
      if (digfreq(i)>=0)
                std::cout << digfreq(i) << " " << contfreq(i) << " "
            << XSIZE(X)*Xmag(i) << " "
     << ABS(discreteTransform(digfreq(i),N1)) << " "
     << ABS(continuousTransform(contfreq(i),T1)/T)
     << std::endl;
  }
    } catch (XmippError XE) {
        std::cout << XE << std::endl;
    }
    return 0;
}
Functions
Format conversions
void	Whole2Half (const MultidimArray< std::complex< double > > &in, MultidimArray< std::complex< double > > &out)

void	Half2Whole (const MultidimArray< std::complex< double > > &in, MultidimArray< std::complex< double > > &out, size_t oridim)

void	Complex2RealImag (const MultidimArray< std::complex< double > > &in, MultidimArray< double > &real, MultidimArray< double > &imag)

void	RealImag2Complex (const MultidimArray< double > &real, const MultidimArray< double > &imag, MultidimArray< std::complex< double > > &out)

Fourier Transforms
The theoretical relationship between the Fourier transform of a discrete signal and the Fourier transform of the continuous signal is X(e^jw)=1/TX_c(jw/T) Xmipp is not computing X(e^jw) but samples from it so that X(e^jw)=NX_XMIPP[k] where N is the length of the signal being transformed and X_XMIPP[k] is the k-th sample. The following program illustrates how the continuous, discrete and Xmipp Fourier transform relate #include <data/matrix1d.h> #include <data/xmipp_fft.h> double discreteTransform(double w, int N1) { if (w==0) return 2N1+1; else return sin(w(N1+0.5))/sin(0.5w); } double continuousTransform(double W, double T1) { if (W==0) return 2T1; else return 2sin(WT1)/W; } int main() { try { MultidimArray<double> x(65); x.setXmippOrigin(); double T=0.5; double T1=6; int N1=(int)CEIL(T1/T); // Fill x with a pulse from -N1 to N1 (-T1 to T1 in continuous) FOR_ALL_ELEMENTS_IN_ARRAY1D(x) if (ABS(i)<=N1) x(i)=1; // Compute the Fourier transform MultidimArray< std::complex<double> > X; MultidimArray<double> Xmag; FourierTransform(x,X); FFT_magnitude(X,Xmag); // Compute the frequency axes MultidimArray<double> contfreq(XSIZE(X)), digfreq(XSIZE(X)); FOR_ALL_ELEMENTS_IN_ARRAY1D(X) FFT_IDX2DIGFREQ(i,XSIZE(X),digfreq(i)); digfreq=2PI; contfreq=digfreq/T; // Show all Fourier transforms FOR_ALL_ELEMENTS_IN_ARRAY1D(X) { if (digfreq(i)>=0) std::cout << digfreq(i) << " " << contfreq(i) << " " << XSIZE(X)*Xmag(i) << " " << ABS(discreteTransform(digfreq(i),N1)) << " " << ABS(continuousTransform(contfreq(i),T1)/T) << std::endl; } } catch (XmippError XE) { std::cout << XE << std::endl; } return 0; }
void	FourierTransform (const MultidimArray< double > &in, MultidimArray< std::complex< double > > &out)

void	InverseFourierTransform (const MultidimArray< std::complex< double > > &in, MultidimArray< double > &out)

void	FourierTransformHalf (const MultidimArray< double > &in, MultidimArray< std::complex< double > > &out)

void	InverseFourierTransformHalf (const MultidimArray< std::complex< double > > &in, MultidimArray< double > &out, int oridim)
Index <--> Frequency, Continuous <--> Discrete
#define	FFT_IDX2DIGFREQ(idx, size, freq)

#define	FFT_IDX2DIGFREQ_DOUBLE(idx, size, freq)

#define	FFT_IDX2DIGFREQ_FAST(idx, size, size_2, isize, freq) freq = ( ((idx) <= (size_2)) ? (idx) : -(size) + (idx) ) * (isize);

#define	DIGFREQ2FFT_IDX(freq, size, idx)

#define	DIGFREQ2FFT_IDX_DOUBLE(freq, size, idx)

template<typename T >
void	FFT_idx2digfreq (T &v, const Matrix1D< int > &idx, Matrix1D< double > &freq)

template<typename T >
void	digfreq2FFT_idx (T &v, const Matrix1D< double > &freq, Matrix1D< int > &idx)

void	digfreq2contfreq (const Matrix1D< double > &digfreq, Matrix1D< double > &contfreq, double pixel_size)

void	contfreq2digfreq (const Matrix1D< double > &contfreq, Matrix1D< double > &digfreq, double pixel_size)
Operations with the Fourier Transforms
#define	SWAP_ARRAY(a, b, n) memcpy(buffer, a, n); memcpy(a, b, n); memcpy(b, buffer, n);

void	centerFFT2 (MultidimArray< double > &v)

template<typename T >
void	CenterFFT (MultidimArray< T > &v, bool forward)

void	ShiftFFT (MultidimArray< std::complex< double > > &v, double xshift)

void	ShiftFFT (MultidimArray< std::complex< double > > &v, double xshift, double yshift)

void	ShiftFFT (MultidimArray< std::complex< double > > &v, double xshift, double yshift, double zshift)

void	CenterOriginFFT (MultidimArray< std::complex< double > > &v, bool forward)

template<typename T >
void	xmipp2PSD (const MultidimArray< T > &input, MultidimArray< T > &output, bool takeLog=true)

void	xmipp2CTF (const MultidimArray< double > &input, MultidimArray< double > &output)