Last active
August 29, 2015 14:01
-
-
Save isqua/5744beda2d8d9e352976 to your computer and use it in GitHub Desktop.
Fourier Transform
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <cstdio> | |
#include <cstdlib> | |
#include <math.h> | |
using namespace std; | |
// Трансформируемая функция | |
float transforming_func(float x) { | |
float bottom = 1 + pow(x, 2.0); | |
return (1 / (M_PI * bottom)); | |
}; | |
// Середина отрезка | |
float get_middle(float t) { | |
return sqrt(fabs(t / M_PI - 1)); | |
}; | |
// Вычитаемое в модуле D | |
float get_first_interval_param(float t, float delta) { | |
return sqrt(t * exp(1) / delta - 1); | |
}; | |
// Интеграл от x_0 до x_n | |
float calculate_trapezium(float x_0, float x_n, int n) { | |
float sum = 0; // То, что в скобках в формуле трапеций | |
float dx = (x_n - x_0) / 2; // Длина маленького отрезочка | |
// Считаем сумму f(x_i) от (x_1 до x_n-1) | |
for (int i = 1; i < n; i++) { | |
sum += transforming_func(x_0 + dx * n); | |
} | |
// Прибавляем к ней полусумму f(x_0) и f(x_n) | |
sum += (transforming_func(x_0) + transforming_func(x_n)) / 2; | |
// Умножаем её на длину маленького отрезочка - это и есть искомый интеграл | |
return sum * dx; | |
}; | |
// Интеграл от x_0 до x_n | |
float calculate_simpson(float x_0, float x_n, int n) { | |
float sum = 0; | |
float coef = (x_n - x_0) / (2 * n); | |
float f_i; | |
float f[2 * n + 1]; | |
for (int i = 0; i <= 2 * n; i++) { | |
f_i = transforming_func(x_0 + i * coef); | |
f[i] = f_i; | |
if ((i >= 2) && (i <= (2 * n - 2))) { | |
sum += f_i; | |
} | |
} | |
sum = f[0] + 4 * (f[1] + sum + f[2 * n - 1]) + 2 * sum + f[2 * n]; | |
return sum * (x_n - x_0) / n / 6; | |
}; | |
float get_transform_x_i(float x_0, float x_n, int n) { | |
float dx = (x_n - x_0) / n; | |
float x_i = x_0; | |
for (int i = 1; i <= n; i++) { | |
printf("F(%.10e) = %.10e\n", x_i, transforming_func(x_i)); | |
x_i += dx; | |
} | |
}; | |
// Y(t) | |
float get_y_t(float t, float delta, int n) { | |
float x_m = get_middle(t); // Середина отрезка | |
float interval_half = fabs(get_first_interval_param(t, delta) - x_m); // Половина интервала | |
float x_0 = x_m - interval_half; | |
float x_n = x_m + interval_half; | |
float trapezium; | |
float simpson; | |
printf("x_m = %.8e\t", x_m); | |
printf("D = %.8e\t", interval_half); | |
printf("\nY1(t) = %.20e", calculate_trapezium(x_0, x_n, n)); | |
printf("\nY2(t) = %.20e\n", calculate_simpson(x_0, x_n, n)); | |
// Посчитать все f(x) на участке от x_0 до x_n: | |
// printf("\nY2(t) = %.20e\n", get_transform_x_i(x_0, x_n, n)); | |
}; | |
// Выводит все Y(t) для t = (10^(-k), ... , 10^(m)) | |
void get_y_t_on_interval(int k, int m, float delta) { | |
float t; | |
// Заменим первое число на отрицательное | |
k = -k; | |
if (k > m) { // Если первое число меньше второго | |
m ^= k ^= m ^= k; // Поменяем их местами | |
} | |
for (int i = k; i <= m; i++) { | |
printf("t = %.1e\t", t); | |
t = pow(10.0, i); | |
get_y_t(t, delta, 100); | |
printf("\n"); | |
} | |
}; | |
int main() | |
{ | |
int k = 3; | |
int m = 9; | |
float delta = 0.000000000000001; | |
get_y_t_on_interval(k, m, delta); | |
return 0; | |
}; |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Чтобы получить все f(xi) на участке от x0 до xn, надо раскомментировать строку про
закомментиовать вызов
get_y_t_on_interval
, а вместо него вызыватьget_y_t(t, delta, N)
вmain
с желаемымt