【嵌入式经验积累】C语言实现傅里叶变换

文章目录

C 语言实现离散傅利叶变换（DFT）代码如下：
C 语言实现快速傅利叶变换（FFT），代码如下：
分裂基快速傅利叶变换代码如下。
实序列快速傅利叶变换（一），代码如下：
实序列快速傅利叶变换（二），代码如下：
Chirp Z 变换，代码如下：

算法笔记

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);
}