Blame view

src/Parameters/fonctions/fourier/DiscreteFourierTransform.cc 5.26 KB
977f8ee1   Menouard AZIB   Creation of DFT
1
2
3
4
#include "DiscreteFourierTransform.hh"
#include <bits/stdc++.h>

template <class T, class E>
07bf1e7c   Menouard AZIB   EveryThing works ...
5
DiscreteFourierTransform<T, E>::DiscreteFourierTransform(vector<T> signal_, double sampleSpacing_)
977f8ee1   Menouard AZIB   Creation of DFT
6
{
977f8ee1   Menouard AZIB   Creation of DFT
7
    signal = signal_;
07bf1e7c   Menouard AZIB   EveryThing works ...
8
    sampleSpacing = sampleSpacing_;
977f8ee1   Menouard AZIB   Creation of DFT
9
10
11
}

template <class T, class E>
07bf1e7c   Menouard AZIB   EveryThing works ...
12
void DiscreteFourierTransform<T, E>::compute(bool computeFFT)
977f8ee1   Menouard AZIB   Creation of DFT
13
{
07bf1e7c   Menouard AZIB   EveryThing works ...
14
    if (computeFFT)
3cf11811   Menouard AZIB   FFT Validated and...
15
    {
07bf1e7c   Menouard AZIB   EveryThing works ...
16
17
18
19
20
21
        const int N = signal.size();
        const int powerOf2 = DiscreteFourierTransform<T, E>::highestPowerof2(N);
        // Zero padding is a simple concept; it simply refers to adding zeros to end of a time-domain signal to increase its length to reach power of 2.
        while (signal.size() < powerOf2)
            signal.push_back((T)0);

3cf11811   Menouard AZIB   FFT Validated and...
22
23
24
        phasors = DiscreteFourierTransform<T, E>::fft(signal);
    }
    else
3cf11811   Menouard AZIB   FFT Validated and...
25
        phasors = DiscreteFourierTransform<T, E>::dft(signal);
3cf11811   Menouard AZIB   FFT Validated and...
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
}

template <class T, class E>
std::vector<E> DiscreteFourierTransform<T, E>::createTestPoints(int N, double sampleSpacing)
{
    std::vector<E> out;
    const E t = (E)sampleSpacing;
    for (int i = 0; i < N; i++)
    {
        out.push_back((E)t * i);
    }
    return out;
}

template <class T, class E>
std::vector<E> DiscreteFourierTransform<T, E>::createTestFunction(std::vector<E> x)
{
    const int N = x.size();
    std::vector<E> out;
    for (int i = 0; i < N; i++)
    {
        const E val = (E)sin(50.0 * 2.0 * M_PI * x[i]) + 0.5 * sin(80.0 * 2.0 * M_PI * x[i]);
        out.push_back(val);
    }
    return out;
}

/**
 * @brief We use Cooley-Tukey FFT Algorithm if the size of the signal  is a power of 2 otherwise you should DFT brut. In The future we will use bluestein algorithm instead of DFT.
 * Cooley-Tukey FFT Algorithms: http://people.scs.carleton.ca/~maheshwa/courses/5703COMP/16Fall/FFT_Report.pdf
 *
 * @tparam T template type of input data
 * @tparam E template type of output data
 * @param sig the signal as an input
 * @return std::vector<std::complex<E>> the phasors (array of complex number)
 */
template <class T, class E>
std::vector<std::complex<E>> DiscreteFourierTransform<T, E>::fft(std::vector<T> sig)
{
    const int N = sig.size();
    if (N == 1)
977f8ee1   Menouard AZIB   Creation of DFT
67
    {
3cf11811   Menouard AZIB   FFT Validated and...
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
        std::vector<std::complex<E>> out;
        const std::complex<E> temp((E)sig[0], (E)0);
        out.push_back(temp);
        return out;
    }

    const std::complex<E> WN = (complex<E>)std::polar(1.0, 2 * M_PI / N);
    std::complex<E> W((E)1, (E)0);

    // divide and conquer:
    // Recurse: all even samples
    std::vector<std::complex<E>>
        x_evens = fft(getEven(sig));

    // Recurse: all odd samples
    std::vector<std::complex<E>> x_odds = fft(getOdd(sig));

    // Now, combine and perform N/2 operations!
    std::complex<E> zeroComplex((E)0, (E)0);
    std::vector<std::complex<E>> x(N, zeroComplex);
    for (int k = 0; k < N / 2; k++)
    {
        x[k] = x_evens[k] + W * x_odds[k];
        x[k + (N / 2)] = x_evens[k] - W * x_odds[k];
        W = W * WN;
    }
    return x;
}

template <class T, class E>
std::vector<std::complex<E>> DiscreteFourierTransform<T, E>::dft(std::vector<T> x)
{
    const int N = x.size();
    std::complex<E> zeroComplex((E)0, (E)0);
    std::vector<std::complex<E>> out(N, zeroComplex);
    for (int k = 0; k < N; k++)
    {
        for (int n = 0; n < N; n++)
977f8ee1   Menouard AZIB   Creation of DFT
106
        {
3cf11811   Menouard AZIB   FFT Validated and...
107
108
            const std::complex<E> expVal = (complex<E>)std::polar(1.0, -2 * k * n * M_PI / N);
            out[k] += ((E)x[n]) * expVal;
977f8ee1   Menouard AZIB   Creation of DFT
109
        }
3cf11811   Menouard AZIB   FFT Validated and...
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
    }
    return out;
}

template <class T, class E>
std::vector<E> DiscreteFourierTransform<T, E>::computeDSP(std::vector<std::complex<E>> x)
{
    const int N = x.size();
    std::vector<E> out;
    for (int k = 0; k < N / 2; k++)
    {
        const E magnitude = (E)pow(abs(x[k]), 1);
        const E dsp = (E)(2.0 / N) * magnitude;
        out.push_back(dsp);
    }
977f8ee1   Menouard AZIB   Creation of DFT
125

3cf11811   Menouard AZIB   FFT Validated and...
126
127
    return out;
}
977f8ee1   Menouard AZIB   Creation of DFT
128

3cf11811   Menouard AZIB   FFT Validated and...
129
template <class T, class E>
b57f2949   Menouard AZIB   Add comments
130
std::vector<E> DiscreteFourierTransform<T, E>::getFreq(std::vector<std::complex<E>> x, double frequency_)
3cf11811   Menouard AZIB   FFT Validated and...
131
132
{
    const int N = x.size();
3cf11811   Menouard AZIB   FFT Validated and...
133
134
135
    std::vector<E> out;
    for (int k = 0; k < N / 2; k++)
    {
b57f2949   Menouard AZIB   Add comments
136
        const E freq = (E)k * frequency_ / N;
3cf11811   Menouard AZIB   FFT Validated and...
137
        out.push_back(freq);
977f8ee1   Menouard AZIB   Creation of DFT
138
    }
3cf11811   Menouard AZIB   FFT Validated and...
139
140

    return out;
977f8ee1   Menouard AZIB   Creation of DFT
141
}
3cf11811   Menouard AZIB   FFT Validated and...
142
143

template <class T, class E>
b57f2949   Menouard AZIB   Add comments
144
std::vector<E> DiscreteFourierTransform<T, E>::getPeriods(std::vector<std::complex<E>> x, double frequency_)
ec9b5ba2   Menouard AZIB   Everything is wor...
145
146
{
    const int N = x.size();
b57f2949   Menouard AZIB   Add comments
147
    const double stepPeriod = N / sampleSpacing_;
ec9b5ba2   Menouard AZIB   Everything is wor...
148
149
150
    std::vector<E> out;
    for (int k = 0; k < N / 2; k++)
    {
b57f2949   Menouard AZIB   Add comments
151
        const E period = (E)k * frequency_;
ec9b5ba2   Menouard AZIB   Everything is wor...
152
153
154
155
156
157
158
        out.push_back(period);
    }

    return out;
}

template <class T, class E>
3cf11811   Menouard AZIB   FFT Validated and...
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
std::vector<T> DiscreteFourierTransform<T, E>::getOdd(std::vector<T> x)
{
    std::vector<T> odd;
    for (int i = 0; i < x.size(); i++)
    {
        if (i % 2 != 0)
            odd.push_back(x[i]);
    }
    return odd;
}

template <class T, class E>
std::vector<T> DiscreteFourierTransform<T, E>::getEven(std::vector<T> x)
{
    std::vector<T> even;
    for (int i = 0; i < x.size(); i++)
    {
        if (i % 2 == 0)
            even.push_back(x[i]);
    }
    return even;
}

template <class T, class E>
bool DiscreteFourierTransform<T, E>::isPowerOfTwo(int N)
{
    return (N & (N - 1)) == 0;
07bf1e7c   Menouard AZIB   EveryThing works ...
186
187
188
189
190
191
192
193
194
195
196
}

template <class T, class E>
int DiscreteFourierTransform<T, E>::highestPowerof2(int n)
{
    int p = n;
    while (!DiscreteFourierTransform<T, E>::isPowerOfTwo(p))
    {
        p += 1;
    }
    return p;
3cf11811   Menouard AZIB   FFT Validated and...
197
}