在上一篇内容
我们将傅里叶级数推导为傅里叶变换,而傅里叶变换计算的时候因为是一个积分,计算机并不是很好计算,所以要把积分换成一种累加形式,也就是本文要讨论的 离散傅里叶变化 DFT。
我们取上一篇的公式(7)
其中
因为傅里叶变化令 从而使一个累加的式子变成了一个积分,而DFT中 会根据输入的信号点数确定具体的值。具体计算公式为:
(注: 的计算方式是因为 的一个周期是 ,N为你采样的点数)
因此我们可以简化公式为
其中:
简单的换元法。将两个公式整合后,我们可以重新定义周期中两个点的距离,从而约掉 得到:
其中(1)为离散傅里叶变换,(2)为离散傅里叶逆变化。
代码实现:
先定义一个复数的结构体
public struct complex//定义复数
{ public double real; public double imag; // 写个函数用于显示 public void show() { if (Math.Abs(real) < 0.0001) real = 0; if (Math.Abs(imag) < 0.0001) imag = 0; if (imag > 0) Debug.Log(string.Format("{0} +{1}i", real, imag)); else if (imag == 0) Debug.Log(string.Format("{0}", real)); else Debug.Log(string.Format("{0} {1}i", real, imag)); } }
注:t为f数组的索引,n为F数组的索引,理清楚了代码很好理解
计算DFT,即已知一个 float 数组求频谱
public static complex[] calDFT(double[] f) //(信号,信号长度) { int N = f.Length; complex[] F = new complex[N]; for (int n = 0; n < N; n++) { // 声明 F[n].real = 0; F[n].imag = 0; for (int t = 0; t < N; t++) { // 计算 exp(-i * 2PI * n / N * t) complex temp; // 欧拉公式 exp(ix) = cos(x) + isin(x) temp.real = Math.Cos(-2 * Math.PI / N * t * n) * f[t]; temp.imag = Math.Sin(-2 * Math.PI / N * t * n) * f[t]; F[n].real = F[n].real + temp.real; F[n].imag = F[n].imag + temp.imag; } } return F; }
计算IDFT,即已知一个 float 数组求频谱
public static complex[] calIDFT(complex[] F) //(频谱) { int N = F.Length; complex[] f = new complex[N]; for (int t = 0; t < N; t++) { // 声明 f[t].real = 0; f[t].imag = 0; for (int n = 0; n < N; n++) { // 计算 exp(i * 2PI * n / N * t) complex temp; // 欧拉公式 exp(ix) = cos(x) + isin(x) double real = Math.Cos(2 * Math.PI * n / N * t); double imag = Math.Sin(2 * Math.PI * n / N * t); // 复数乘法 temp.real = F[n].real * real - F[n].imag * imag; temp.imag = F[n].real * imag + F[n].imag * real; f[t].real = f[t].real + temp.real; f[t].imag = f[t].imag + temp.imag; } f[t].real = f[t].real / N; f[t].imag = f[t].imag / N; } return f; }
随便输入一个float的数组做一下实验
double[] signal = new double[] { 14,12,10,8,6,10}; Debug.Log("----------计算离散傅里叶变化-------------"); var rate = calDFT(signal); foreach (var item in rate) { item.show(); } Debug.Log("----------计算反离散傅里叶变化------------"); var irate = calIDFT(rate); foreach (var item in irate) { item.show(); }
结果如下:
学过线性代数的会觉得有点像是 某个 维数很高的向量乘以 一个对应的矩阵,然后在乘以一个逆矩阵...转回来的过程。
我们记
得到DFT的矩阵:
以及IDFT的矩阵:
编辑于 2019-08-14