fourier.cs

C#写的复数和傅立叶变换算法

			Fourier.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( 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;
					}
				}
			}

		}
		
		/// <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;
					}
				}
			}

		}

		/// <summary>
		/// Compute a 1D real-symmetric fast fourier transform.
		/// </summary>
		/// <param name="data"></param>
		/// <param name="direction"></param>
		public static void	RFFT( float[] data, FourierDirection direction ) {
			if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			Fourier.RFFT( data, data.Length, direction );
		}
		
		/// <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 );
			}
		}
		
		/// <summary>
		/// Compute a 2D fast fourier transform on a data set of complex numbers (represented as pairs of floats)
		/// </summary>
		/// <param name="data"></param>
		/// <param name="xLength"></param>
		/// <param name="yLength"></param>
		/// <param name="direction"></param>
		public static void	FFT2( float[] data, int xLength, int yLength, FourierDirection direction ) {
			if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			if( data.Length < xLength*yLength*2 ) {
				throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength * 2' parameter" );
			}
			if( Fourier.IsPowerOf2( xLength ) == false ) {
				throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );
			}
			if( Fourier.IsPowerOf2( yLength ) == false ) {
				throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );
			}

			int xInc = 1;
			int yInc = xLength;

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

		/// <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 ) {
			if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			if( data.Length < xLength*yLength ) {
				throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength' parameter" );
			}
			if( Fourier.IsPowerOf2( xLength ) == false ) {
				throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );
			}
			if( Fourier.IsPowerOf2( yLength ) == false ) {
				throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );
			}

			int xInc = 1;
			int yInc = xLength;

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

		/// <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( Complex[] data, int xLength, int yLength, FourierDirection direction ) {
			if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			if( data.Length < xLength*yLength ) {
				throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength' parameter" );
			}
			if( Fourier.IsPowerOf2( xLength ) == false ) {
				throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );
			}
			if( Fourier.IsPowerOf2( yLength ) == false ) {
				throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );
			}

			int xInc = 1;
			int yInc = xLength;

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

		/// <summary>
		/// Compute a 3D 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="zLength"></param>
		/// <param name="direction"></param>
		public static void	FFT3( ComplexF[] data, int xLength, int yLength, int zLength, FourierDirection direction ) {
			if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			if( data.Length < xLength*yLength*zLength ) {
				throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength * zLength' parameter" );
			}
			if( Fourier.IsPowerOf2( xLength ) == false ) {
				throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );
			}
			if( Fourier.IsPowerOf2( yLength ) == false ) {
				throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );
			}
			if( Fourier.IsPowerOf2( zLength ) == false ) {
				throw new ArgumentOutOfRangeException( "zLength", zLength, "must be a power of 2" );
			}

			int xInc = 1;
			int yInc = xLength;
			int zInc = xLength * yLength;

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

		/// <summary>
		/// Compute a 3D 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="zLength"></param>
		/// <param name="direction"></param>
		public static void	FFT3( Complex[] data, int xLength, int yLength, int zLength, FourierDirection direction ) {
			if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			if( data.Length < xLength*yLength*zLength ) {
				throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength * zLength' parameter" );
			}
			if( Fourier.IsPowerOf2( xLength ) == false ) {
				throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );
			}
			if( Fourier.IsPowerOf2( yLength ) == false ) {
				throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );
			}
			if( Fourier.IsPowerOf2( zLength ) == false ) {
				throw new ArgumentOutOfRangeException( "zLength", zLength, "must be a power of 2" );
			}

			int xInc = 1;
			int yInc = xLength;
			int zInc = xLength * yLength;

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