init

gaoweidong

2026-01-13 c5eeada2c735b0209a061a233483c5cfc29c2230

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
#include <cmath>
#include "FFTHelper.h"
#include "const.h"
 
void FFTHelper::bitrp(double* xreal, double* ximag, int n)
{
    // Î»·´×ªÖû» Bit-reversal Permutation
    int i, j, a, b, p;
    for (i = 1, p = 0; i < n; i *= 2)
    {
        p++;
    }
    for (i = 0; i < n; i++)
    {
        a = i;
        b = 0;
        for (j = 0; j < p; j++)
        {
            b = b * 2 + a % 2;
            a = a / 2;
        }
        if (b > i)
        {
            double t = xreal[i];
            xreal[i] = xreal[b];
            xreal[b] = t;
            t = ximag[i];
            ximag[i] = ximag[b];
            ximag[b] = t;
        }
    }
}
 
 
int FFTHelper::FFT(double* xreal, double* ximag, int len)
{
    //nֵΪ2µÄN´Î·½
    int n = 2;
    while (n <= len)
    {
        n *= 2;
    }
    n /= 2;
    // ¿ìËÙ¸µÁ¢Ò¶±ä»»£¬½«¸´Êý x ±ä»»ºóÈÔ±£´æÔÚ x ÖУ¬xreal, ximag ·Ö±ðÊÇ x µÄʵ²¿ºÍÐ鲿
    double* wreal = new double[n / 2];
    double* wimag = new double[n / 2];
    double treal, timag, ureal, uimag, arg;
    int m, k, j, t, index1, index2;
 
    bitrp(xreal, ximag, n);
    // ¼ÆËã 1 µÄÇ° n / 2 ¸ö n ´Î·½¸ùµÄ¹²éÊý W'j = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
    arg = -2 * pi / n;
    treal = cos(arg);
    timag = sin(arg);
    wreal[0] = 1.0f;
    wimag[0] = 0.0f;
    for (j = 1; j < n / 2; j++)
    {
        wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
        wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
    }
    for (m = 2; m <= n; m *= 2)
    {
        for (k = 0; k < n; k += m)
        {
            for (j = 0; j < m / 2; j++)
            {
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;    // ÐýתÒò×Ó w µÄʵ²¿ÔÚ wreal [] ÖеÄϱêΪ t
                treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                ureal = xreal[index1];
                uimag = ximag[index1];
                xreal[index1] = ureal + treal;
                ximag[index1] = uimag + timag;
                xreal[index2] = ureal - treal;
                ximag[index2] = uimag - timag;
            }
        }
    }
    delete[] wreal;
    delete[] wimag;
    return n;
}
 
int FFTHelper::IFFT(double* xreal, double* ximag,int len)
{
    //nֵΪ2µÄN´Î·½
    int n = 2;
    while (n <= len)
    {
        n *= 2;
    }
    n /= 2;
    // ¿ìËÙ¸µÁ¢Ò¶Äæ±ä»»
    double* wreal = new double[n / 2];
    double* wimag = new double[n / 2];
    double treal, timag, ureal, uimag, arg;
    int m, k, j, t, index1, index2;
    bitrp(xreal, ximag, n);
    // ¼ÆËã 1 µÄÇ° n / 2 ¸ö n ´Î·½¸ù Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
    arg = 2 * pi / n;
    treal = cos(arg);
    timag = sin(arg);
    wreal[0] = 1.0;
    wimag[0] = 0.0;
    for (j = 1; j < n / 2; j++)
    {
        wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
        wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
    }
    for (m = 2; m <= n; m *= 2)
    {
        for (k = 0; k < n; k += m)
        {
            for (j = 0; j < m / 2; j++)
            {
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;    // ÐýתÒò×Ó w µÄʵ²¿ÔÚ wreal [] ÖеÄϱêΪ t
                treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                ureal = xreal[index1];
                uimag = ximag[index1];
                xreal[index1] = ureal + treal;
                ximag[index1] = uimag + timag;
                xreal[index2] = ureal - treal;
                ximag[index2] = uimag - timag;
            }
        }
    }
    for (j = 0; j < n; j++)
    {
        xreal[j] /= n;
        ximag[j] /= n;
    }
 
    delete[] wreal;
    delete[] wimag;
    return n;
}