C语言实现快速傅里叶变换的算法

C语言实现快速傅里叶变换的算法

快速傅立叶/余弦/正弦变换
尺寸：一个
数据长度：2的幂
抽取：频率
基数：4，2

函数：

cdft：复数离散傅立叶变换
rdft：实离散傅里叶变换
ddct：离散余弦变换
ddst：离散正弦变换
dfct：RDFT（实对称DFT）的余弦变换
dfst：RDFT（实反对称DFT）的正弦变换

函数原型：

void cdft(int, int, double *, int *, double *);      // 复数离散傅立叶变换
void rdft(int, int, double *, int *, double *);      // 实离散傅里叶变换
void ddct(int, int, double *, int *, double *);      // 离散余弦变换
void ddst(int, int, double *, int *, double *);      // 离散正弦变换
void dfct(int, double *, double *, int *, double *); // 实对称DFT
void dfst(int, double *, double *, int *, double *); // 实反对称DFT

void cdft(int n, int isgn, double *a, int *ip, double *w)
{if (n > (ip[0] << 2)) {makewt(n >> 2, ip, w);}if (n > 4) {if (isgn >= 0) {bitrv2(n, ip + 2, a);cftfsub(n, a, w);} else {bitrv2conj(n, ip + 2, a);cftbsub(n, a, w);}} else if (n == 4) {cftfsub(n, a, w);}
}void rdft(int n, int isgn, double *a, int *ip, double *w)
{int nw, nc;double xi;nw = ip[0];if (n > (nw << 2)) {nw = n >> 2;makewt(nw, ip, w);}nc = ip[1];if (n > (nc << 2)) {nc = n >> 2;makect(nc, ip, w + nw);}if (isgn >= 0) {if (n > 4) {bitrv2(n, ip + 2, a);cftfsub(n, a, w);rftfsub(n, a, nc, w + nw);} else if (n == 4) {cftfsub(n, a, w);}xi = a[0] - a[1];a[0] += a[1];a[1] = xi;} else {a[1] = 0.5 * (a[0] - a[1]);a[0] -= a[1];if (n > 4) {rftbsub(n, a, nc, w + nw);bitrv2(n, ip + 2, a);cftbsub(n, a, w);} else if (n == 4) {cftfsub(n, a, w);}}
}void ddct(int n, int isgn, double *a, int *ip, double *w)
{int j, nw, nc;double xr;nw = ip[0];if (n > (nw << 2)) {nw = n >> 2;makewt(nw, ip, w);}nc = ip[1];if (n > nc) {nc = n;makect(nc, ip, w + nw);}if (isgn < 0) {xr = a[n - 1];for (j = n - 2; j >= 2; j -= 2) {a[j + 1] = a[j] - a[j - 1];a[j] += a[j - 1];}a[1] = a[0] - xr;a[0] += xr;if (n > 4) {rftbsub(n, a, nc, w + nw);bitrv2(n, ip + 2, a);cftbsub(n, a, w);} else if (n == 4) {cftfsub(n, a, w);}}dctsub(n, a, nc, w + nw);if (isgn >= 0) {if (n > 4) {bitrv2(n, ip + 2, a);cftfsub(n, a, w);rftfsub(n, a, nc, w + nw);} else if (n == 4) {cftfsub(n, a, w);}xr = a[0] - a[1];a[0] += a[1];for (j = 2; j < n; j += 2) {a[j - 1] = a[j] - a[j + 1];a[j] += a[j + 1];}a[n - 1] = xr;}
}void ddst(int n, int isgn, double *a, int *ip, double *w)
{int j, nw, nc;double xr;nw = ip[0];if (n > (nw << 2)) {nw = n >> 2;makewt(nw, ip, w);}nc = ip[1];if (n > nc) {nc = n;makect(nc, ip, w + nw);}if (isgn < 0) {xr = a[n - 1];for (j = n - 2; j >= 2; j -= 2) {a[j + 1] = -a[j] - a[j - 1];a[j] -= a[j - 1];}a[1] = a[0] + xr;a[0] -= xr;if (n > 4) {rftbsub(n, a, nc, w + nw);bitrv2(n, ip + 2, a);cftbsub(n, a, w);} else if (n == 4) {cftfsub(n, a, w);}}dstsub(n, a, nc, w + nw);if (isgn >= 0) {if (n > 4) {bitrv2(n, ip + 2, a);cftfsub(n, a, w);rftfsub(n, a, nc, w + nw);} else if (n == 4) {cftfsub(n, a, w);}xr = a[0] - a[1];a[0] += a[1];for (j = 2; j < n; j += 2) {a[j - 1] = -a[j] - a[j + 1];a[j] -= a[j + 1];}a[n - 1] = -xr;}
}void dfct(int n, double *a, double *t, int *ip, double *w)
{int j, k, l, m, mh, nw, nc;double xr, xi, yr, yi;nw = ip[0];if (n > (nw << 3)) {nw = n >> 3;makewt(nw, ip, w);}nc = ip[1];if (n > (nc << 1)) {nc = n >> 1;makect(nc, ip, w + nw);}m = n >> 1;yi = a[m];xi = a[0] + a[n];a[0] -= a[n];t[0] = xi - yi;t[m] = xi + yi;if (n > 2) {mh = m >> 1;for (j = 1; j < mh; j++) {k = m - j;xr = a[j] - a[n - j];xi = a[j] + a[n - j];yr = a[k] - a[n - k];yi = a[k] + a[n - k];a[j] = xr;a[k] = yr;t[j] = xi - yi;t[k] = xi + yi;}t[mh] = a[mh] + a[n - mh];a[mh] -= a[n - mh];dctsub(m, a, nc, w + nw);if (m > 4) {bitrv2(m, ip + 2, a);cftfsub(m, a, w);rftfsub(m, a, nc, w + nw);} else if (m == 4) {cftfsub(m, a, w);}a[n - 1] = a[0] - a[1];a[1] = a[0] + a[1];for (j = m - 2; j >= 2; j -= 2) {a[2 * j + 1] = a[j] + a[j + 1];a[2 * j - 1] = a[j] - a[j + 1];}l = 2;m = mh;while (m >= 2) {dctsub(m, t, nc, w + nw);if (m > 4) {bitrv2(m, ip + 2, t);cftfsub(m, t, w);rftfsub(m, t, nc, w + nw);} else if (m == 4) {cftfsub(m, t, w);}a[n - l] = t[0] - t[1];a[l] = t[0] + t[1];k = 0;for (j = 2; j < m; j += 2) {k += l << 2;a[k - l] = t[j] - t[j + 1];a[k + l] = t[j] + t[j + 1];}l <<= 1;mh = m >> 1;for (j = 0; j < mh; j++) {k = m - j;t[j] = t[m + k] - t[m + j];t[k] = t[m + k] + t[m + j];}t[mh] = t[m + mh];m = mh;}a[l] = t[0];a[n] = t[2] - t[1];a[0] = t[2] + t[1];} else {a[1] = a[0];a[2] = t[0];a[0] = t[1];}
}void dfst(int n, double *a, double *t, int *ip, double *w)
{int j, k, l, m, mh, nw, nc;double xr, xi, yr, yi;nw = ip[0];if (n > (nw << 3)) {nw = n >> 3;makewt(nw, ip, w);}nc = ip[1];if (n > (nc << 1)) {nc = n >> 1;makect(nc, ip, w + nw);}if (n > 2) {m = n >> 1;mh = m >> 1;for (j = 1; j < mh; j++) {k = m - j;xr = a[j] + a[n - j];xi = a[j] - a[n - j];yr = a[k] + a[n - k];yi = a[k] - a[n - k];a[j] = xr;a[k] = yr;t[j] = xi + yi;t[k] = xi - yi;}t[0] = a[mh] - a[n - mh];a[mh] += a[n - mh];a[0] = a[m];dstsub(m, a, nc, w + nw);if (m > 4) {bitrv2(m, ip + 2, a);cftfsub(m, a, w);rftfsub(m, a, nc, w + nw);} else if (m == 4) {cftfsub(m, a, w);}a[n - 1] = a[1] - a[0];a[1] = a[0] + a[1];for (j = m - 2; j >= 2; j -= 2) {a[2 * j + 1] = a[j] - a[j + 1];a[2 * j - 1] = -a[j] - a[j + 1];}l = 2;m = mh;while (m >= 2) {dstsub(m, t, nc, w + nw);if (m > 4) {bitrv2(m, ip + 2, t);cftfsub(m, t, w);rftfsub(m, t, nc, w + nw);} else if (m == 4) {cftfsub(m, t, w);}a[n - l] = t[1] - t[0];a[l] = t[0] + t[1];k = 0;for (j = 2; j < m; j += 2) {k += l << 2;a[k - l] = -t[j] - t[j + 1];a[k + l] = t[j] - t[j + 1];}l <<= 1;mh = m >> 1;for (j = 1; j < mh; j++) {k = m - j;t[j] = t[m + k] + t[m + j];t[k] = t[m + k] - t[m + j];}t[0] = t[m + mh];m = mh;}a[l] = t[0];}a[0] = 0;
}