一、理论基础
任何周期函数都可以看成是不同振幅，不同相位的正弦波的叠加。
1. 三角函数系的正交性：
三角函数系：
三角函数系中任意两个不同的函数之积在[-π,π]上的积分为0。
2. 将周期函数展开为三角级数
假设有周期函数 f(x)=f(x+2π)，周期为2*PI；则函数可以展开为下面的三角级数：
书中展示的公式为： 
以上两个公式本质上是形同的,下面给给出求解上式参数的过程：
2 求取参数an,bn
3. 离散傅里叶变换：
已知，连续的傅里叶变换公式为：
4. 快速傅里叶变化（FFT）：
离散傅里叶变换的计算量是巨大的。对一幅 MxN 的图像，他需要M^2 N^2 次的复数乘法和加法。因此，一般使用离散傅里叶变换的快速算法——快速傅里叶变换。
4.1 算法思想：
复数域中将一般形式转化为指数形式的公式推导：
2. n次单位根：
n次单位根在复数域相当于是一个：
幅角为： 2π/n
长度为： 1
的虚数。
3. n次单位根的性质：
4.2.多项式系数表示法：
4.3.多项式点值表示法：
因为n次函数可由n+1个点确定下来（可以将每一个点列一个n次方程）。
由上面的式子可以唯一的确定一组多项式系数，a0,a1,a2,…an,因此多项式可以使用下面的点值表示法：
点值表示法的乘法运算：
两个多项式P,Q分别取点(x，y1)和(x，y2) ，P*Q就会去点（x，y1*y2）；
DFT(离散傅里叶变换):
求多项式的点值表示法：
想要求出一个多项式的点值表示法，需要选出n+1个数分别带入到多项式里面。带入一个数的复杂度是Θ(n) Θ(n)的，那么总复杂度就是Θ(n2) 。因为单位根有下面的性质，所以可以尝试去n次单位根来作为n+1个数。
FFT算法：
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace DFT与FFT算法
{
/// <summary>
/// 实现复数运算的类
/// </summary>
class complex
{
//复数为a+bi
//其中real为实部字段
private double real;
/// <summary>
/// 复数的实部（属性）
/// </summary>
public double Real
{
get { return real; }
set { real = value; }
}
//image为虚部字段。
private double image;
/// <summary>
/// 复数的虚部（属性）
/// </summary>
public double Image
{
get { return image; }
set { image = value; }
}
/// <summary>
/// 默认构造函数
/// </summary>
public complex ()
:this(0,0)
{
}
/// <summary>
/// 只有实部的构造函数
/// </summary>
/// <param name="real"></param>
public complex(double real)
: this(real, 0)
{
}
/// <summary>
/// 由实部和虚部创建的构造函数
/// </summary>
/// <param name="real"></param>
/// <param name="image"></param>
public complex(double real,double image)
{
this.real = real;
this.image = image;
}
/// <summary>
/// 重载加法符号
/// </summary>
/// <param name="c1"></param>
/// <param name="c2"></param>
/// <returns></returns>
public static complex operator +(complex c1,complex c2)
{
return new complex(c1.real + c2.real, c1.image + c2.image);
}
/// <summary>
/// 重载减法
/// </summary>
/// <param name="c1"></param>
/// <param name="c2"></param>
/// <returns></returns>
public static complex operator -(complex c1, complex c2)
{
return new complex(c1.real - c2.real, c1.image - c2.image);
}
/// <summary>
/// 重载乘法
/// </summary>
/// <param name="c1"></param>
/// <param name="c2"></param>
/// <returns></returns>
public static complex operator *(complex c1, complex c2)
{
return new complex(c1.real * c2.real - c1.image * c2.image, c1.image * c2.real + c1.real * c2.image);
}
/// <summary>
/// 复数的模
/// </summary>
/// <returns></returns>
public double ToModul()
{
return Math.Sqrt(real * real + image * image);
}
/// <summary>
/// 重载ToString 方法
/// </summary>
/// <returns></returns>
public override string ToString()
{
if(Real==0 && Image ==0)
{
return string.Format("{0}", 0);
}
if(Real==0 && (Image !=1 && Image !=-1))
{
return string.Format("{0} i", Image);
}
if(Image==0)
{
return string.Format("{0}", Real);
}
if(Image==1)
{
return string.Format("i");
}
if (Image == -1)
{
return string.Format("- i");
}
if (Image < 0)
{
return string.Format("{0} - {1} i", Real, -Image);
}
return string.Format("{0} + {1} i", Real, Image);
}
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace DFT与FFT算法
{
/// <summary>
/// 实现运算
/// </summary>
class DFT
{
/// <summary>
/// 求复数数组的模数组
/// </summary>
/// <param name="input"></param>
/// <returns></returns>
public static double[] Cmp2Mdl(complex[] input)
{
//确定模数组的长度
double[] output = new double[input.Length];
//对复数求模
for(int i=0;i<input.Length;i++)
{
if(input[i].Real>0)
{
output[i] = -input[i].ToModul();
}
else
{
output[i] = input[i].ToModul();
}
}
return output;
}
/// <summary>
/// 傅里叶变换以及傅里叶反变换
/// </summary>
/// <param name="input"></param>
/// <param name="invert，false FFT变换，true IFFT变换"></param>
/// <returns></returns>
public static complex[] FFT(double[] input,bool invert)
{
//由输入序列确定输出序列的长度
complex[] output = new complex[input.Length];
///将输入的实数数组转换为复数数组
for(int i=0;i<input.Length;i++)
{
output[i]=new complex(input[i]);
}
//返回FFT或者IFFT序列
return output = FFT(output, invert);
}
/// <summary>
/// 傅里叶变换以及傅里叶反变换,递归蝶形运算
/// </summary>
/// <param name="input"></param>
/// <param name="invert，false FFT变换，true IFFT变换"></param>
/// <returns></returns>
public static complex[] FFT(complex[] input,bool invert)
{
///如果输入序列只有一个元素，输出这个元素并返回
if(input.Length==1)
{
return new complex[] { input[0]};
}
//输入序列的长度
int length = input.Length;
//输入序列长度的一半
int half = length / 2;
//由输入序列的长度确定输出序列的长度
complex[] output = new complex[length];
//正变换旋转因子的基数
double fac = -2.0 * Math.PI / length;
//反变换旋转因子的基数是正变换的相反数
if(invert)
{
fac = -fac;
}
//序列中下标中为偶数的点
complex[] evens = new complex[half];
for(int i=0;i<half;i++)
{
evens[i] = input[2 * i];
}
//求偶数点FFT或IFFT的结果，递归实现多级蝶形运算
complex[] evenResult = FFT(evens,invert);
//序列中下标为奇数的点
complex[] odds=new complex[half];
for(int i=0;i<half;i++)
{
odds[i] = input[2 * i + 1];
}
//求奇数点FFT或IFFT的结果，递归实现多级蝶形运算
complex[] oddResult = FFT(odds, invert);
for(int k=0;k<half;k++)
{
///旋转因子
double fack = fac * k;
//进行蝶形运算
complex oddpart = oddResult[k] * new complex(Math.Cos(fack), Math.Sin(fack));
output[k] = evenResult[k] + oddpart;
output[k + half] = evenResult[k] - oddpart;
}
//返回FFT或IFFT的结果
return output;
}
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace DFT与FFT算法
{
class Program
{
static void Main(string[] args)
{
double[] input = {
2,3,4,5,6,7,8,9
//0.0928,0.0353,0.6124,0.6085,0.0158,0.0164,0.1901,0.5869
};
complex[] output_complex = new complex[input.Length];
output_complex = DFT.FFT(input, false);//正变换
Console.WriteLine();
Console.WriteLine("i\treal\timage");
Console.WriteLine();
for (int i = 0; i < input.Length; i++)
{
Console.WriteLine("{0}\t{1}\t{2}", i, output_complex[i].Real, output_complex[i].Image);
}
Console.WriteLine("------------------------------------------");
output_complex = DFT.FFT(output_complex, true);//反变换
Console.WriteLine();
Console.WriteLine("i\treal\timage");
Console.WriteLine();
for (int i = 0; i < input.Length; i++)
{
Console.WriteLine("{0}\t{1}\t{2}", i, output_complex[i].Real, output_complex[i].Image);
}
input = DFT.Cmp2Mdl(output_complex);
Console.WriteLine("------------------------------------------");
//反变换完了除以数据长度
for (int i = 0; i < input.Length; i++)
{
input[i] = -input[i] / input.Length;
}
foreach (double a in input)
{
Console.WriteLine("{0}", a);
}
Console.ReadKey();
}
}
}
5657
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
