快速傅里叶变换
快速傅里叶变换
1、快速傅里叶变换
功能
计算复数序列的快速傅里叶变换
方法简介
序列\(x(n)\)的离散傅里叶变换为：
\[X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{nk},k=0,1,...N-1\]将序列\(x(n)\)按序号n的奇偶分为两组，即：
\[\begin{cases} x_1(n)=x(2n), & \text{if }n\text{ is even} \\ x_2(n)=x(2n+1), & \text{if }n\text{ is odd} \end{cases}\]因此，x(n)的傅里叶变换根据奇数偶数可以写成为：
\[X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{nk}+\sum_{n=0}^{N-1}x(n)W_{N}^{nk}\] \[=\sum_{n=0}^{N/2-1}x(2n)W_{N}^{2nk}+\sum_{n=0}^{N/2-1}x(2n+1)W_{N}^{(2n+1)k}\] \[=\sum_{n=0}^{N/2-1}x_1(n)W_{N/2}^{nk}+W_{N}^{k}\sum_{n=0}^{N/2-1}x_2(n)W_{N/2}^{nk}\]由此可得\(X(k)=X_1(k)+W_{N}^{k}X_2(k),k=0,1,...\frac{N}{2}-1\).
通过上面的分解，可以将一个N点的DFT分解为两个N/2点的DFT,每个N/2点的DFT又可以分为两个N/4点的DFT，依次类推，当N为2的整数次幂的时候\(N=2^M\)，由于每分解一次降低一阶幂次，所以通过M次的分解，最后全部成为一系列2点的DFT运算，也就是按时间抽取的快速傅里叶变化（FFT）算法。
快速傅里叶变换
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
/**
* 复数进行快速傅里叶变换
* x ： 长度为n的double数组，[in]存放要变换数据的实部,
* [out]存放变换完成数据的实部
* y ： 长度为n的double数组，[in]存放要变换数据的虚部,
* [out]存放变换完成数据的虚部
* n ： int数据长度必须是2的整数次幂,n=2^m
* sign ： 当sign = 1时,fft()计算的是离散傅里叶正变换(DFT),
* 当sign = -1 时fft()计算的是离散傅里叶逆变换(IDFT)
*/
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;
}
}
}
/*****************************
输入：
double x[8]={1,2,1,3,4,2,5,6};
double y[8]={0};
fft：
0:24.000000,0.000000
1:-0.878680,6.121320
2:-1.000000,5.000000
3:-5.121320,-1.878680
4:-2.000000,0.000000
5:-5.121320,1.878680
6:-1.000000,-5.000000
7:-0.878680,-6.121320
ifft：
0:1.000000,0.000000
1:2.000000,0.000000
2:1.000000,-0.000000
3:3.000000,-0.000000
4:4.000000,0.000000
5:2.000000,0.000000
6:5.000000,0.000000
7:6.000000,0.000000
*****************************/
其他：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#ifndef _FFT_H
#define _FFT_H
#define PI 3.14159265358979
// #define pi 3.14159265358979
typedef struct
{
float real;
float imag;
} Complex;
enum
{
FFT,
IFFT
};
int fft(Complex *x, int N);
int ifft(Complex *x, int N);
#endif // _FFT_H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include "fft.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/*
* FFT Algorithm
* === Inputs ===
* x : complex numbers
* N : nodes of FFT. for example:64,128,256,512,1024,etc.
* === Output ===
* the @x contains the result of FFT algorithm, so the original data
* in @x is destroyed, please store them before using FFT.
*/
int fft(Complex *x, int N)
{
int i, j, l, k, ip;
static int M = 0;
static int le, le2;
static float sR, sI, tR, tI, uR, uI;
M = (int)(log(N) / log(2));
/*
* bit reversal sorting
*/
l = N / 2;
j = l;
// BitReverse(x,x,N,M);
for (i = 1; i <= N - 2; i++)
{
if (i < j)
{
tR = x[j].real;
tI = x[j].imag;
x[j].real = x[i].real;
x[j].imag = x[i].imag;
x[i].real = tR;
x[i].imag = tI;
}
k = l;
while (k <= j)
{
j = j - k;
k = k / 2;
}
j = j + k;
}
/*
* For Loops
*/
for (l = 1; l <= M; l++)
{ /* loop for ceil{log2(N)} */
le = (int)pow(2, l);
le2 = (int)(le / 2);
uR = 1;
uI = 0;
sR = cos(PI / le2);
sI = -sin(PI / le2);
for (j = 1; j <= le2; j++)
{ /* loop for each sub DFT */
// jm1 = j - 1;
for (i = j - 1; i <= N - 1; i += le)
{ /* loop for each butterfly */
ip = i + le2;
tR = x[ip].real * uR - x[ip].imag * uI;
tI = x[ip].real * uI + x[ip].imag * uR;
x[ip].real = x[i].real - tR;
x[ip].imag = x[i].imag - tI;
x[i].real += tR;
x[i].imag += tI;
} /* Next i */
tR = uR;
uR = tR * sR - uI * sI;
uI = tR * sI + uI * sR;
} /* Next j */
} /* Next l */
return 0;
}
/*
* Inverse FFT Algorithm
* === Inputs ===
* x : complex numbers
* N : nodes of FFT. @N should be power of 2, that is 2^(*)
* === Output ===
* the @x contains the result of FFT algorithm, so the original data
* in @x is destroyed, please store them before using FFT.
*/
int ifft(Complex *x, int N)
{
int k = 0;
for (k = 0; k <= N - 1; k++)
{
x[k].imag = -x[k].imag;
}
fft(x, N); /* using FFT */
for (k = 0; k <= N - 1; k++)
{
x[k].real = x[k].real / N;
x[k].imag = -x[k].imag / N;
}
return 0;
}
2、基4快速傅里叶变换
若\(N=4^M\)，则将序列x(n)分为四个N/4点的序列，这样，就将一个N点的DFT转换为四个N/4点的DFT来计算。依此类推，直至分解到最后一级。以上就是按频率抽取的基4快速傅立叶变换算法。与基2FFT相比，基4FFT的乘法量约减少25%，加法量也略有减少。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
/****error:该函数计算存在错误****/
/**
* 基4复数快速傅里叶变换
* x ： 长度为n的double数组，[in]存放要变换数据的实部,
* [out]存放变换完成数据的实部
* y ： 长度为n的double数组，[in]存放要变换数据的虚部,
* [out]存放变换完成数据的虚部
* n ： int数据长度必须是4的整数次幂,n=4^m
*/
void fft4(double *x, double *y, int n)
{
int i, j, k, m, i1, i2, i3, n1, n2;
double a, b, c, e, r1, r2, r3, r4, s1, s2, s3, s4;
double co1, co2, co3, si1, si2, si3;
for (j = 1, i = 1; i < 10; i++)
{
m = i;
j = 4 * j;
if (j == n)
{
break;
}
}
n2 = n;
for (k = 1; k <= m; k++)
{
n1 = n2;
n2 = n2 / 4;
e = 6.28318530718 / n1;
a = 0;
for (j = 0; j < n2; j++)
{
b = a + a;
c = a + b;
co1 = cos(a);
co2 = cos(b);
co3 = cos(c);
si1 = sin(a);
si2 = sin(b);
si3 = sin(c);
a = (j + 1) * e;
for (i = j; i < n; i = i + n1)
{
i1 = i + n2;
i2 = i1 + n2;
i3 = i2 + n2;
r1 = x[i] + x[i2];
r3 = x[i] - x[i2];
s1 = y[i] + y[i2];
s3 = y[i] - y[i2];
r2 = x[i1] + x[i3];
r4 = x[i1] - x[i3];
s2 = y[i1] + y[i3];
s4 = y[i1] - y[i3];
x[i] = r1 - r2;
r2 = r1 - r2;
r1 = r3 - s4;
r3 = r3 + s4;
y[i] = s1 + s2;
s2 = s1 - s2;
s1 = s3 + r4;
s3 = s3 - r4;
x[i1] = co1 * r3 + si1 * s3;
y[i1] = co1 * s3 - si1 * r3;
x[i2] = co2 * r2 + si2 * s2;
y[i2] = co2 * s2 - si2 * r2;
x[i3] = co3 * r1 + si3 * s1;
y[i3] = co3 * s1 - si3 * r1;
}
}
}
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 / 4;
while (3 * k < (j + 1))
{
j = j - 3 * k;
k = k / 4;
}
j = j + k;
}
}
其他：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#ifndef _FFT_H
#define _FFT_H
#define PI 3.14159265358979
#define pi 3.14159265358979
typedef struct
{
float real;
float imag;
} Complex;
enum
{
FFT,
IFFT
};
void fft4(Complex *in, int log4_N);
void ifft4(Complex *in, int log4_N);
#endif // _FFT_H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
#include "fft.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
int BitReverse(int src, int size)
{
int tmp = src;
int des = 0;
int i;
for (i = size - 1; i >= 0; i--)
{
des = ((tmp & 0x3) << (i * 2)) | des;
tmp = tmp >> 2;
}
return des;
}
void reverse_idx(Complex *in, int log4_N)
{
int i;
int N = 1 << (log4_N * 2);
Complex *temp;
temp = malloc(N * sizeof(Complex));
if (!temp)
{
printf("malloc failed!\n");
exit(0);
}
for (i = 0; i < N; i++)
{
int idx;
idx = BitReverse(i, log4_N);
temp[idx].real = in[i].real;
temp[idx].imag = in[i].imag;
}
for (i = 0; i < N; i++)
{
in[i].real = temp[i].real;
in[i].imag = temp[i].imag;
}
free(temp);
}
void fft_ifft_4_common(Complex *in, Complex *win, int log4_N, int reverse)
{
int N = (1 << log4_N * 2);
int i, j, k;
int span = 1;
int n = N >> 2;
int widx;
Complex temp1, temp2, temp3, temp4;
int idx1, idx2, idx3, idx4;
for (i = 0; i < log4_N; i++)
{
for (j = 0; j < n; j++)
{
widx = 0;
idx1 = j * span * 4;
idx2 = idx1 + span;
idx3 = idx2 + span;
idx4 = idx3 + span;
for (k = 0; k < span; k++)
{
temp1.real = in[idx1 + k].real;
temp1.imag = in[idx1 + k].imag;
temp2.real = win[widx].real * in[idx2 + k].real - win[widx].imag * in[idx2 + k].imag;
temp2.imag = win[widx].imag * in[idx2 + k].real + win[widx].real * in[idx2 + k].imag;
temp3.real = win[widx * 2].real * in[idx3 + k].real - win[widx * 2].imag * in[idx3 + k].imag;
temp3.imag = win[widx * 2].imag * in[idx3 + k].real + win[widx * 2].real * in[idx3 + k].imag;
temp4.real = win[widx * 3].real * in[idx4 + k].real - win[widx * 3].imag * in[idx4 + k].imag;
temp4.imag = win[widx * 3].imag * in[idx4 + k].real + win[widx * 3].real * in[idx4 + k].imag;
in[idx1 + k].real = temp1.real + temp3.real;
in[idx1 + k].imag = temp1.imag + temp3.imag;
in[idx2 + k].real = temp1.real - temp3.real;
in[idx2 + k].imag = temp1.imag - temp3.imag;
in[idx3 + k].real = temp2.real + temp4.real;
in[idx3 + k].imag = temp2.imag + temp4.imag;
in[idx4 + k].real = temp2.real - temp4.real;
in[idx4 + k].imag = temp2.imag - temp4.imag;
temp1.real = in[idx1 + k].real + in[idx3 + k].real;
temp1.imag = in[idx1 + k].imag + in[idx3 + k].imag;
if (reverse == 0)
{
temp2.real = in[idx2 + k].real + in[idx4 + k].imag;
temp2.imag = in[idx2 + k].imag - in[idx4 + k].real;
}
else
{
temp2.real = in[idx2 + k].real - in[idx4 + k].imag;
temp2.imag = in[idx2 + k].imag + in[idx4 + k].real;
}
temp3.real = in[idx1 + k].real - in[idx3 + k].real;
temp3.imag = in[idx1 + k].imag - in[idx3 + k].imag;
if (reverse == 0)
{
temp4.real = in[idx2 + k].real - in[idx4 + k].imag;
temp4.imag = in[idx2 + k].imag + in[idx4 + k].real;
}
else
{
temp4.real = in[idx2 + k].real + in[idx4 + k].imag;
temp4.imag = in[idx2 + k].imag - in[idx4 + k].real;
}
in[idx1 + k].real = temp1.real;
in[idx1 + k].imag = temp1.imag;
in[idx2 + k].real = temp2.real;
in[idx2 + k].imag = temp2.imag;
in[idx3 + k].real = temp3.real;
in[idx3 + k].imag = temp3.imag;
in[idx4 + k].real = temp4.real;
in[idx4 + k].imag = temp4.imag;
widx += n;
}
}
n >>= 2;
span <<= 2;
}
}
/**
* @brief
*
* @param in :complex number input
* @param log4_N :64=3,256=4,1024=5
*/
void fft4(Complex *in, int log4_N)
{
Complex *win;
int N = 1 << (log4_N * 2);
int i;
win = malloc((3 * N / 4 - 2) * sizeof(Complex));
if (!win)
{
printf("malloc failed!\n");
exit(0);
}
reverse_idx(in, log4_N);
for (i = 0; i < (3 * N / 4 - 2); i++)
{
win[i].real = cos(2 * pi * i / (float)N);
win[i].imag = -sin(2 * pi * i / (float)N);
}
fft_ifft_4_common(in, win, log4_N, 0);
free(win);
}
void ifft4(Complex *in, int log4_N)
{
int N = 1 << (log4_N * 2);
Complex *win;
int i;
win = malloc((3 * N / 4 - 2) * sizeof(Complex));
if (!win)
{
printf("malloc failed!\n");
exit(0);
}
reverse_idx(in, log4_N);
for (i = 0; i < (3 * N / 4 - 2); i++)
{
win[i].real = cos(2 * pi * i / (float)N);
win[i].imag = sin(2 * pi * i / (float)N);
}
fft_ifft_4_common(in, win, log4_N, 1);
for (i = 0; i < N; i++)
{
in[i].real /= (float)N;
in[i].imag /= (float)N;
}
free(win);
}
3、分裂基快速傅里叶变换
将X(k)按序号k的奇偶分为两组，当k为偶数的时候，进行基2频率抽取分解，当k为奇数的时候，进行基4频率抽取分解，这样之后，一个N 点的DFT可以分解为一个N/2点的DFT和两个N/4 点的DFT。这种分解既有基2的部分，又有基4的部分，因此称为分裂基分解。上面的N/2点DFT又可分解为一个N/4点的DFT和两个N/8点的DFT，而两个N/4点的DFT也分别可以分解为-个N/8点的DFT和两个N/16点的DFT。依此类推，直至分解到最后一级为止。这就是按频率抽取的分裂基快速傅立叶变换算法。。 分裂基快速算法是将基2和基4分解组合而成。在基2类快速算法中，分裂基算法具有最少的运算量，且仍保留结构规则、原位计算等优点。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
/**
* 分裂基快速傅里叶变换
* x ： 长度为n的double数组，[in]存放要变换数据的实部,
* [out]存放变换完成数据的实部
* y ： 长度为n的double数组，[in]存放要变换数据的虚部,
* [out]存放变换完成数据的虚部
* n ： int数据长度必须是2的整数次幂,n=2^m
*/
void srfft(double *x, double *y, int n)
{
int i, j, k, m, i1, i2, i3, n1, n2, n4, id, is;
double a, b, c, e, a3, r1, r2, r3, r4;
double c1, c3, 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 = (j + 1) * 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;
}
}
/*****
输入：
double x2[16]={1,2,1,3,2,5,6,3,7,8,2,4,5,8,3,2};
double y2[16]={0};
输出：
rcfft:
0:62.000000,0.000000
1:-14.530217,7.194722
2:-2.535534,6.707107
3:-7.698116,-1.325550
4:3.000000,-11.000000
5:1.354970,7.502877
6:4.535534,-5.292893
7:-3.126637,4.023149
8:-8.000000,0.000000
9:-3.126637,-4.023149
10:4.535534,5.292893
11:1.354970,-7.502877
12:3.000000,11.000000
13:-7.698116,1.325550
14:-2.535534,-6.707107
15:-14.530217,-7.194722
*****/
4、实数快速傅里叶变换
上式可用复序列FFT算法进行计算。但考虑到x(n)是实数，为进一步提高计算效率，需要对按时间抽取的基2复序列FFT算法进行一定的修改。 如果序列x(n)是实数，那么其傅立叶变换X(k)一般是复数，但其实部是偶对称，虚部是奇对称，即X(k)具有如下共轭对称性：X(0)和X(N/2)都是实数，且有： \(X(k)=X(N-k),1\leq k \leq \frac{N}{2}-1\) 在计算离散傅立叶变换时，利用这种共轭对称性，我们就可以不必计算与存储X(k) (N/2 +1≤k≤N-1 )以及X(0)和X(N/2)的虚部，而仅需计算X(0)到X(N/2)即可。此处我们选择的是计算X(0)到X(N/4)和X(N/2)到X(3N/4)，这样做可以恰好利用复序列FFT算法的前(N/4)+1个复数蝶形。这就是按时间抽取的基2实序列FFT算法，它比复序列FFT算法大约可减少一半的运算量和存储量。
针对在单片机上运行的代码，这些FFT都不快，只能运行基于查表法或汇编语言的FFT库
正变换：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
/**
* 实数傅里叶变换
* x ： 长度为n的double数组，[in]存放要变换数据的实部,
* [out]存放变换完成数据的实部和虚部，顺序为：[Re(0),Re(1),...,Re(n/2),Im(n/2-1),...,Im(1)]
* n ： int数据长度必须是2的整数次幂,n=2^m
*/
void rfft(double *x, int 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;
}
}
}
}
/*****
输入：
double x2[16]={1,2,1,3,2,5,6,3,7,8,2,4,5,8,3,2};
输出：
rfft
0:62.000000
1:-14.530217
2:-2.535534
3:-7.698116
4:3.000000
5:1.354970
6:4.535534
7:-3.126637
8:-8.000000
9:4.023149
10:-5.292893
11:7.502877
12:-11.000000
13:-1.325550
14:6.707107
15:7.194722
*****/
反变换：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
/**
* 实数傅里叶逆变换
* 开始时存放具有共轭对称性的复序列X的前n/2+1个值,
* x ： 长度为n的double数组，[in]存放要变换数据的实部,
* 顺序为：[Re(0),Re(1),...,Re(n/2),Im(n/2-1),...,Im(1)]
* [out]存放变换完成数据的实部，x(i),(i=0,1,...,n-1)
* n ：int数据长度必须是2的整数次幂,n=2^m
*/
void irfft(double *x, int n)
{
int i, j, k, m, i1, i2, i3, i4, i5, i6, i7, i8, n2, n4, n8, id, is;
double a, e, a3, t1, t2, t3, t4, t5, 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++)
{
is = 0;
id = n2;
n2 = n2 / 2;
n4 = n2 / 4;
n8 = n4 / 2;
e = 6.28318530718 / n2;
do
{
for (i = is; i < n; i += id)
{
i1 = i;
i2 = i1 + n4;
i3 = i2 + n4;
i4 = i3 + n4;
t1 = x[i1] - x[i3];
x[i1] = x[i1] + x[i3];
x[i2] = 2 * x[i2];
x[i3] = t1 - 2 * x[i4];
x[i4] = t1 + 2 * x[i4];
if (n4 == 1)
continue;
i1 += n8;
i2 += n8;
i3 += n8;
i4 += n8;
t1 = (x[i2] - x[i1]) / sqrt(2.0);
t2 = (x[i4] + x[i3]) / sqrt(2.0);
x[i1] = x[i1] + x[i2];
x[i2] = x[i4] - x[i3];
x[i3] = 2 * (-t2 - t1);
x[i4] = 2 * (-t2 + t1);
}
is = 2 * id - n2;
id = 4 * id;
} while (is < (n - 1));
a = e;
for (j = 1; j < n8; j++)
{
a3 = 3 * a;
cc1 = cos(a);
ss1 = sin(a);
cc3 = cos(a3);
ss3 = sin(a3);
a = (j + 1) * e;
is = 0;
id = 2 * n2;
do
{
for (i = is; i <= (n - 1); i = i + id)
{
i1 = i + j;
i2 = i1 + n4;
i3 = i2 + n4;
i4 = i3 + n4;
i5 = i + n4 - j;
i6 = i5 + n4;
i7 = i6 + n4;
i8 = i7 + n4;
t1 = x[i1] - x[i6];
x[i1] = x[i1] + x[i6];
t2 = x[i5] - x[i2];
x[i5] = x[i2] + x[i5];
t3 = x[i8] + x[i3];
x[i6] = x[i8] - x[i3];
t4 = x[i4] + x[i7];
x[i2] = x[i4] - x[i7];
t5 = t1 - t4;
t1 = t1 + t4;
t4 = t2 - t3;
t2 = t2 + t3;
x[i3] = t5 * cc1 + t4 * ss1;
x[i7] = -t4 * cc1 + t5 * ss1;
x[i4] = t1 * cc3 - t2 * ss3;
x[i8] = t2 * cc3 + t1 * ss3;
}
is = 2 * id - n2;
id = 4 * id;
} while (is < (n - 1));
}
}
is = 0;
id = 4;
do
{
for (i = is; i < n; i = i + id)
{
i1 = i + 1;
t1 = x[i];
x[i] = t1 + x[i1];
x[i1] = t1 - x[i1];
}
is = 2 * id - 2;
id = 4 * id;
} while (is < (n - 1));
for (j = 0, i = 0; i < (n - 1); i++)
{
if (i < j)
{
t1 = x[j];
x[j] = x[i];
x[i] = t1;
}
k = n / 2;
while (k < (j + 1))
{
j = j - k;
k = k / 2;
}
j = j + k;
}
for (i = 0; i < n; i++)
{
x[i] = x[i] / n;
}
}
5、Chirp-Z小片段高精度计算
在Z平面单位圆上计算有限长序列x(n)的Z变换的采样值。序列 x(n)的 DFT 实质上是其 Z-变换 X(z)在 Z 平面单位圆上等间隔采样点\(z_k=\frac{2\pi}{N}\),k处的采样值。但实际应用中，往往只需要对信号的一小段频带进行高分辨率分析，这时通常使用Chirp Z-变换算法。
设x(n)是长度为N的序列，其傅立叶变换为\(X(e^{j\omega})\)。在单位圆任意一段弧上等间隔地取M个频率点 \(\omega _k =\omega _0 + k \Delta \omega\) 式中\(\omega _0\)为起始频率，\(\Delta \omega\)为频率增量。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
/**
* chirp-z变换算法，对一小段频带进行高分辨率分析
* xr---double型数组，长度大于等于(n+m-1),且是2的整数次幂,
* [in]:输入数据的实部,[out]:输出数据的实部
* xi---double型数组，长度大于等于(n+m-1),且是2的整数次幂,
* [in]:输入数据的虚部,[out]:输出数据的虚部
* n---整形变量，输入数据的长度。
* m---整形变量，输出信号的长度，即频率采样点数
* f1---double变量，起始数字频率，单位为Hz-s
* f2---double变量，终止数字频率，单位为Hz-s
*/
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 = M_PI;
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);
}
测试：
假如信号的采样率为16khz，一个正弦波信号为1200hz，采样128个点的数据，做128个点的FFT，则有：
1
2
3
4
5
6
7
8
delta_f=16000/128=125hz
1200/125=9.6
正常的fft无法精确的计算出1200hz处的能量强度
结果如下：
1000.00, 13.24
1125.00, 33.36
1250.00, 47.42
1375.00, 12.87
对该信号做czt处理
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#define N (128)
#define M (256)
int main(int argc, char **argv)
{
double x[M] = {0}, y[M] = {0}, z[M] = {0};
int i = 0;
FILE *fp = fopen("fft.txt", "w");
// 初始化数据x
double fs = 16000;
double f = 1200;
for (i = 0; i < N; i++)
{
x[i] = sin(2 * M_PI * f * i / fs);
}
fft(x, y, N, 1);
user_abs(x, y, z, N);
for (i = 0; i < N / 2; i++)
{
fprintf(fp, "%8.2f,%8.2f\n", i * (fs / N), z[i]);
}
fclose(fp);
fp = fopen("czt.txt", "w");
memset(x, 0, M * sizeof(double));
memset(y, 0, M * sizeof(double));
memset(z, 0, M * sizeof(double));
for (i = 0; i < N; i++)
{
x[i] = sin(2 * M_PI * f * i / fs);
}
/*精确计算数字频率0.05~0.1之间的频点*/
double f1 = 0.05;
double f2 = 0.1;
double xf;
// 这里的有效输出结果为N个点，也就是函数的入口参数m
czt(x, y, N, N, f1, f2);
user_abs(x, y, z, N);
for (i = 0; i < N; i++)
{
xf = (f1 + i * (f2 - f1) / (N - 1)) * fs;
fprintf(fp, "%8.2f,%8.2f\n", xf, z[i]);
}
fclose(fp);
printf("process done\n");
return 0;
}
对比fft结果和czt结果：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
FFT结果：
750.00, 6.65
875.00, 8.64
1000.00, 13.24
1125.00, 33.36
1250.00, 47.42
1375.00, 12.87
1500.00, 7.15
1625.00, 4.82
1750.00, 3.56
CZT结果：
1102.36, 17.65
1108.66, 22.00
1114.96, 26.39
1121.26, 30.78
1127.56, 35.10
1133.86, 39.30
1140.16, 43.31
1146.46, 47.09
1152.76, 50.58
1159.06, 53.72
1165.35, 56.48
1171.65, 58.81
1177.95, 60.69
1184.25, 62.09
1190.55, 62.98
1196.85, 63.36
1203.15, 63.23
1209.45, 62.57
1215.75, 61.42
1222.05, 59.78
1228.35, 57.67
1234.65, 55.14
1240.94, 52.22
1247.24, 48.95
1253.54, 45.38
1259.84, 41.56
1266.14, 37.54
1272.44, 33.39
1278.74, 29.15
1285.04, 24.88
1291.34, 20.65
1297.64, 16.49
1303.94, 12.48
本文由作者按照 CC BY 4.0 进行授权