Fourier Transform

fourier_transform.cpp

#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;
};

Чтобы получить все f(x_i) на участке от x₀ до x_n, надо раскомментировать строку про

printf("\nY2(t) = %.20e\n", get_transform_x_i(x_0, x_n, n));

закомментиовать вызов get_y_t_on_interval, а вместо него вызывать get_y_t(t, delta, N) в main с желаемым t