public class FFT
{
public static Boolean IsPowerOfTwo(UInt32 x)
{
return ((x != 0) && (x & (x - 1)) == 0);
}
public static UInt32 NextPowerOfTwo(UInt32 x)
{
x = x - 1;
x = x | (x >> 1);
x = x | (x >> 2);
x = x | (x >> 4);
x = x | (x >> 8);
x = x | (x >> 16);
return x + 1;
}
public static UInt32 NumberOfBitsNeeded(UInt32 PowerOfTwo)
{
if (PowerOfTwo > 0)
{
for (UInt32 i = 0, mask = 1; ; i++, mask <<= 1)
{
if ((PowerOfTwo & mask) != 0)
return i;
}
}
return 0;
}
public static UInt32 ReverseBits(UInt32 index, UInt32 NumBits)
{
UInt32 i, rev;
for (i = rev = 0; i < NumBits; i++)
{
rev = (rev << 1) | (index & 1);
index >>= 1;
}
return rev;
}
public struct Complex
{
public double Re, Im;
public Complex(double x, double y)
{
this.Re = x; this.Im = y;
}
public Complex Mul(Complex C2)
{
Complex aux = new Complex();
aux.Re = this.Re * C2.Re - this.Im * C2.Im;
aux.Im = this.Re * C2.Im + this.Im * C2.Re;
return aux;
}
public Complex Add(Complex C2)
{
Complex aux = new Complex();
aux.Re = this.Re + C2.Re;
aux.Im = this.Im + C2.Im;
return aux;
}
public Complex Sub(Complex C2)
{
Complex aux = new Complex();
aux.Re = this.Re - C2.Re;
aux.Im = this.Im - C2.Im;
return aux;
}
public double Magnitude()
{
double Len = Math.Sqrt(this.Re * this.Re + this.Im * this.Im);
double bb = Math.Abs(Math.Tan(this.Re / Len)) / Math.Tan(this.Re / Len);
return (Len * bb );
}
public Complex Rev()
{
Complex aux = new Complex();
aux.Re = this.Re;
aux.Im = -1.0d * this.Im;
return aux;
}
}
public List<Complex> Cos_Fun(int N)
{
List<Complex> cc = new List<Complex>(N);
Complex aux;
double c = 2.0d * DDC_PI / (double)N;
for (int i = 0; i < N; i++)
{
aux = new Complex(
Math.Cos((double)(3 * i) * c)
, Math.Sin((double)(6 * i * i) * c)
);
cc.Add(aux);
}
return cc;
}
public const Double DDC_PI = 3.14159265358979323846;
UInt32 Block_Size = 2;
public Complex Get_W(int block_size, int index)
{
Complex c = new Complex(1.0d, 1.0d);
c.Re = c.Re * Math.Cos((-2 * DDC_PI / block_size) * index);
c.Im = c.Im * Math.Sin((-2 * DDC_PI / block_size) * index);
return c;
}
public List<Complex> Init_FFT(List<Complex> x)
{
Complex auxx = new Complex(0.0d, 0.0d);
List<Complex> aux = new List<Complex>(x.Count);
UInt32 count;
if (IsPowerOfTwo((UInt32)x.Count))
{
count = (UInt32)x.Count;
}
else
{
count = NextPowerOfTwo((UInt32)x.Count);
}
UInt32 b = NumberOfBitsNeeded(count);
aux = new List<Complex>((int)count);
Complex cccc= new Complex();
for (int i = 0; i < (int)count; i++)
{
aux.Add(cccc);
}
for (int i = x.Count; i < count; i++)
{
x.Add(auxx);
}
for (int i = 0; i < count; i++)
{
aux[i] = (x[(int)ReverseBits((UInt32)i, b)]);
}
x.Clear();
return aux;
}
public List<Complex> FFT_Calc(List<Complex> Xx)
{
Block_Size = 2; int index = 0; int indexx = 0;
int Block_num = Xx.Count / (int)Block_Size;
Complex auxx = new Complex(0.0d, 0.0d);
List<Complex> aux = new List<Complex>(Xx.Count);
for (int i = 0; i < Xx.Count; i++)
aux.Add(auxx);
int test = 0;
while (Block_Size <= Xx.Count)
{
index = 0;
Block_num = Xx.Count / (int)Block_Size;
for (int i = 0; i < Block_num; i++)
{
indexx = 0; index = i * (int)Block_Size;
for (int j = 0; j < (int)Block_Size / 2; j++)
{
auxx = Get_W((int)Block_Size, indexx);
aux[index] = Xx[index];
auxx = auxx.Mul(Xx[index + (int)Block_Size / 2]);
aux[index] = aux[index].Add(auxx);
test++;
auxx = Get_W((int)Block_Size, indexx + (int)Block_Size / 2);
aux[index + (int)Block_Size / 2] = Xx[index];
auxx = auxx.Mul(Xx[index + (int)Block_Size / 2]);
aux[index + (int)Block_Size / 2] = aux[index + (int)Block_Size / 2].Add(auxx);
test++;
indexx++;
index++;
}
}
for (int k = 0; k < aux.Count; k++)
Xx[k] = aux[k];
Block_Size <<= 1;
}
return aux;
}
public List<Complex> Inv_FFT_Calc(List<Complex> Xx)
{
Block_Size = 2; int index = 0; int indexx = 0;
int Block_num = Xx.Count / (int)Block_Size;
Complex auxx = new Complex(0.0d, 0.0d);
List<Complex> aux = new List<Complex>(Xx.Count);
for (int i = 0; i < Xx.Count; i++)
aux.Add(auxx);
int test = 0;
while (Block_Size <= Xx.Count)
{
index = 0;
Block_num = Xx.Count / (int)Block_Size;
for (int i = 0; i < Block_num; i++)
{
indexx = 0; index = i * (int)Block_Size;
for (int j = 0; j < (int)Block_Size / 2; j++)
{
auxx = Get_W((int)Block_Size, -indexx);
aux[index] = Xx[index];
auxx = auxx.Mul(Xx[index + (int)Block_Size / 2]);
aux[index] = aux[index].Add(auxx);
test++;
auxx = Get_W((int)Block_Size, -indexx + (int)Block_Size / 2);
aux[index + (int)Block_Size / 2] = Xx[index];
auxx = auxx.Mul(Xx[index + (int)Block_Size / 2]);
aux[index + (int)Block_Size / 2] = aux[index + (int)Block_Size / 2].Add(auxx);
test++;
indexx++;
index++;
}
}
for (int k = 0; k < aux.Count; k++)
Xx[k] = aux[k];
Block_Size <<= 1;
}
Complex c = new Complex(1.0d / aux.Count, 0.0d);
List<Complex> auxxx = new List<Complex>(aux.Count);
for (int k = 0; k < aux.Count; k++)
{
c.Re = aux[k].Re / (double)aux.Count; c.Im = 0.0d;
auxxx.Add(c);
}
return auxxx;
}
public System.Drawing.Bitmap Draw_Spectral(float[] Percent, int From, int To, Color c, System.Drawing.Size Bitmap_Size)
{
Bitmap bb = new Bitmap(Bitmap_Size.Width, Bitmap_Size.Height);
if (Percent != null)
{
System.Drawing.Graphics g = Graphics.FromImage(bb);
g.SmoothingMode = System.Drawing.Drawing2D.SmoothingMode.HighQuality;
g.Clear(Color.White);
float max = 0.0f;
for (int i = From; i < To; i++)
{
if (max < Percent[i])
max = Percent[i];
}
float xStep = ((float)bb.Width - 10.0f) / (float)(To - From);
System.Drawing.SolidBrush sb1 = new SolidBrush(Color.LightBlue);
System.Drawing.Pen p1 = new Pen(c, 1.5f);
float yStep = 1.0f * ((float)bb.Height - 10.0f) / max;
float x = 5.0f, y = bb.Height - 5.0f;
for (int i = From; i < To; i++)
{
g.DrawLine(p1
, x
, bb.Height - 5.0f
, x
, y - (Percent[i] * yStep));
x += xStep;
}
}
return bb;
}
public void Draw_Wave(List<Complex> Spectral, System.Windows.Forms.Panel p, Color c, float D_Width)
{
if (Spectral != null)
{
System.Drawing.Graphics g = p.CreateGraphics();
g.SmoothingMode = System.Drawing.Drawing2D.SmoothingMode.HighSpeed;
g.Clear(Color.White);
float max = 0.0f;
for (int i = 0; i < Spectral.Count; i++)
{
if (max < Spectral[i].Magnitude())
max = (float)Spectral[i].Magnitude();
}
float xStep = ((float)p.Width - 10.0f) / Spectral.Count;
System.Drawing.Pen p1 = new Pen(c, D_Width);
float yStep = 1.0f * ((float)p.Height - 10.0f) / max;
float x = 5.0f, y = p.Height - 5.0f;
for (int i = 0; i < Spectral.Count; i++)
{
g.DrawLine(p1
, x
, p.Height - 5.0f
, x
, y - ((float)Spectral[i].Magnitude() * yStep));
x += xStep;
}
}
}
public enum Complex_Part
{
Real,
Image,
Magenda
}
public void Draw_Spectral(List<Complex> Spectral, System.Windows.Forms.Panel p,Complex_Part cp, Color c, float D_Width)
{
if (Spectral != null)
{
System.Drawing.Graphics g = p.CreateGraphics();
g.SmoothingMode = System.Drawing.Drawing2D.SmoothingMode.HighSpeed;
g.Clear(Color.White);
float max = 0.0f; float min = 0.0f; ;
int freq_step = (Spectral.Count) * 10 / 100;
int f_s = freq_step;
if (cp == Complex_Part.Image)
{
for (int i = 0; i < (Spectral.Count); i++)
{
if (max < Spectral[i].Im)
max = (float)Spectral[i].Im;
if (min > Spectral[i].Im)
min = (float)Spectral[i].Im;
}
}
else if (cp == Complex_Part.Real)
{
for (int i = 0; i < (Spectral.Count); i++)
{
if (max < Spectral[i].Re)
max = (float)Spectral[i].Re;
if (min > Spectral[i].Re)
min = (float)Spectral[i].Re;
}
}
else if (cp == Complex_Part.Magenda)
{
for (int i = 0; i < (Spectral.Count); i++)
{
if (max < Spectral[i].Magnitude())
max = (float)Spectral[i].Magnitude();
if (min > Spectral[i].Magnitude())
min = (float)Spectral[i].Magnitude();
}
}
float xStep = ((float)p.Width - 10.0f) / (Spectral.Count);
System.Drawing.Pen p1 = new Pen(c, D_Width);
System.Drawing.Pen p2 = new Pen(Color.Black, 1.0f);
float yStep = 1.0f * ((float)p.Height - 20.0f) / (max - min);
float x = 5.0f, y = p.Height - (10.0f - yStep * min);
Font f = new Font("Tahoma", 6.5f, FontStyle.Regular);
SolidBrush b = new SolidBrush(Color.Black);
for (int i = 0; i < (Spectral.Count); i++)
{
if (i == f_s)
{
g.DrawString(f_s + "", f, b, x - 25.0f, 5.0f);
f_s += freq_step;
}
g.DrawLine(p2
, x - xStep
, p.Height - (10.0f - yStep * min)
, x + xStep
, p.Height - (10.0f - yStep * min));
if (cp == Complex_Part.Image)
{
g.DrawLine(p1
, x
, p.Height - (10.0f - yStep * min)
, x
, y - (((float)Spectral[i].Im) * yStep));
}
else if (cp == Complex_Part.Real)
{
g.DrawLine(p1
, x
, p.Height - (10.0f - yStep * min)
, x
, y - (((float)Spectral[i].Re) * yStep));
}
else if (cp == Complex_Part.Magenda)
{
g.DrawLine(p1
, x
, p.Height - (10.0f - yStep * min)
, x
, y - (((float)Spectral[i].Magnitude()) * yStep));
}
x += xStep;
}
}
}
public void Draw_Spectral(List<Complex> Spectral, System.Windows.Forms.Panel p, Color c, float D_Width)
{
if (Spectral != null)
{
System.Drawing.Graphics g = p.CreateGraphics();
g.SmoothingMode = System.Drawing.Drawing2D.SmoothingMode.HighSpeed;
g.Clear(Color.White);
float max = 0.0f; float min = 0.0f; ;
int freq_step = (Spectral.Count / 2) * 10 / 100;
int f_s = freq_step;
for (int i = 0; i < (Spectral.Count/2); i++)
{
if (max < Spectral[i].Magnitude())
max = (float)Spectral[i].Magnitude();
if (min > Spectral[i].Magnitude())
min = (float)Spectral[i].Magnitude();
}
float xStep = ((float)p.Width - 10.0f) / (Spectral.Count / 2);
System.Drawing.Pen p1 = new Pen(c, D_Width);
System.Drawing.Pen p2 = new Pen(Color.Black, 1.0f);
float yStep = 1.0f * ((float)p.Height - 20.0f) / (max - min);
float x = 5.0f, y = p.Height - (10.0f - yStep * min);
Font f = new Font("Tahoma", 6.5f, FontStyle.Regular);
SolidBrush b = new SolidBrush(Color.Black);
for (int i = 0; i < (Spectral.Count / 2); i++)
{
if (i == f_s)
{
g.DrawString(f_s + "", f, b, x - 25.0f, 5.0f);
f_s += freq_step;
}
g.DrawLine(p1
, x
, p.Height - (10.0f - yStep * min)
, x
, y - (((float)Spectral[i].Magnitude()) * yStep ));
g.DrawLine(p2
, x - xStep
, p.Height - (10.0f - yStep * min)
, x + xStep
, p.Height - (10.0f - yStep * min));
x += xStep;
}
}
}
public void Draw_Spectral(double[] Spectral, System.Windows.Forms.Panel p, Color c, float D_Width)
{
if (Spectral != null)
{
System.Drawing.Graphics g = p.CreateGraphics();
g.SmoothingMode = System.Drawing.Drawing2D.SmoothingMode.HighSpeed;
g.Clear(Color.White);
float max = 0.0f;
int freq_step = (Spectral.Length) * 5 / 100;
int f_s = freq_step;
for (int i = 0; i < Spectral.Length; i++)
{
if (max < Math.Abs(Spectral[i]))
max = (float)Math.Abs(Spectral[i]);
}
float xStep = ((float)p.Width - 10.0f) / (Spectral.Length);
System.Drawing.Pen p1 = new Pen(c, D_Width);
float yStep = 1.0f * ((float)p.Height - 10.0f) / max;
float x = 5.0f, y = p.Height - 5.0f;
Font f = new Font("Tahoma", 6.5f, FontStyle.Regular);
SolidBrush b = new SolidBrush(Color.Black);
for (int i = 0; i < (Spectral.Length); i++)
{
if (i == f_s)
{
g.DrawString(f_s + "", f, b, x - 25.0f, 5.0f);
f_s += freq_step;
}
g.DrawLine(p1
, x
, p.Height - 5.0f
, x
, y - ((float)Math.Abs(Spectral[i]) * yStep));
x += xStep;
}
}
}
public List<Complex> InitRSine(int n)
{
List<Complex> Data_In = new List<Complex>();
Complex aux = new Complex();
for (int i = 0; i < n; i++)
{
aux = new Complex();
aux.Re =
Math.Cos((double)(45.0) * i * (0.5 * Math.PI / n))
+ 4.0d * Math.Cos((double)(60.0) * i * (0.1 * Math.PI / n))
;
aux.Im = 0.0d;
Data_In.Add(aux);
}
return Data_In;
}
}