算法笔记
C 语言实现离散傅利叶变换（DFT）代码如下：
/*********************************************************************************************************
* 函数名称： DFT
* 函数功能： 离散傅利叶变换（DFT）
* 输入参数： x：双精度实型一维数组，长度为 n。存放要变换数据的实部。
* y：双精度实型一维数组，长度为 n。存放要变换数据的虚部。
* a：双精度实型一维数组，长度为 n。存放变换结果的实部。
* b：双精度实型一维数组，长度为 n。存放变换结果的虚部。
* n：整形变量。数据长度
* sign：整形变量。当 sign=1 时，子函数 DFT() 计算离散傅利叶正变换；当 sign=-1 时，做反变换
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年08月14日
* 注 意：
*********************************************************************************************************/
void DFT(double* x, double* y, double* a, double* b, int n, int sign)
{
int i, k;
double c, d, q, w, s;
q = 6.28318530718 / (double)n;
for (k = 0; k < n; k++)
{
w = k * q;
a[k] = b[k] = 0.0;
for (i = 0; i < n; i++)
{
d = i * w;
c = cos(d);
s = sin(d) * sign;
a[k] += c * x[i] + s * y[i];
b[k] += c * y[i] - s * x[i];
}
}
if (sign == -1)
{
c = 1.0 / (double)n;
for (k = 0; k < n; k++)
{
a[k] = c * a[k];
b[k] = c * b[k];
}
}
}C 语言实现快速傅利叶变换（FFT），代码如下：
/*********************************************************************************************************
* 函数名称： FFT
* 函数功能： 快速傅利叶变换
* 输入参数： x：双精度实型一维数组，长度为 n。开始时存放要变换数据的实部，最后存放变换结果的实部。
* y：双精度实型一维数组，长度为 n。开始时存放要变换数据的虚部，最后存放变换结果的虚部。
* n：整形变量。数据长度。必须是 2 的整数次幂，即 n=2^m 。
* sign：整形变量。当 sign=1 时，子函数 FFT() 计算离散傅利叶正变换；当 sign=-1 时，做反变换
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年08月14日
* 注 意：
*********************************************************************************************************/
void FFT(double* x, double* y, int n, int sign)
{
int i, j, k, l, m, n1, n2;
double c, c1, e, s, s1, t, tr, ti;
for (j = 1, i = 1; i < 16; i++)
{
m = i;
j = 2 * j;
if (j == n)
{
break;
}
}
n1 = n - 1;
for (j = 0, i = 0; i < n1; i++)
{
if (i < j)
{
tr = x[j];
ti = y[j];
x[j] = x[i];
y[j] = y[i];
x[i] = tr;
y[i] = ti;
}
k = n / 2;
while (k < (j + 1))
{
j = j - k;
k = k / 2;
}
j = j + k;
}
n1 = 1;
for (l = 1; l <= m; l++)
{
n1 = 2 * n1;
n2 = n1 / 2;
e = 3.14159265359 / n2;
c = 1.0;
s = 0.0;
c1 = cos(e);
s1 = -sign * sin(e);
for (j = 0; j < n2; j++)
{
for (i = j; i < n; i += n1)
{
k = i + n2;
tr = c * x[k] - s * y[k];
ti = c * y[k] + s * x[k];
x[k] = x[i] - tr;
y[k] = y[i] - ti;
x[i] = x[i] + tr;
y[i] = y[i] + ti;
}
t = c;
c = c * c1 - s * s1;
s = t * s1 + s * c1;
}
}
if (sign == -1)
{
for (i = 0; i < n; i++)
{
x[i] /= n;
y[i] /= n;
}
}
}
分裂基快速傅利叶变换代码如下。
/*********************************************************************************************************
* 函数名称： SRFFT
* 函数功能： 分裂基快速傅利叶变换
* 输入参数： x：双精度实型一维数组，长度为 n。开始时存放要变换数据的实部，最后存放变换结果的实部。
* y：双精度实型一维数组，长度为 n。开始时存放要变换数据的虚部，最后存放变换结果的虚部。
* n：整形变量。数据长度。必须是 2 的整数次幂，即 n=2^m 。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年08月18日
* 注 意：
*********************************************************************************************************/
void SRFFT(double* x, double* y, int n)
{
int i, j, k, m, i1, i2, i3, n1, n2, n4, id, is;
double a, e, a3, r1, r2;
double s1, s2, s3, cc1, cc3, ss1, ss3;
for (j = 1, i = 1; i < 16; i++)
{
m = i;
j = 2 * j;
if (j == n)
{
break;
}
}
n2 = 2 * n;
for (k = 1; k < m; k++)
{
n2 = n2 / 2;
n4 = n2 / 4;
e = 6.28318530718 / n2;
a = 0;
for (j = 0; j < n4; j++)
{
a3 = 3 * a;
cc1 = cos(a);
ss1 = sin(a);
cc3 = cos(a3);
ss3 = sin(a3);
a = ((double)j + 1.0) * e;
is = j;
id = 2 * n2;
do
{
for (i = is; i < (n - 1); i = i + id)
{
i1 = i + n4;
i2 = i1 + n4;
i3 = i2 + n4;
r1 = x[i] - x[i2];
x[i] = x[i] + x[i2];
r2 = x[i1] - x[i3];
x[i1] = x[i1] + x[i3];
s1 = y[i] - y[i2];
y[i] = y[i] + y[i2];
s2 = y[i1] - y[i3];
y[i1] = y[i1] + y[i3];
s3 = r1 - s2;
r1 = r1 + s2;
s2 = r2 - s1;
r2 = r2 + s1;
x[i2] = r1 * cc1 - s2 * ss1;
y[i2] = -s2 * cc1 - r1 * ss1;
x[i3] = s3 * cc3 + r2 * ss3;
y[i3] = r2 * cc3 - s3 * ss3;
}
is = 2 * id - n2 + j;
id = 4 * id;
}while(is < (n - 1));
}
}
is = 0;
id = 4;
do
{
for (i = is; i < n; i = i + id)
{
i1 = i + 1;
r1 = x[i];
r2 = y[i];
x[i] = r1 + x[i1];
y[i] = r2 + y[i1];
x[i1] = r1 - x[i1];
y[i1] = r2 - y[i1];
}
is = 2 * id - 2;
id = 4 * id;
}while(is < (n - 1));
n1 = n - 1;
for (j = 0, i = 0; i < n1; i++)
{
if (i < j)
{
r1 = x[j];
s1 = y[j];
x[j] = x[i];
y[j] = y[i];
x[i] = r1;
y[i] = s1;
}
k = n / 2;
while (k < (j + 1))
{
j = j - k;
k = k / 2;
}
j = j + k;
}
}实序列快速傅利叶变换（一），代码如下：
/*********************************************************************************************************
* 函数名称： RFFT
* 函数功能： 实序列快速傅利叶变换（一）
* 输入参数： x：双精度实型一维数组，长度为 n。开始时存放要变换的实数据，最后存放变换结果的前 (n/2)+1 个值。
* 存储顺序为：[Re(0),Re(1),...,Re(n/2),Im((n/2)-1),...,Im(1)]。
* 其中 Re(0)=X(0)，Re(n/2)=X(n/2)。根据 X(k) 的共轭对称性，很容易写出后半部分的值
* n：整形变量。数据长度。必须是 2 的整数次幂，即 n=2^m 。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年08月18日
* 注 意：
*********************************************************************************************************/
void RFFT(double* x, double n)
{
int i, j, k, m, i1, i2, i3, i4, n1, n2, n4;
double a, e, cc, ss, xt, t1, t2;
for (j = 1, i = 1; i < 16; i++)
{
m = i;
j = 2 * j;
if (j == n)
{
break;
}
}
n1 = n - 1;
for (j = 0, i = 0; i < n1; i++)
{
if (i < j)
{
xt = x[j];
x[j] = x[i];
x[i] = xt;
}
k = n / 2;
while (k < (j + 1))
{
j = j - k;
k = k / 2;
}
j = j + k;
}
for (i = 0; i < n; i += 2)
{
xt = x[i];
x[i] = xt + x[i + 1];
x[i + 1] = xt - x[i + 1];
}
n2 = 1;
for (k = 2; k <= m; k++)
{
n4 = n2;
n2 = 2 * n4;
n1 = 2 * n2;
e = 6.28318530718 / n1;
for (i = 0; i < n; i += n1)
{
xt = x[i];
x[i] = xt + x[i + n2];
x[i + n2] = xt - x[i + n2];
x[i + n2 + n4] = -x[i + n2 + n4];
a = e;
for (j = 1; j <= (n4 - 1); j++)
{
i1 = i + j;
i2 = i - j + n2;
i3 = i + j + n2;
i4 = i - j + n1;
cc = cos(a);
ss = sin(a);
a = a + e;
t1 = cc * x[i3] + ss * x[i4];
t2 = ss * x[i3] - cc * x[i4];
x[i4] = x[i2] - t2;
x[i3] = -x[i2] - t2;
x[i2] = x[i1] - t1;
x[i1] = x[i1] + t1;
}
}
}
}实序列快速傅利叶变换（二），代码如下：
/*********************************************************************************************************
* 函数名称： RLFFT
* 函数功能： 实序列快速傅利叶变换（二）
* 输入参数： x：双精度实型一维数组，长度为 n。开始时存放要变换的实数据，最后存放变换结果的前 (n/2)+1 个值。
* 存储顺序为：[Re(0),Re(1),...,Re(n/2),Im((n/2)-1),...,Im(1)]。
* 其中 Re(0)=X(0)，Re(n/2)=X(n/2)。根据 X(k) 的共轭对称性，很容易写出后半部分的值
* n：整形变量。数据长度。必须是 2 的整数次幂，即 n=2^m 。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年08月18日
* 注 意：
*********************************************************************************************************/
void RLFFT(double* x, int n)
{
int i, n1;
double a, c, e, s, fr, fi, gr, gi, *f, *g;
f = malloc((n / 2) * sizeof(double));
g = malloc((n / 2) * sizeof(double));
n1 = n/ 2;
e = 3.141592653589793 / n1;
for (i = 0; i < n1; i++)
{
f[i] = x[2 * i];
g[i] = x[2 * i + 1];
}
FFT(f, g, n1, 1);
x[0] = f[0] + g[0];
x[n1] = f[0] - g[0];
for (i = 1; i < n1; i++)
{
fr = (f[i] + f[n1 - i]) / 2;
fi = (g[i] - g[n1 - i]) / 2;
gr = (g[n1 - i] + g[i]) / 2;
gi = (f[n1 - i] - f[i]) / 2;
a = i * e;
c = cos(a);
s = sin(a);
x[i] = fr + c * gr + s * gi;
x[n - i] = fi + c * gi - s * gr;
}
free(f);
free(g);
}Chirp Z 变换，代码如下：
#include "math.h"
#include "stdlib.h"
#include "malloc.h"
/*********************************************************************************************************
* 函数名称： CZT
* 函数功能： Chirp Z 变换
* 输入参数： xr：双精度实型一维数组，其长度大于或等于（n + m - 1），且是 2 的 整数次幂。
* 开始时存放输入数据的实部，最后存放变换结果的实部。
* xi：双精度实型一维数组，其长度大于或等于（n + m - 1），且是 2 的 整数次幂。
* 开始时存放输入数据的虚部，最后存放变换结果的虚部。
* n：整形变量。输入数据的长度。
* m：整形变量。输出数据的长度，即频率采样点数。
* f1：双精度实型变量。起始数字频率，单位为 Hz-s。
* f2：双精度实型变量。终止数字频率。单位为 Hz-s。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年08月29日
* 注 意：
*********************************************************************************************************/
void CZT(double* xr, double* xi, int n, int m, double f1, double f2)
{
int i, j, n1, n2, len;
double e, t, ar, ai, ph, pi, tr, ti, *wr, *wr1, *wi, *wi1;
len = n + m - 1;
for (j = 1, i = 1; i < 16; i++)
{
j = 2 * j;
if (j >= len)
{
len = j;
break;
}
}
wr = malloc(len * sizeof(double));
wi = malloc(len * sizeof(double));
wr1 = malloc(len * sizeof(double));
wi1 = malloc(len * sizeof(double));
pi = 3.14159265358979;
ph = 2.0 * pi * (f2 - f1)/(m - 1);
n1 = (n >= m)? n:m;
for (i = 0; i < n1; i++)
{
e = ph * i * i / 2.0;
wr[i] = cos(e);
wi[i] = sin(e);
wr1[i] = wr[i];
wi1[i] = -wi[i];
}
n2 = len - n + 1;
for (i = m; i < n2; i++)
{
wr[i] = 0.0;
wi[i] = 0.0;
}
for (i = n2; i < len; i++)
{
j = len - i;
wr[i] = wr[j];
wi[i] = wi[j];
}
FFT(wr, wi, len, 1);
ph = -2.0 * pi * f1;
for (i = 0; i < n; i++)
{
e = ph * i;
ar = cos(e);
ai = sin(e);
tr = ar * wr1[i] - ai * wi1[i];
ti = ai * wr1[i] + ar * wi1[i];
t = xr[i] * tr - xi[i] * ti;
xi[i] = xr[i] * ti + xi[i] * tr;
xr[i] = t;
}
for (i = n; i < len; i++)
{
xr[i] = 0.0;
xi[i] = 0.0;
}
FFT(xr, xi, len, 1);
for (i = 0; i < len; i++)
{
tr = xr[i] * wr[i] - xi[i] * wi[i];
xi[i] = xr[i] * wi[i] + xi[i] * wr[i];
xr[i] = tr;
}
FFT(xr, xi, len, -1);
for (i = 0; i < m; i++)
{
tr = xr[i] * wr1[i] - xi[i] * wi1[i];
xi[i] = xr[i] * wi1[i] + xi[i] * wr1[i];
xr[i] = tr;
}
free(wr);
free(wi);
free(wr1);
free(wi1);
}