fourier.cs

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

			// copy to buffer
			float[]	buffer = null;
			LockBufferF( length * 2, ref buffer );
			int j = start;
			for( int i = 0; i < length * 2; i ++ ) {
				buffer[ i ] = data[ j ];
				j += inc;
			}

			FFT( buffer, length, direction );

			// copy from buffer
			j = start;
			for( int i = 0; i < length; i ++ ) {
				data[ j ] = buffer[ i ];
				j += inc;
			}
			UnlockBufferF( ref buffer );
		}

		private static void	LinearFFT_Quick( float[] data, int start, int inc, int length, FourierDirection direction ) {
			/*Debug.Assert( data != null );
			Debug.Assert( start >= 0 );
			Debug.Assert( inc >= 1 );
			Debug.Assert( length >= 1 );
			Debug.Assert( ( start + inc * ( length - 1 ) ) * 2 < data.Length );*/
			
			// copy to buffer
			float[]	buffer = null;
			LockBufferF( length * 2, ref buffer );
			int j = start;
			for( int i = 0; i < length * 2; i ++ ) {
				buffer[ i ] = data[ j ];
				j += inc;
			}

			FFT_Quick( buffer, length, direction );

			// copy from buffer
			j = start;
			for( int i = 0; i < length; i ++ ) {
				data[ j ] = buffer[ i ];
				j += inc;
			}
			UnlockBufferF( ref buffer );
		}

		//======================================================================================
		//======================================================================================
		
		static private bool			_bufferCFLocked	= false;
		static private ComplexF[]	_bufferCF		= new ComplexF[ 0 ];

		static private void		LockBufferCF( int length, ref ComplexF[] buffer ) {
			Debug.Assert( length >= 0 );
			Debug.Assert( _bufferCFLocked == false );
			
			_bufferCFLocked = true;
			if( length != _bufferCF.Length ) {
				_bufferCF	= new ComplexF[ length ];
			}
			buffer =	_bufferCF;
		}
		static private void		UnlockBufferCF( ref ComplexF[] buffer ) {
			Debug.Assert( _bufferCF == buffer );
			Debug.Assert( _bufferCFLocked == true );
			
			_bufferCFLocked = false;
			buffer = null;
		}

		private static void	LinearFFT( ComplexF[] data, int start, int inc, int length, FourierDirection direction ) {
			Debug.Assert( data != null );
			Debug.Assert( start >= 0 );
			Debug.Assert( inc >= 1 );
			Debug.Assert( length >= 1 );
			Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );
			
			// copy to buffer
			ComplexF[]	buffer = null;
			LockBufferCF( length, ref buffer );
			int j = start;
			for( int i = 0; i < length; i ++ ) {
				buffer[ i ] = data[ j ];
				j += inc;
			}

			FFT( buffer, length, direction );

			// copy from buffer
			j = start;
			for( int i = 0; i < length; i ++ ) {
				data[ j ] = buffer[ i ];
				j += inc;
			}
			UnlockBufferCF( ref buffer );
		}

		private static void	LinearFFT_Quick( ComplexF[] data, int start, int inc, int length, FourierDirection direction ) {
			/*Debug.Assert( data != null );
			Debug.Assert( start >= 0 );
			Debug.Assert( inc >= 1 );
			Debug.Assert( length >= 1 );
			Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );	*/
			
			// copy to buffer
			ComplexF[]	buffer = null;
			LockBufferCF( length, ref buffer );
			int j = start;
			for( int i = 0; i < length; i ++ ) {
				buffer[ i ] = data[ j ];
				j += inc;
			}

			FFT( buffer, length, direction );

			// copy from buffer
			j = start;
			for( int i = 0; i < length; i ++ ) {
				data[ j ] = buffer[ i ];
				j += inc;
			}
			UnlockBufferCF( ref buffer );
		}

		//======================================================================================
		//======================================================================================
		
		static private bool			_bufferCLocked	= false;
		static private Complex[]	_bufferC		= new Complex[ 0 ];

		static private void		LockBufferC( int length, ref Complex[] buffer ) {
			Debug.Assert( length >= 0 );
			Debug.Assert( _bufferCLocked == false );
			
			_bufferCLocked = true;
			if( length >= _bufferC.Length ) {
				_bufferC	= new Complex[ length ];
			}
			buffer =	_bufferC;
		}
		static private void		UnlockBufferC( ref Complex[] buffer ) {
			Debug.Assert( _bufferC == buffer );
			Debug.Assert( _bufferCLocked == true );
			
			_bufferCLocked = false;
			buffer = null;
		}

		private static void	LinearFFT( Complex[] data, int start, int inc, int length, FourierDirection direction ) {
			Debug.Assert( data != null );
			Debug.Assert( start >= 0 );
			Debug.Assert( inc >= 1 );
			Debug.Assert( length >= 1 );
			Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );
			
			// copy to buffer
			Complex[]	buffer = null;
			LockBufferC( length, ref buffer );
			int j = start;
			for( int i = 0; i < length; i ++ ) {
				buffer[ i ] = data[ j ];
				j += inc;
			}

			FFT( buffer, length, direction );

			// copy from buffer
			j = start;
			for( int i = 0; i < length; i ++ ) {
				data[ j ] = buffer[ i ];
				j += inc;
			}
			UnlockBufferC( ref buffer );
		}

		private static void	LinearFFT_Quick( Complex[] data, int start, int inc, int length, FourierDirection direction ) {
			/*Debug.Assert( data != null );
			Debug.Assert( start >= 0 );
			Debug.Assert( inc >= 1 );
			Debug.Assert( length >= 1 );
			Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );*/
			
			// copy to buffer
			Complex[]	buffer = null;
			LockBufferC( length, ref buffer );
			int j = start;
			for( int i = 0; i < length; i ++ ) {
				buffer[ i ] = data[ j ];
				j += inc;
			}

			FFT_Quick( buffer, length, direction );

			// copy from buffer
			j = start;
			for( int i = 0; i < length; i ++ ) {
				data[ j ] = buffer[ i ];
				j += inc;
			}
			UnlockBufferC( ref buffer );
		}

		//======================================================================================
		//======================================================================================

		/// <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( Fourier.IsPowerOf2( length ) == true );

			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;

				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>
		/// 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_Quick( float[] data, int length, FourierDirection direction ) {
			/*Debug.Assert( data != null );
			Debug.Assert( data.Length >= length*2 );
			Debug.Assert( Fourier.IsPowerOf2( length ) == true );

			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;

				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>
		/// 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( ComplexF[] 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;

				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 even = j; even < length; even += N ) {
						int odd	 = even + M;
						
						float	r = data[ odd ].Re;
						float	i = data[ odd ].Im;

						float	odduR = r * uR - i * uI;
						float	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( ComplexF[] 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;

				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 even = j; even < length; even += N ) {
						int odd	 = even + M;
						
						float	r = data[ odd ].Re;
						float	i = data[ odd ].Im;

						float	odduR = r * uR - i * uI;
						float	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="direction"></param>
		public static void	FFT( ComplexF[] data, FourierDirection direction ) {
			if( data == null ) {
				throw new ArgumentNullException( "data" );
			}