6 changed files with 204 additions and 15 deletions
@ -0,0 +1,51 @@ |
|||||
|
#include <fftw3.h> |
||||
|
#include <cmath> |
||||
|
#include <libplotcpp/plotcpp.hpp> |
||||
|
|
||||
|
double pi() { |
||||
|
return 3.1415; |
||||
|
} |
||||
|
|
||||
|
std::vector<double> tfd2vect(fftw_complex* tfd, int N) { |
||||
|
std::vector<double> res; |
||||
|
auto it = tfd; |
||||
|
for (int i = 0; i != N; ++i) { |
||||
|
fftw_complex c = {*it[0], *it[1]}; |
||||
|
res.push_back(sqrt(c[0]*c[0] + c[1]*c[1])); |
||||
|
it++; |
||||
|
} |
||||
|
|
||||
|
return res; |
||||
|
} |
||||
|
|
||||
|
int main(int argc, char** argv) { |
||||
|
QApplication app(argc, argv); |
||||
|
fftw_complex *in, *out; |
||||
|
fftw_plan p; |
||||
|
|
||||
|
int N = 500; |
||||
|
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); |
||||
|
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); |
||||
|
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_MEASURE); |
||||
|
|
||||
|
std::vector<double> xx; |
||||
|
for (int i = 0; i != N; ++i) { |
||||
|
xx.push_back(i); |
||||
|
} |
||||
|
|
||||
|
for (int i = 0; i != N; ++i) { |
||||
|
in[i][0] = sin(2*pi()*50*i/N); |
||||
|
} |
||||
|
|
||||
|
fftw_execute(p); /* repeat as needed */ |
||||
|
|
||||
|
std::vector<double> res = tfd2vect(out, N); |
||||
|
PlotCpp g; |
||||
|
g.plot(xx, res); |
||||
|
g.draw(); |
||||
|
|
||||
|
|
||||
|
fftw_destroy_plan(p); |
||||
|
fftw_free(in); fftw_free(out); |
||||
|
return app.exec(); |
||||
|
} |
||||
@ -0,0 +1,111 @@ |
|||||
|
#pragma once |
||||
|
|
||||
|
#include <vector> |
||||
|
#include <complex> |
||||
|
#include <opencv2/opencv.hpp> |
||||
|
#include <iterator> |
||||
|
#include <cmath> |
||||
|
|
||||
|
namespace math { |
||||
|
|
||||
|
using complex = std::complex<float>; |
||||
|
using signal = std::vector<double>; |
||||
|
using csignal = std::vector<complex>; |
||||
|
using contour = std::vector<cv::Point>; |
||||
|
constexpr double pi() {return std::atan(1)*4;} |
||||
|
|
||||
|
//TODO implémenter la fft
|
||||
|
csignal fft(const signal& input) { |
||||
|
//TODO: s'assurer que le signal est bien formé (i.e. bonne taille)
|
||||
|
return fft_rec(input); |
||||
|
}; |
||||
|
|
||||
|
csignal fft_rec(const signal& input) { |
||||
|
int size = input.size(); |
||||
|
|
||||
|
if (size == 1) { |
||||
|
return input; |
||||
|
} else { |
||||
|
signal odd; |
||||
|
signal even; |
||||
|
std::back_insert_iterator<signal> odd_back_it(odd); |
||||
|
std::back_insert_iterator<signal> even_back_it(even); |
||||
|
bool insert_in_even = false; |
||||
|
|
||||
|
for (auto it = input.begin(); it != input.end(); ++it) { |
||||
|
if (insert_in_even) { |
||||
|
*even_back_it++ = *it; |
||||
|
insert_in_even = false; |
||||
|
} else { |
||||
|
*odd_back_it++ = *it; |
||||
|
insert_in_even = true; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
signal odd_fft = fft_rec(odd); |
||||
|
signal even_fft = fft_rec(even); |
||||
|
signal res; |
||||
|
res.reserve(size); |
||||
|
|
||||
|
for (int k = 0; k<size/2; ++k) { |
||||
|
complex t = std::exp(complex(0, -2*pi()*k/size))*odd[k]; |
||||
|
res[k] = even[k] + t; |
||||
|
res[size/2+k] = even[k] - t; |
||||
|
} |
||||
|
return res; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
complex mean(const signal& sig) { |
||||
|
complex res = 0; |
||||
|
for (auto x: sig) { |
||||
|
res += x; |
||||
|
} |
||||
|
return complex(res.real()/sig.size(), res.imag()/sig.size()); |
||||
|
}; |
||||
|
|
||||
|
signal cont2sig(const contour& cont) { |
||||
|
signal sig; |
||||
|
auto sig_it = sig.begin(); |
||||
|
auto cont_it = cont.begin(); |
||||
|
|
||||
|
for (auto cont_it = cont.begin(); cont_it != cont.end(); ++cont_it) { |
||||
|
*sig_it = complex((*cont_it).x, (*cont_it).y); |
||||
|
sig_it++; |
||||
|
} |
||||
|
return sig; |
||||
|
}; |
||||
|
|
||||
|
contour coef2cont(const signal& tfd, complex mean, int size, int cmax) { |
||||
|
contour cont; |
||||
|
auto tf_it = tfd.begin(); |
||||
|
auto cont_it = cont.begin(); |
||||
|
|
||||
|
for (auto tf_it = tfd.begin(); tf_it != tfd.end(); ++tf_it) { |
||||
|
//TODO retrouver la formule
|
||||
|
//*cont_it = mean;
|
||||
|
} |
||||
|
return cont; |
||||
|
}; |
||||
|
|
||||
|
int max_cont(const std::vector<contour>& contours) { |
||||
|
int max = 0; |
||||
|
int id = 0; |
||||
|
for (int i=0; i<contours.size(); ++i) { |
||||
|
if (contours[i].size() > max) { |
||||
|
max = contours[i].size(); |
||||
|
id = i; |
||||
|
} |
||||
|
} |
||||
|
return id; |
||||
|
}; |
||||
|
|
||||
|
contour simplify_contour(const contour& cont, int cmax) { |
||||
|
contour res; |
||||
|
signal z = cont2sig(cont); |
||||
|
complex zm = mean(z); |
||||
|
signal tfd = fft(z); |
||||
|
res = coef2cont(tfd, zm, 0, cmax); |
||||
|
return res; |
||||
|
}; |
||||
|
} |
||||
@ -0,0 +1,13 @@ |
|||||
|
#include <math.hpp> |
||||
|
#include <cmath> |
||||
|
|
||||
|
int main(int argc, char** argv) { |
||||
|
math::signal s; |
||||
|
|
||||
|
for (int i=0; i<100; ++i) { |
||||
|
s.push_back(std::sin(2*math::pi()*50*i/100)); |
||||
|
} |
||||
|
math::fft(s); |
||||
|
|
||||
|
return 0; |
||||
|
} |
||||
Loading…
Reference in new issue