快速傅立叶变换(FFT)的C#代码
2023-09-14 08:56:49 时间
这个代码是从《快速傅立叶变换(FFT)的C++实现与Matlab实验》这篇文章里的源代码转换而来,请注意查看原文。
在这里自己转换成了C#代码,并作了一些改动,主要是对N值的确定,原文的N值为常量1024,自己通过对输入的数组的长度来确定N值,N值的确定符合2的n次方,函数返回N值。通过作者提供的测试变量进行测试,能得到相应的结果。代码如下:
FFT代码:
using System;
using System.Collections.Generic;
using System.Text; namespace ConsoleApplication1
{
/// summary
/// 快速傅立叶变换(Fast Fourier Transform)。
/// /summary
public class TWFFT
{
private TWFFT()
{
} private static void bitrp(float[] xreal, float[] ximag, int n)
{
// 位反转置换 Bit-reversal Permutation
int i, j, a, b, p; for (i = 1, p = 0; i i *= 2)
{
p++;
}
for (i = 0; i i++)
{
a = i;
b = 0;
for (j = 0; j j++)
{
b = b * 2 + a % 2;
a = a / 2;
}
if (b i)
{
float t = xreal[i];
xreal[i] = xreal[b];
xreal[b] = t; t = ximag[i];
ximag[i] = ximag[b];
ximag[b] = t;
}
}
} public static int FFT(float[] xreal, float[] ximag)
{
//n值为2的N次方
int n = 2;
while (n = xreal.Length)
{
n *= 2;
}
n /= 2; // 快速傅立叶变换,将复数 x 变换后仍保存在 x 中,xreal, ximag 分别是 x 的实部和虚部
float[] wreal = new float[n / 2];
float[] wimag = new float[n / 2];
float treal, timag, ureal, uimag, arg;
int m, k, j, t, index1, index2; bitrp(xreal, ximag, n); // 计算 1 的前 n / 2 个 n 次方根的共轭复数 Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
arg = (float)(-2 * Math.PI / n);
treal = (float)Math.Cos(arg);
timag = (float)Math.Sin(arg);
wreal[0] = 1.0f;
wimag[0] = 0.0f;
for (j = 1; j n / 2; j++)
{
wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
} for (m = 2; m m *= 2)
{
for (k = 0; k k += m)
{
for (j = 0; j m / 2; j++)
{
index1 = k + j;
index2 = index1 + m / 2;
t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中的下标为 t
treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
ureal = xreal[index1];
uimag = ximag[index1];
xreal[index1] = ureal + treal;
ximag[index1] = uimag + timag;
xreal[index2] = ureal - treal;
ximag[index2] = uimag - timag;
}
}
} return n;
} public static int IFFT(float[] xreal, float[] ximag)
{
//n值为2的N次方
int n = 2;
while (n = xreal.Length)
{
n *= 2;
}
n /= 2; // 快速傅立叶逆变换
float[] wreal = new float[n / 2];
float[] wimag = new float[n / 2];
float treal, timag, ureal, uimag, arg;
int m, k, j, t, index1, index2; bitrp(xreal, ximag, n); // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
arg = (float)(2 * Math.PI / n);
treal = (float)(Math.Cos(arg));
timag = (float)(Math.Sin(arg));
wreal[0] = 1.0f;
wimag[0] = 0.0f;
for (j = 1; j n / 2; j++)
{
wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
} for (m = 2; m m *= 2)
{
for (k = 0; k k += m)
{
for (j = 0; j m / 2; j++)
{
index1 = k + j;
index2 = index1 + m / 2;
t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中的下标为 t
treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
ureal = xreal[index1];
uimag = ximag[index1];
xreal[index1] = ureal + treal;
ximag[index1] = uimag + timag;
xreal[index2] = ureal - treal;
ximag[index2] = uimag - timag;
}
}
} for (j = 0; j j++)
{
xreal[j] /= n;
ximag[j] /= n;
} return n;
}
}
}
这里只给出了函数FFT()的,函数IFFT()可以进行相应测试。测试代码: using System;
using System.Collections.Generic;
using System.Text; namespace ConsoleApplication1
{
class Program
{
static void Main(string[] args)
{
float[] a = {
0.5751f,0.4514f,0.0439f,0.0272f,0.3127f,0.0129f,0.3840f,0.6831f,
0.0928f,0.0353f,0.6124f ,0.6085f,0.0158f,0.0164f,0.1901f,0.5869f}; float[] b = {
0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,
0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}; int n = TWFFT.FFT(a,b); Console.WriteLine("FFT:"); Console.WriteLine("FFT的返回值N = {0}",n);
Console.WriteLine(); Console.WriteLine("i\ta\tb"); Console.WriteLine();
for (int i = 0; i i++)
{
Console.WriteLine("{0}\t{1}\t{2}",i,a[i],b[i]);
} Console.ReadLine();
}
}
【MATLAB】离散余弦变换滤波算法(DCT) 之前介绍的所有滤波算法都是空间域滤波算法(即2D滤波算法)。离散余弦变换滤波算法(DCT)属于频率域滤波算法(即3D滤波算法)。
using System.Collections.Generic;
using System.Text; namespace ConsoleApplication1
{
/// summary
/// 快速傅立叶变换(Fast Fourier Transform)。
/// /summary
public class TWFFT
{
private TWFFT()
{
} private static void bitrp(float[] xreal, float[] ximag, int n)
{
// 位反转置换 Bit-reversal Permutation
int i, j, a, b, p; for (i = 1, p = 0; i i *= 2)
{
p++;
}
for (i = 0; i i++)
{
a = i;
b = 0;
for (j = 0; j j++)
{
b = b * 2 + a % 2;
a = a / 2;
}
if (b i)
{
float t = xreal[i];
xreal[i] = xreal[b];
xreal[b] = t; t = ximag[i];
ximag[i] = ximag[b];
ximag[b] = t;
}
}
} public static int FFT(float[] xreal, float[] ximag)
{
//n值为2的N次方
int n = 2;
while (n = xreal.Length)
{
n *= 2;
}
n /= 2; // 快速傅立叶变换,将复数 x 变换后仍保存在 x 中,xreal, ximag 分别是 x 的实部和虚部
float[] wreal = new float[n / 2];
float[] wimag = new float[n / 2];
float treal, timag, ureal, uimag, arg;
int m, k, j, t, index1, index2; bitrp(xreal, ximag, n); // 计算 1 的前 n / 2 个 n 次方根的共轭复数 Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
arg = (float)(-2 * Math.PI / n);
treal = (float)Math.Cos(arg);
timag = (float)Math.Sin(arg);
wreal[0] = 1.0f;
wimag[0] = 0.0f;
for (j = 1; j n / 2; j++)
{
wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
} for (m = 2; m m *= 2)
{
for (k = 0; k k += m)
{
for (j = 0; j m / 2; j++)
{
index1 = k + j;
index2 = index1 + m / 2;
t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中的下标为 t
treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
ureal = xreal[index1];
uimag = ximag[index1];
xreal[index1] = ureal + treal;
ximag[index1] = uimag + timag;
xreal[index2] = ureal - treal;
ximag[index2] = uimag - timag;
}
}
} return n;
} public static int IFFT(float[] xreal, float[] ximag)
{
//n值为2的N次方
int n = 2;
while (n = xreal.Length)
{
n *= 2;
}
n /= 2; // 快速傅立叶逆变换
float[] wreal = new float[n / 2];
float[] wimag = new float[n / 2];
float treal, timag, ureal, uimag, arg;
int m, k, j, t, index1, index2; bitrp(xreal, ximag, n); // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
arg = (float)(2 * Math.PI / n);
treal = (float)(Math.Cos(arg));
timag = (float)(Math.Sin(arg));
wreal[0] = 1.0f;
wimag[0] = 0.0f;
for (j = 1; j n / 2; j++)
{
wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
} for (m = 2; m m *= 2)
{
for (k = 0; k k += m)
{
for (j = 0; j m / 2; j++)
{
index1 = k + j;
index2 = index1 + m / 2;
t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中的下标为 t
treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
ureal = xreal[index1];
uimag = ximag[index1];
xreal[index1] = ureal + treal;
ximag[index1] = uimag + timag;
xreal[index2] = ureal - treal;
ximag[index2] = uimag - timag;
}
}
} for (j = 0; j j++)
{
xreal[j] /= n;
ximag[j] /= n;
} return n;
}
}
}
这里只给出了函数FFT()的,函数IFFT()可以进行相应测试。测试代码: using System;
using System.Collections.Generic;
using System.Text; namespace ConsoleApplication1
{
class Program
{
static void Main(string[] args)
{
float[] a = {
0.5751f,0.4514f,0.0439f,0.0272f,0.3127f,0.0129f,0.3840f,0.6831f,
0.0928f,0.0353f,0.6124f ,0.6085f,0.0158f,0.0164f,0.1901f,0.5869f}; float[] b = {
0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,
0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}; int n = TWFFT.FFT(a,b); Console.WriteLine("FFT:"); Console.WriteLine("FFT的返回值N = {0}",n);
Console.WriteLine(); Console.WriteLine("i\ta\tb"); Console.WriteLine();
for (int i = 0; i i++)
{
Console.WriteLine("{0}\t{1}\t{2}",i,a[i],b[i]);
} Console.ReadLine();
}
}
【MATLAB】离散余弦变换滤波算法(DCT) 之前介绍的所有滤波算法都是空间域滤波算法(即2D滤波算法)。离散余弦变换滤波算法(DCT)属于频率域滤波算法(即3D滤波算法)。
相关文章
- c# mysql executenonquery_C#与数据库访问技术总结(八)之ExecuteNonQuery方法
- C#交换两个变量值的几种方法总结分享
- 跨语言调用C#代码的新方式-DllExport
- C#实现图片分割方法与代码
- asp.net(c#)利用构造器链的代码
- C#系统热键注册实现代码
- C#jpg缩略图函数代码
- C#利用ODP.net连接Oracle数据库的操作方法
- C#实现简单打印的实例代码
- C#打印出正等腰三角形实例代码
- C#实现Web文件上传的两种方法实例代码
- 浅析C#与C++相关概念的比较
- c#list部分操作实现代码
- c#代码自动修改解决方案下任意文件实例
- C#计算两个时间差的方法代码分享
- C#实现word文件下载的代码
- C#图片按比例缩放的实现代码
- C#中Ilist与list的区别小结
- c#的params参数使用示例
- C#遍历DataSet控件实例总结
- C#利用Openxml读取Excel数据实例