C# FourierDirection类代码示例

		public static void DFT(Complex[] data, FourierDirection direction) {
			int length = data.Length;
			Complex[] complexArray = new Complex[length];
			for (int i = 0; i < length; i++) {
				complexArray[i] = Complex.Zero;
				double num2 = (((((double)-((int)direction)) * 2.0) * 3.1415926535897931) * i) / ((double)length);
				for (int j = 0; j < length; j++) {
					double num3 = Math.Cos(j * num2);
					double num4 = Math.Sin(j * num2);
					complexArray[i].Re += (float)((data[j].Re * num3) - (data[j].Im * num4));
					complexArray[i].Im += (float)((data[j].Re * num4) + (data[j].Im * num3));
				}
			}
			if (direction == FourierDirection.Forward) {
				for (int k = 0; k < length; k++) {
					data[k].Re = complexArray[k].Re / ((float)length);
					data[k].Im = complexArray[k].Im / ((float)length);
				}
			} else {
				for (int m = 0; m < length; m++) {
					data[m].Re = complexArray[m].Re;
					data[m].Im = complexArray[m].Im;
				}
			}
		}

 /// <summary>
 /// Performs a out-of-place fourier transformation. The original values are kept.
 /// </summary>
 /// <param name="inputarr">The data to transform.</param>
 /// <param name="direction">Specify forward or reverse transformation here.</param>
 /// <param name="outputarr">. On output, contains the fourier transformed data.</param>
 public void Transform(double[] inputarr, FourierDirection direction, double[] outputarr)
 {
   if(inputarr.Length!=_numberOfData)
     throw new ArgumentException(string.Format("Length of array inputarr ({0}) is different from the length specified at construction ({1})",inputarr.Length,_numberOfData),"inputarr");
   if(outputarr.Length!=_numberOfData)
     throw new ArgumentException(string.Format("Length of array outputarr ({0}) is different from the length specified at construction ({1})",outputarr.Length,_numberOfData),"outputarr");
   
   Array.Copy(inputarr,0,outputarr,0,inputarr.Length);
   Transform(outputarr,direction);
 }

		public static void DFT2(Complex[,] data, FourierDirection direction) {
			double num3;
			double num4;
			double num5;
			int length = data.GetLength(0);
			int num2 = data.GetLength(1);
			Complex[] complexArray = new Complex[Math.Max(length, num2)];
			for (int i = 0; i < length; i++) {
				for (int k = 0; k < num2; k++) {
					complexArray[k] = Complex.Zero;
					num3 = (((((double)-((int)direction)) * 2.0) * 3.1415926535897931) * k) / ((double)num2);
					for (int m = 0; m < num2; m++) {
						num4 = Math.Cos(m * num3);
						num5 = Math.Sin(m * num3);
						complexArray[k].Re += (float)((data[i, m].Re * num4) - (data[i, m].Im * num5));
						complexArray[k].Im += (float)((data[i, m].Re * num5) + (data[i, m].Im * num4));
					}
				}
				if (direction == FourierDirection.Forward) {
					for (int n = 0; n < num2; n++) {
						data[i, n].Re = complexArray[n].Re / ((float)num2);
						data[i, n].Im = complexArray[n].Im / ((float)num2);
					}
				} else {
					for (int num10 = 0; num10 < num2; num10++) {
						data[i, num10].Re = complexArray[num10].Re;
						data[i, num10].Im = complexArray[num10].Im;
					}
				}
			}
			for (int j = 0; j < num2; j++) {
				for (int num12 = 0; num12 < length; num12++) {
					complexArray[num12] = Complex.Zero;
					num3 = (((((double)-((int)direction)) * 2.0) * 3.1415926535897931) * num12) / ((double)length);
					for (int num13 = 0; num13 < length; num13++) {
						num4 = Math.Cos(num13 * num3);
						num5 = Math.Sin(num13 * num3);
						complexArray[num12].Re += (float)((data[num13, j].Re * num4) - (data[num13, j].Im * num5));
						complexArray[num12].Im += (float)((data[num13, j].Re * num5) + (data[num13, j].Im * num4));
					}
				}
				if (direction == FourierDirection.Forward) {
					for (int num14 = 0; num14 < length; num14++) {
						data[num14, j].Re = complexArray[num14].Re / ((float)length);
						data[num14, j].Im = complexArray[num14].Im / ((float)length);
					}
				} else {
					for (int num15 = 0; num15 < length; num15++) {
						data[num15, j].Re = complexArray[num15].Re;
						data[num15, j].Im = complexArray[num15].Im;
					}
				}
			}
		}

        //======================================================================================
        //======================================================================================
        /// <summary>
        ///   Compute a 1D fast Fourier transform of a dataset of complex numbers (as pairs of float's).
        /// </summary>
        /// <param name = "data"></param>
        /// <param name = "length"></param>
        /// <param name = "direction"></param>
        public static void FFT(float[] data, int length, FourierDirection direction)
        {
            Debug.Assert(data != null);
            Debug.Assert(data.Length >= length*2);
            Debug.Assert(IsPowerOf2(length));

            SyncLookupTableLength(length);

            int ln = Log2(length);

            // reorder array
            ReorderArray(data);

            // successive doubling
            int N = 1;
            int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
            for (int level = 1; level <= ln; level++)
            {
                int M = N;
                N <<= 1;

                float[] uRLookup = _uRLookupF[level, signIndex];
                float[] uILookup = _uILookupF[level, signIndex];

                for (int j = 0; j < M; j++)
                {
                    float uR = uRLookup[j];
                    float uI = uILookup[j];

                    for (int evenT = j; evenT < length; evenT += N)
                    {
                        int even = evenT << 1;
                        int odd = (evenT + M) << 1;

                        float r = data[odd];
                        float i = data[odd + 1];

                        float odduR = r*uR - i*uI;
                        float odduI = r*uI + i*uR;

                        r = data[even];
                        i = data[even + 1];

                        data[even] = r + odduR;
                        data[even + 1] = i + odduI;

                        data[odd] = r - odduR;
                        data[odd + 1] = i - odduI;
                    }
                }
            }
        }

 /// <summary>
 ///     Fourier transform
 /// </summary>
 /// <param name="array">Array</param>
 /// <param name="direction">Fourier direction</param>
 public static void Fourier(Array array, FourierDirection direction)
 {
     var handle = GCHandle.Alloc(array, GCHandleType.Pinned);
     FftwLock.WaitOne();
     var plan = dft(array.Rank, Enumerable.Range(0, array.Rank).Select(array.GetLength).ToArray(),
         handle.AddrOfPinnedObject(), handle.AddrOfPinnedObject(),
         (fftw_direction) direction,
         fftw_flags.Estimate);
     execute(plan);
     destroy_plan(plan);
     FftwLock.ReleaseMutex();
     handle.Free();
 }

			public ChirpNativeFFTStorage(int msize, int arrsize, FourierDirection direction)
			{
				_msize = msize;
				_arrSize = arrsize;
				_direction = direction;

				_xjfj_real = new double[msize];
				_xjfj_imag = new double[msize];
				_fserp_real = new double[msize];
				_fserp_imag = new double[msize];
				_resarray_real = new double[msize];
				_resarray_imag = new double[msize];
				_chirpfactors_real = new double[arrsize];
				_chirpfactors_imag = new double[arrsize];

				PreCompute_ChirpFactors(); // Precompute the factors for An: Exp(sign * I * Pi * i^2/N)
				Precompute_Fouriertransformed_ChirpFactorsConjugate(); // Pre-compute fserp using Pre-computed chirp-factors
			}

        // One dimensional Discrete Fourier Transform
        public static void DFT( Complex[] data, FourierDirection direction )
        {
            int			n = data.Length;
            double		arg, cos, sin;
            Complex[]	dst = new Complex[n];

            // for each destination element
            for ( int i = 0; i < n; i++ )
            {
                dst[i] = Complex.Zero;

                arg = - (int) direction * 2.0 * System.Math.PI * (double) i / (double) n;

                // sum source elements
                for ( int j = 0; j < n; j++ )
                {
                    cos = System.Math.Cos( j * arg );
                    sin = System.Math.Sin( j * arg );

                    dst[i].Re += (float) ( data[j].Re * cos - data[j].Im * sin );
                    dst[i].Im += (float) ( data[j].Re * sin + data[j].Im * cos );
                }
            }

            // copy elements
            if ( direction == FourierDirection.Forward )
            {
                // devide also for forward transform
                for ( int i = 0; i < n; i++ )
                {
                    data[i].Re = dst[i].Re / n;
                    data[i].Im = dst[i].Im / n;
                }
            }
            else
            {
                for ( int i = 0; i < n; i++ )
                {
                    data[i].Re = dst[i].Re;
                    data[i].Im = dst[i].Im;
                }
            }
        }

    /// <summary>
    /// Convolves or deconvolves a splitted complex-valued data set data[] (including any
    /// user supplied zero padding) with a response function response[].
    /// The result is returned in the splitted complex arrays resultre[] and resultim[]. All arrays including
    /// the scratch[] arrays must have the same dimensions (or larger).
    /// The data set (and of course the other arrays) can be either one-dimensional,
    /// two-dimensional, or three-dimensional, d = 1,2,3.  Each dimension must be 
    /// of the form n = (2**p) * (3**q) * (5**r), because of the underlying FFT.
    /// </summary>
    /// <param name="datare">
    /// The splitted complex-valued data set. Note, that you have to
    /// care for end effects by zero padding. This means, 
    /// that you have to pad the data with a number of zeros
    /// on one end equal to the maximal positive duration
    /// or maximal negative duration of the response function,
    /// whichever is larger!!</param>
    /// <param name="dataim">The imaginary part of the data array.</param>
    /// <param name="responsere">
    ///  The response function must be stored in wrap-around
    ///  order. This means that the first half of the array
    ///  response[] (in each dimension) contains the impulse
    ///  response function at positive times, while the second
    ///  half of the array contains the impulse response
    ///  function at negative times, counting down from the
    ///  element with the highest index. The array must have 
    ///  at least the size of the data array.
    /// </param>
    /// <param name="responseim">The imaginary part of the response array.</param>
    /// <param name="resultre">
    /// The real part of the result array. It must have 
    /// at least the size of the data array.
    /// </param>
    /// <param name="resultim">The imaginary part of the result array.</param>
    /// <param name="scratchre">
    ///  A work array. If a NULL pointer is passed the
    /// work array is allocated and freed auotomatically.
    /// If the array is given by the user it must have 
    /// at least the size of the data array.
    /// </param>
    /// <param name="scratchim">
    ///  A work array. If a NULL pointer is passed the
    /// work array is allocated and freed auotomatically.
    /// If the array is given by the user it must have 
    /// at least the size of the data array.
    /// </param>
    /// <param name="isign">
    /// If isign == forward a convolution is performed. 
    /// If isign == inverse then a deconvolution is performed.
    /// </param>
    /// <returns>
    /// In the case of a convolution (isign == forward) the value "true" is returned
    /// always. In the case of deconvolution (isign == inverse) the value "false" is
    /// returned if the FFT transform of the response function is exactly zero for 
    /// some value. This indicates that the original convolution has lost all 
    /// information at this particular frequency, so that a reconstruction is not
    /// possible. If the transform of the response function is non-zero everywhere
    /// the deconvolution can be performed and the value "true" is returned.
    /// </returns>
 
    public bool Convolute (
      double[] datare, double[] dataim, 
      double[] responsere, double[] responseim,
      double[] resultre, double[] resultim,
      double[] scratchre, double[] scratchim,
      FourierDirection isign)
    {
      // return status
      bool status = true;

      // get total size of data array
      int size = 0;
      if (ndim == 0) 
      {
        throw new ArithmeticException("Convolute: no dimensions have been specified");
      } 
      else if (ndim == 1) 
      {
        size = dim[0];
      } 
      else if (ndim == 2) 
      {
        size = row_order ? (dim[0] * id) : (id * dim[1]);
      } 
      else if (ndim == 3) 
      {
        size = row_order ? (dim[0] * dim[1] * id) : (id * dim[1] * dim[2]);
      }

      // allocate the scratch array
      //bool auto_scratch = false;
      if ( null==scratchre )
        scratchre = new  double[size];
      if ( null==scratchim )
        scratchim = new  double[size];

      //---------------------------------------------------------------------------//
      //  1-dimensional convolution (original data not are overwritten)
      //---------------------------------------------------------------------------//

      if (ndim == 1) 
//.........这里部分代码省略.........

        /// <summary>
        /// Voert een 2 dimensionale Fourier analyse uit op een afbeelding
        /// </summary>
        /// <param name="invoer">De complexe waarden van de afbeelding</param>
        /// <param name="dir">FourierDirection.Forward or FourierDirection.Backward</param>
        /// <returns>Fourier resultset</returns>
        private static ComplexGetal[,] FFT2D(ComplexGetal[,] invoer, FourierDirection dir)
        {
            float[] reëel = new float[invoer.GetLength(1)];
            float[] imaginair = new float[invoer.GetLength(1)];
            ComplexGetal[,] resultaat = invoer;

            //rijen
            for (int i = 0; i < invoer.GetLength(0); i++) {
                // Plaats de waarden in een 1 dimensionale array
                for (int j = 0; j < invoer.GetLength(1); j++) {
                    reëel[j] = invoer[i, j].Reëel;
                    imaginair[j] = invoer[i, j].Imaginair;
                }

                // Voer een 1 dimensionale fourier analyse uit op de huidige rij
                FFT1D(dir, reëel, imaginair);

                // Plaats de waarden in het resultaat object
                for (int j = 0; j < resultaat.GetLength(1); j++) {
                    resultaat[i, j].Reëel = reëel[j];
                    resultaat[i, j].Imaginair = imaginair[j];
                }
            }

            //kolommen - gebruik nu het resultaat object, anders gaan de berekeningen van de rijen verloren
            for (int i = 0; i < resultaat.GetLength(1); i++) {
                // Plaats de waarden in een 1 dimensionale array
                for (int j = 0; j < resultaat.GetLength(0); j++) {
                    reëel[j] = resultaat[j, i].Reëel;
                    imaginair[j] = resultaat[j, i].Imaginair;
                }

                // Voer een 1 dimensionale fourier analyse uit op de huidige kolom
                FFT1D(dir, reëel, imaginair);

                // Plaats de waarden in het resultaat object
                for (int j = 0; j < resultaat.GetLength(0); j++) {
                    resultaat[j, i].Reëel = reëel[j];
                    resultaat[j, i].Imaginair = imaginair[j];
                }
            }

            return resultaat;
        }

        /// <summary>
        /// This function computes an in-place complex-to-complex FFT
        ///    x and y are the real and imaginary arrays of 2^m points.
        ///    Both arrays must have equal size and must be a power of 2. 
        ///    dir = 1 gives forward transform
        ///    dir = -1 gives reverse transform
        ///    Formula: forward
        ///             N-1
        ///              ---
        ///            1 \         - j k 2 pi n / N
        ///    X(K) = --- > x(n) e                  = Forward transform
        ///            N /                            n=0..N-1
        ///              ---
        ///             n=0
        ///    Formula: reverse
        ///             N-1
        ///             ---
        ///             \          j k 2 pi n / N
        ///    X(n) =    > x(k) e                  = Inverse transform
        ///             /                             k=0..N-1
        ///             ---
        ///             k=0
        /// The result is found in the arrays x and y. 
        /// </summary>
        /// <param name="dir">FourierDirection.Forward or FourierDirection.Backward</param>
        /// <param name="x">real values</param>
        /// <param name="y">imaginary values</param>
        private static void FFT1D(FourierDirection dir, float[] x, float[] y)
        {
            long nn, i, i1, j, k, i2, l, l1, l2;
            float c1, c2, tx, ty, t1, t2, u1, u2, z;

            if (x.Length != y.Length)
                throw new FormatException("Real values and imaginary values arrays lengths do not match");

            int m = (int)Math.Log(x.Length, 2);
            if (m != Math.Log(x.Length, 2))
                throw new FormatException("Data arrays lenght is no power of two");

            /* Calculate the number of points */
            nn = 1;
            for (i = 0; i < m; i++)
                nn *= 2;
            /* Do the bit reversal */
            i2 = nn >> 1;
            j = 0;
            for (i = 0; i < nn - 1; i++)
            {
                if (i < j)
                {
                    tx = x[i];
                    ty = y[i];
                    x[i] = x[j];
                    y[i] = y[j];
                    x[j] = tx;
                    y[j] = ty;
                }
                k = i2;
                while (k <= j)
                {
                    j -= k;
                    k >>= 1;
                }
                j += k;
            }
            /* Compute the FFT */
            c1 = -1f;
            c2 = 0f;
            l2 = 1;
            for (l = 0; l < m; l++)
            {
                l1 = l2;
                l2 <<= 1;
                u1 = 1;
                u2 = 0;
                for (j = 0; j < l1; j++)
                {
                    for (i = j; i < nn; i += l2)
                    {
                        i1 = i + l1;
                        t1 = u1 * x[i1] - u2 * y[i1];
                        t2 = u1 * y[i1] + u2 * x[i1];
                        x[i1] = x[i] - t1;
                        y[i1] = y[i] - t2;
                        x[i] += t1;
                        y[i] += t2;
                    }
                    z = u1 * c1 - u2 * c2;
                    u2 = u1 * c2 + u2 * c1;
                    u1 = z;
                }
                c2 = (float)Math.Sqrt((1f - c1) / 2f);
                if (dir == FourierDirection.Forward)
                    c2 = -c2;
                c1 = (float)Math.Sqrt((1f + c1) / 2f);
            }
            /* Scaling for forward transform */
            if (dir == FourierDirection.Forward)
            {
                for (i = 0; i < nn; i++)
//.........这里部分代码省略.........

		/// <summary>
		/// Compute a 2D fast fourier transform on a data set of complex numbers
		/// </summary>
		/// <param name="data"></param>
		/// <param name="xLength"></param>
		/// <param name="yLength"></param>
		/// <param name="direction"></param>
		public static void	FFT2( ComplexF[] data, int xLength, int yLength, FourierDirection direction ) {

		int xInc = 1;
			int yInc = xLength;

			if( xLength > 1 ) {
				SyncLookupTableLength( xLength );
				for( int y = 0; y < yLength; y ++ ) {
					int xStart = y * yInc;
					LinearFFT_Quick( data, xStart, xInc, xLength, direction );
				}
			}
		
			if( yLength > 1 ) {
				SyncLookupTableLength( yLength );
				for( int x = 0; x < xLength; x ++ ) {
					int yStart = x * xInc;
					LinearFFT_Quick( data, yStart, yInc, yLength, direction );
				}
			}
		}

 /// <summary>
 /// Compute a 1D fast Fourier transform of a dataset of complex numbers.
 /// </summary>
 /// <param name="data"></param>
 /// <param name="direction"></param>
 public static void FFT( ComplexF[] data, FourierDirection direction )
 {
     if( data == null ) {
         throw new ArgumentNullException( "data" );
     }
     FourierBase.FFT(data, data.Length, direction);
 }

        /// <summary>
        /// Compute a 1D fast Fourier transform of a dataset of complex numbers.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void FFT_Quick( Complex[] data, int length, FourierDirection direction )
        {
            /*if( data == null ) {
                throw new ArgumentNullException( "data" );
            }
            if( data.Length < length ) {
                throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );
            }
            if( Fourier.IsPowerOf2( length ) == false ) {
                throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );
            }

            Fourier.SyncLookupTableLength( length );   */

            int ln = Fourier.Log2( length );

            // reorder array
            Fourier.ReorderArray( data );

            // successive doubling
            int N = 1;
            int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;

            for( int level = 1; level <= ln; level ++ ) {
                int M = N;
                N <<= 1;

                double[] uRLookup = _uRLookup[ level, signIndex ];
                double[] uILookup = _uILookup[ level, signIndex ];

                for( int j = 0; j < M; j ++ ) {
                    double uR = uRLookup[j];
                    double uI = uILookup[j];

                    for( int even = j; even < length; even += N ) {
                        int odd	 = even + M;

                        double	r = data[ odd ].Re;
                        double	i = data[ odd ].Im;

                        double	odduR = r * uR - i * uI;
                        double	odduI = r * uI + i * uR;

                        r = data[ even ].Re;
                        i = data[ even ].Im;

                        data[ even ].Re	= r + odduR;
                        data[ even ].Im	= i + odduI;

                        data[ odd ].Re	= r - odduR;
                        data[ odd ].Im	= i - odduI;
                    }
                }
            }
        }

        // Get rotation of complex number
        private static Complex[] GetComplexRotation( int numberOfBits, FourierDirection direction )
        {
            int		directionIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;

            // check if the array is already calculated
            if ( complexRotation[numberOfBits - 1, directionIndex] == null )
            {
                int			n = 1 << (numberOfBits - 1);
                float		uR = 1.0f;
                float		uI = 0.0f;
                double		angle = System.Math.PI / n * (int) direction;
                float		wR = (float) System.Math.Cos( angle );
                float		wI = (float) System.Math.Sin( angle );
                float		t;
                Complex[]	rotation = new Complex[n];

                for ( int i = 0; i < n; i++ )
                {
                    rotation[i] = new Complex(uR, uI);
                    t = uR * wI + uI * wR;
                    uR = uR * wR - uI * wI;
                    uI = t;
                }

                complexRotation[numberOfBits - 1, directionIndex] = rotation;
            }
            return complexRotation[numberOfBits - 1, directionIndex];
        }

        /// <summary>
        /// Compute a 1D real-symmetric fast fourier transform.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void RFFT( float[] data, int length, FourierDirection direction )
        {
            if( data == null ) {
                throw new ArgumentNullException( "data" );
            }
            if( data.Length < length ) {
                throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );
            }
            if( Fourier.IsPowerOf2( length ) == false ) {
                throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );
            }

            float	c1 = 0.5f, c2;
            float	theta	= (float) Math.PI / (length/2);

            if( direction == FourierDirection.Forward ) {
                c2 = -0.5f;
                FFT( data, length/2, direction );
            }
            else {
                c2 = 0.5f;
                theta = - theta;
            }

            float	wtemp = (float) Math.Sin( 0.5*theta );
            float	wpr = -2 * wtemp*wtemp;
            float	wpi	=(float)  Math.Sin( theta );
            float	wr	= 1 + wpr;
            float	wi	= wpi;

            // do / undo packing
            for( int i = 1; i < length/4; i ++ ) {
                int a = 2*i;
                int b = length - 2*i;
                float	h1r	= c1 * ( data[ a ] + data[ b ] );
                float	h1i	= c1 * ( data[ a+1 ] - data[ b+1 ] );
                float	h2r	= -c2 * ( data[ a+1 ] + data[ b+1 ] );
                float	h2i	= c2 * ( data[ a ] - data[ b ] );
                data[ a ]	= h1r + wr*h2r - wi*h2i;
                data[ a+1 ]	= h1i + wr*h2i + wi*h2r;
                data[ b ]	= h1r - wr*h2r + wi*h2i;
                data[ b+1 ]	= -h1i + wr*h2i + wi*h2r;
                wr = (wtemp = wr) * wpr - wi * wpi + wr;
                wi = wi * wpr + wtemp * wpi + wi;
            }

            if( direction == FourierDirection.Forward ) {
                float  hir = data[0];
                data[0] = hir + data[1];
                data[1] = hir - data[1];
            }
            else {
                float  hir = data[0];
                data[0] = c1 * ( hir + data[1] );
                data[1] = c1 * ( hir - data[1] );
                Fourier.FFT( data, length/2, direction );
            }
        }