課程名稱：線性代數 CS233 B
授課教師：簡廷因
發放時間：2020-12-01
截止時間：2020-12-22 23:59:59

題目說明

給一些訊號，求做完傅立葉轉換後的結果。

解題思路

我完全不會，我是抄 這篇 的。

參考解法

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
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <cmath>

#define SIZE(c) int(c.size())
const double PI = 3.14159;

using namespace std;

// reference: https://ppt.cc/f1yLAx

static auto __ = []
{
    ios_base::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    return 0;
}();

struct signal
{
    double real;
    double image;
};

void split(string& str); // 資料分割
void bitrp();            // 位反轉置換 Bit-reversal Permutation
void FFT();

vector<signal> signals;

int main(int argc, char* argv[])
{
    int Case = 0;
    ifstream inFile(argv[1]);
    string str;
    while (getline(inFile, str))
    {
        signals.clear();
        split(str);
        FFT();
        
        cout << "Data " << ++Case << ":\n";
        for (auto& [r, i] : signals)
            cout << setprecision(4) << fixed << r << " " << i << '\n';
        cout << '\n';
    }
}

void split(string& str)
{
    stringstream ss(str);
    while (getline(ss, str, ',')) signals.push_back({ stod(str), 0 });
}

void bitrp()
{
    int p = 0;

    for (int i = 1; i < SIZE(signals); i *= 2) ++p;

    for (int i = 0; i < SIZE(signals); ++i)
    {
        int a = i, b = 0;
        for (int j = 0; j < p; ++j) b = b * 2 + a % 2, a /= 2;
        if (b > i) swap(signals[i], signals[b]);
    }
}

void FFT()
{
    vector<signal> w(SIZE(signals));

    bitrp();

    double arg = -2 * PI / SIZE(signals);
    double treal = cos(arg);
    double timage = sin(arg);
    w[0] = { 1.0, 0.0 };

    for (int i = 1; i < SIZE(signals) / 2; ++i)
    {
        w[i].real = w[i - 1].real * treal - w[i - 1].image * timage;
        w[i].image = w[i - 1].real * timage + w[i - 1].image * treal;
    }

    for (int m = 2; m <= SIZE(signals); m *= 2)
    {
        for (int k = 0; k < SIZE(signals); k += m)
        {
            for (int j = 0; j < m / 2; ++j)
            {
                int i1 = k + j;
                int i2 = i1 + m / 2;
                int t = SIZE(signals) * j / m;
                treal = w[t].real * signals[i2].real - w[t].image * signals[i2].image;
                timage = w[t].real * signals[i2].image + w[t].image * signals[i2].real;
                double ureal = signals[i1].real;
                double uimage = signals[i1].image;
                signals[i1] = { ureal + treal, uimage + timage };
                signals[i2] = { ureal - treal, uimage - timage };
            }
        }
    }
}

參考資料

快速傅立葉變換（FFT）的C++實現 - 开发者知识库