📄 fourier.cs
字号:
/*
* BSD Licence:
* Copyright (c) 2001, 2002 Ben Houston [ ben@exocortex.org ]
* Exocortex Technologies [ www.exocortex.org ]
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. Neither the name of the <ORGANIZATION> nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*/
using System;
using System.Diagnostics;
using System.IO;
using System.Text;
//using Exocortex.Imaging;
namespace Exocortex.DSP {
// Comments? Questions? Bugs? Tell Ben Houston at ben@exocortex.org
// Version: May 4, 2002
/// <summary>
/// <p>Static functions for doing various Fourier Operations.</p>
/// </summary>
public class Fourier {
//======================================================================================
private Fourier() {
}
//======================================================================================
static private void Swap( ref float a, ref float b ) {
float temp = a;
a = b;
b = temp;
}
static private void Swap( ref double a, ref double b ) {
double temp = a;
a = b;
b = temp;
}
static private void Swap( ref ComplexF a, ref ComplexF b ) {
ComplexF temp = a;
a = b;
b = temp;
}
static private void Swap( ref Complex a, ref Complex b ) {
Complex temp = a;
a = b;
b = temp;
}
//-------------------------------------------------------------------------------------
private const int cMaxLength = 4096;
private const int cMinLength = 1;
private const int cMaxBits = 12;
private const int cMinBits = 0;
static private bool IsPowerOf2( int x ) {
return (x & (x - 1)) == 0;
//return ( x == Pow2( Log2( x ) ) );
}
static private int Pow2( int exponent ) {
if( exponent >= 0 && exponent < 31 ) {
return 1 << exponent;
}
return 0;
}
static private int Log2( int x ) {
if( x <= 65536 ) {
if( x <= 256 ) {
if( x <= 16 ) {
if( x <= 4 ) {
if( x <= 2 ) {
if( x <= 1 ) {
return 0;
}
return 1;
}
return 2;
}
if( x <= 8 )
return 3;
return 4;
}
if( x <= 64 ) {
if( x <= 32 )
return 5;
return 6;
}
if( x <= 128 )
return 7;
return 8;
}
if( x <= 4096 ) {
if( x <= 1024 ) {
if( x <= 512 )
return 9;
return 10;
}
if( x <= 2048 )
return 11;
return 12;
}
if( x <= 16384 ) {
if( x <= 8192 )
return 13;
return 14;
}
if( x <= 32768 )
return 15;
return 16;
}
if( x <= 16777216 ) {
if( x <= 1048576 ) {
if( x <= 262144 ) {
if( x <= 131072 )
return 17;
return 18;
}
if( x <= 524288 )
return 19;
return 20;
}
if( x <= 4194304 ) {
if( x <= 2097152 )
return 21;
return 22;
}
if( x <= 8388608 )
return 23;
return 24;
}
if( x <= 268435456 ) {
if( x <= 67108864 ) {
if( x <= 33554432 )
return 25;
return 26;
}
if( x <= 134217728 )
return 27;
return 28;
}
if( x <= 1073741824 ) {
if( x <= 536870912 )
return 29;
return 30;
}
// since int is unsigned it can never be higher than 2,147,483,647
// if( x <= 2147483648 )
// return 31;
// return 32;
return 31;
}
//-------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------
static private int ReverseBits( int index, int numberOfBits ) {
Debug.Assert( numberOfBits >= cMinBits );
Debug.Assert( numberOfBits <= cMaxBits );
int reversedIndex = 0;
for( int i = 0; i < numberOfBits; i ++ ) {
reversedIndex = ( reversedIndex << 1 ) | ( index & 1 );
index = ( index >> 1 );
}
return reversedIndex;
}
//-------------------------------------------------------------------------------------
static private int[][] _reversedBits = new int[ cMaxBits ][];
static private int[] GetReversedBits( int numberOfBits ) {
Debug.Assert( numberOfBits >= cMinBits );
Debug.Assert( numberOfBits <= cMaxBits );
if( _reversedBits[ numberOfBits - 1 ] == null ) {
int maxBits = Fourier.Pow2( numberOfBits );
int[] reversedBits = new int[ maxBits ];
for( int i = 0; i < maxBits; i ++ ) {
int oldBits = i;
int newBits = 0;
for( int j = 0; j < numberOfBits; j ++ ) {
newBits = ( newBits << 1 ) | ( oldBits & 1 );
oldBits = ( oldBits >> 1 );
}
reversedBits[ i ] = newBits;
}
_reversedBits[ numberOfBits - 1 ] = reversedBits;
}
return _reversedBits[ numberOfBits - 1 ];
}
//-------------------------------------------------------------------------------------
static private void ReorderArray( float[] data ) {
Debug.Assert( data != null );
int length = data.Length / 2;
Debug.Assert( Fourier.IsPowerOf2( length ) == true );
Debug.Assert( length >= cMinLength );
Debug.Assert( length <= cMaxLength );
int[] reversedBits = Fourier.GetReversedBits( Fourier.Log2( length ) );
for( int i = 0; i < length; i ++ ) {
int swap = reversedBits[ i ];
if( swap > i ) {
Fourier.Swap( ref data[ (i<<1) ], ref data[ (swap<<1) ] );
Fourier.Swap( ref data[ (i<<1) + 1 ], ref data[ (swap<<1) + 1 ] );
}
}
}
static private void ReorderArray( double[] data ) {
Debug.Assert( data != null );
int length = data.Length / 2;
Debug.Assert( Fourier.IsPowerOf2( length ) == true );
Debug.Assert( length >= cMinLength );
Debug.Assert( length <= cMaxLength );
int[] reversedBits = Fourier.GetReversedBits( Fourier.Log2( length ) );
for( int i = 0; i < length; i ++ ) {
int swap = reversedBits[ i ];
if( swap > i ) {
Fourier.Swap( ref data[ i<<1 ], ref data[ swap<<1 ] );
Fourier.Swap( ref data[ i<<1 + 1 ], ref data[ swap<<1 + 1 ] );
}
}
}
static private void ReorderArray( Complex[] data ) {
Debug.Assert( data != null );
int length = data.Length;
Debug.Assert( Fourier.IsPowerOf2( length ) == true );
Debug.Assert( length >= cMinLength );
Debug.Assert( length <= cMaxLength );
int[] reversedBits = Fourier.GetReversedBits( Fourier.Log2( length ) );
for( int i = 0; i < length; i ++ ) {
int swap = reversedBits[ i ];
if( swap > i ) {
Complex temp = data[ i ];
data[ i ] = data[ swap ];
data[ swap ] = temp;
}
}
}
static private void ReorderArray( ComplexF[] data ) {
Debug.Assert( data != null );
int length = data.Length;
Debug.Assert( Fourier.IsPowerOf2( length ) == true );
Debug.Assert( length >= cMinLength );
Debug.Assert( length <= cMaxLength );
int[] reversedBits = Fourier.GetReversedBits( Fourier.Log2( length ) );
for( int i = 0; i < length; i ++ ) {
int swap = reversedBits[ i ];
if( swap > i ) {
ComplexF temp = data[ i ];
data[ i ] = data[ swap ];
data[ swap ] = temp;
}
}
}
//======================================================================================
private static int[][] _reverseBits = null;
private static int _ReverseBits( int bits, int n ) {
int bitsReversed = 0;
for( int i = 0; i < n; i ++ ) {
bitsReversed = ( bitsReversed << 1 ) | ( bits & 1 );
bits = ( bits >> 1 );
}
return bitsReversed;
}
private static void InitializeReverseBits( int levels ) {
_reverseBits = new int[levels + 1][];
for( int j = 0; j < ( levels + 1 ); j ++ ) {
int count = (int) Math.Pow( 2, j );
_reverseBits[j] = new int[ count ];
for( int i = 0; i < count; i ++ ) {
_reverseBits[j][i] = _ReverseBits( i, j );
}
}
}
private static int _lookupTabletLength = -1;
private static double[,][] _uRLookup = null;
private static double[,][] _uILookup = null;
private static float[,][] _uRLookupF = null;
private static float[,][] _uILookupF = null;
private static void SyncLookupTableLength( int length ) {
Debug.Assert( length < 1024*10 );
Debug.Assert( length >= 0 );
if( length > _lookupTabletLength ) {
int level = (int) Math.Ceiling( Math.Log( length, 2 ) );
Fourier.InitializeReverseBits( level );
Fourier.InitializeComplexRotations( level );
//_cFFTData = new Complex[ Math2.CeilingBase( length, 2 ) ];
//_cFFTDataF = new ComplexF[ Math2.CeilingBase( length, 2 ) ];
_lookupTabletLength = length;
}
}
private static int GetLookupTableLength() {
return _lookupTabletLength;
}
private static void ClearLookupTables() {
_uRLookup = null;
_uILookup = null;
_uRLookupF = null;
_uILookupF = null;
_lookupTabletLength = -1;
}
private static void InitializeComplexRotations( int levels ) {
int ln = levels;
//_wRLookup = new float[ levels + 1, 2 ];
//_wILookup = new float[ levels + 1, 2 ];
_uRLookup = new double[ levels + 1, 2 ][];
_uILookup = new double[ levels + 1, 2 ][];
_uRLookupF = new float[ levels + 1, 2 ][];
_uILookupF = new float[ levels + 1, 2 ][];
int N = 1;
for( int level = 1; level <= ln; level ++ ) {
int M = N;
N <<= 1;
//float scale = (float)( 1 / Math.Sqrt( 1 << ln ) );
// positive sign ( i.e. [M,0] )
{
double uR = 1;
double uI = 0;
double angle = (double) Math.PI / M * 1;
double wR = (double) Math.Cos( angle );
double wI = (double) Math.Sin( angle );
_uRLookup[level,0] = new double[ M ];
_uILookup[level,0] = new double[ M ];
_uRLookupF[level,0] = new float[ M ];
_uILookupF[level,0] = new float[ M ];
for( int j = 0; j < M; j ++ ) {
_uRLookupF[level,0][j] = (float)( _uRLookup[level,0][j] = uR );
_uILookupF[level,0][j] = (float)( _uILookup[level,0][j] = uI );
double uwI = uR*wI + uI*wR;
uR = uR*wR - uI*wI;
uI = uwI;
}
}
{
// negative sign ( i.e. [M,1] )
double uR = 1;
double uI = 0;
double angle = (double) Math.PI / M * -1;
double wR = (double) Math.Cos( angle );
double wI = (double) Math.Sin( angle );
_uRLookup[level,1] = new double[ M ];
_uILookup[level,1] = new double[ M ];
_uRLookupF[level,1] = new float[ M ];
_uILookupF[level,1] = new float[ M ];
for( int j = 0; j < M; j ++ ) {
_uRLookupF[level,1][j] = (float)( _uRLookup[level,1][j] = uR );
_uILookupF[level,1][j] = (float)( _uILookup[level,1][j] = uI );
double uwI = uR*wI + uI*wR;
uR = uR*wR - uI*wI;
uI = uwI;
}
}
}
}
//======================================================================================
//======================================================================================
static private bool _bufferFLocked = false;
static private float[] _bufferF = new float[ 0 ];
static private void LockBufferF( int length, ref float[] buffer ) {
Debug.Assert( _bufferFLocked == false );
_bufferFLocked = true;
if( length >= _bufferF.Length ) {
_bufferF = new float[ length ];
}
buffer = _bufferF;
}
static private void UnlockBufferF( ref float[] buffer ) {
Debug.Assert( _bufferF == buffer );
Debug.Assert( _bufferFLocked == true );
_bufferFLocked = false;
buffer = null;
}
private static void LinearFFT( 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 );
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -