diff --git a/C/FP_math b/C/FP_math new file mode 100755 index 0000000..b8efcff Binary files /dev/null and b/C/FP_math differ diff --git a/C/FP_math.c b/C/FP_math.c new file mode 100644 index 0000000..2b28509 --- /dev/null +++ b/C/FP_math.c @@ -0,0 +1,152 @@ +#include +#include +#include +#include +#include + +#define DATA_L 1024 +#define FP_acc 1000 + +#define PI 3.141592653589793238462643383279502884197169399375105820974944 + + + +int cos_delta_FP_arr[DATA_L] = {0,}; +int sin_delta_FP_arr[DATA_L] = {0,}; + +void FFT_fp(uint32_t* inp, uint32_t inp_L, uint32_t* buf){ + + // buf имеет длину inp_L * 2 (Re, Im, Re, Im, ...) + // inp содержит inp_L значений uint32_t + + uint32_t i, j, bit; +// uint32_t N = inp_L / 2; // длина комплексного массива (inp содержит только Re-часть) + uint32_t N = inp_L; // длина комплексного массива (inp содержит только Re-часть) + + // --- копирование входных данных в буфер --- + for (i = 0; i < inp_L; i++) { + buf[i * 2] = (int)inp[i]; // Re + buf[i * 2 + 1] = 0; // Im = 0 + } + + // --- bit-reversal перестановка --- + j = 0; + for (i = 1; i < N; i++) { + bit = N >> 1; + while (j & bit) { + j ^= bit; + bit >>= 1; + } + j |= bit; + if (i < j) { + int tmp_re = buf[i * 2]; + int tmp_im = buf[i * 2 + 1]; + buf[i * 2] = buf[j * 2]; + buf[i * 2 + 1] = buf[j * 2 + 1]; + buf[j * 2] = tmp_re; + buf[j * 2 + 1] = tmp_im; + } + } + + // --- уровни бабочек --- + uint32_t m = 2; + while (m <= N) { + uint32_t half = m >> 1; + + int c_delta_FP = cos_delta_FP_arr[m - 2]; + int s_delta_FP = sin_delta_FP_arr[m - 2]; + + for (uint32_t start = 0; start < N; start += m) { + int wr = FP_acc; + int wi = 0; + + for (uint32_t k = 0; k < half; k++) { + int u_re = buf[(start + k) * 2]; + int u_im = buf[(start + k) * 2 + 1]; + int v_re = buf[(start + k + half) * 2]; + int v_im = buf[(start + k + half) * 2 + 1]; + + // t = w * v (в фиксированной точке) + int64_t t_re = ((int64_t)wr * v_re - (int64_t)wi * v_im) / FP_acc; + int64_t t_im = ((int64_t)wr * v_im + (int64_t)wi * v_re) / FP_acc; + + // верх/низ + buf[(start + k) * 2] = u_re + (int)t_re; + buf[(start + k) * 2 + 1] = u_im + (int)t_im; + buf[(start + k + half) * 2] = u_re - (int)t_re; + buf[(start + k + half) * 2 + 1] = u_im - (int)t_im; + + // w *= w_m (обновление twiddle) + int64_t wr_new = ((int64_t)wr * c_delta_FP - (int64_t)wi * s_delta_FP) / FP_acc; + int64_t wi_new = ((int64_t)wr * s_delta_FP + (int64_t)wi * c_delta_FP) / FP_acc; + wr = (int)wr_new; + wi = (int)wi_new; + } + } + m <<= 1; + } + for (uint32_t i = 0; i < inp_L; ++i){ + //buf[i*2] /= inp_L; + //buf[i*2+1] /= inp_L; + printf("re,im: %d, %d\n", buf[i*2], buf[i*2 +1]); + } + +} + +void imagery2abs(uint32_t* inp, uint32_t inp_L, uint32_t *out){ + printf("+++++++++++++++++++++\n"); + printf("calculating abs^2\n"); + for (uint32_t i = 0; i < inp_L; ++i){ + int val = (inp[i*2]*inp[i*2] + inp[i*2 +1]*inp[i*2 +1])/FP_acc; + printf("%d\n", val); + out[i] = val; + } +} + + + +void sin_cos_filler(int* sin_arr, int* cos_arr, uint32_t L){ + //[int(-sin(2.0 * pi / m) * FP_acc) for m in range(2,inp_L//2 + 2)] + for (uint32_t i = 0; i < L; ++i){ + int m = i + 2; + double angle = 2.0* PI/m; + sin_arr[i] = lround(-sin(angle) * FP_acc); + cos_arr[i] = lround(cos(angle) * FP_acc); + printf("i, angle, sin, cos: %d %g %d %d\n", i, angle, sin_arr[i], cos_arr[i]); + } +} + +void data_generator(uint32_t* X, uint32_t* Y, uint32_t N){ + for (int i = 0; i < N; ++i){ + X[i] = (int) i; + Y[i] = lround(sin(((double)i)*2.0* PI/(N/10))*FP_acc); +// Y[I] += 10*FP_acc; + } +} + +void main(){ + char* logfilename = "tmp"; + uint32_t data_X[DATA_L] = {0,}; + uint32_t data_Y[DATA_L] = {0,}; + uint32_t data_res_FFT_imag[DATA_L*2] = {0,}; + uint32_t data_res_FFT_abs[DATA_L] = {0,}; + + sin_cos_filler(sin_delta_FP_arr, cos_delta_FP_arr, DATA_L); + + data_generator(data_X, data_Y, DATA_L); + + + FFT_fp(data_Y, DATA_L, data_res_FFT_imag); + + + imagery2abs(data_res_FFT_imag, DATA_L, data_res_FFT_abs); + + FILE *logfile_ptr; + logfile_ptr = fopen(logfilename, "w"); + fprintf(logfile_ptr, "X, Y, FFT\n"); + for (uint32_t i = 0; i < DATA_L; ++i){ + fprintf(logfile_ptr, "%d, %d, %d \n", data_X[i], data_Y[i],data_res_FFT_abs[i]); + } + fclose(logfile_ptr); +} + diff --git a/C/Makefile b/C/Makefile new file mode 100644 index 0000000..1047e84 --- /dev/null +++ b/C/Makefile @@ -0,0 +1,5 @@ +all: + gcc FP_math.c -o FP_math -lm + +clean: + -rm FP_math diff --git a/C/plotter.py b/C/plotter.py new file mode 100755 index 0000000..509986e --- /dev/null +++ b/C/plotter.py @@ -0,0 +1,56 @@ +#!/usr/bin/pypy3 +import plotly.graph_objs as go +from decimal import * +from sys import argv +import csv + + +def main(): + chart = go.Figure() + chart.add_trace(go.Scatter(x=[1,2,3], y=[1,2,3])) + chart.show() + + +if __name__ == "__main__": + if len(argv) == 1: + main() + else: + f = open(argv[1], "rt") + #csv_reader = csv.reader(f) + data = {} +# for line in csv_reader: + headers = f.readline() + headers = headers.split(",") + headers[-1] = headers[-1][:-1] + for i in range(len(headers)): + headers[i] = headers[i].strip() + for header in headers: + data[header] = [] + for line in f: + try: + line_data = line.split(",") + #line_data = line + #print(line_data) + for i in range(len(line_data)): + val = float(line_data[i]) + if (i < len(headers)): + header = headers[i] + data[header].append(val) + else: + print("too much values!") + print("headers:", headers) + print("values:", line_data) + except ValueError: + pass + f.close() +# print(data) + chart = go.Figure() + keys = [k for k in data.keys()] +# print("data.keys:", keys) + X_name = keys[0] + print("data samples:",len(data[X_name])) + for key, val in data.items(): + if key != X_name: + chart.add_trace(go.Scatter(x=data[X_name], y=val, name=key, mode="markers+lines")) + chart.show() + diff --git a/C/tmp b/C/tmp new file mode 100644 index 0000000..1c3f27f --- /dev/null +++ b/C/tmp @@ -0,0 +1,1025 @@ +X, Y, FFT +0, 0, 136 +1, 62, 269 +2, 123, 873 +3, 184, 1729 +4, 244, 3568 +5, 303, 6902 +6, 361, 16312 +7, 418, 28544 +8, 473, 73996 +9, 526, 321369 +10, 578, 2295409 +11, 627, 450446 +12, 674, 118120 +13, 718, 55293 +14, 759, 32601 +15, 798, 21779 +16, 834, 16211 +17, 866, 12292 +18, 895, 9817 +19, 921, 7996 +20, 943, 6774 +21, 962, 5788 +22, 977, 10472 +23, 988, 4232 +24, 996, 3741 +25, 1000, 3303 +26, 1000, 3957 +27, 996, 2643 +28, 988, 2440 +29, 977, 2226 +30, 962, 2101 +31, 943, 1834 +32, 921, 1804 +33, 895, 1684 +34, 866, 1547 +35, 834, 1439 +36, 798, 1362 +37, 759, 1258 +38, 718, 1004 +39, 674, 1088 +40, 627, 1039 +41, 578, 966 +42, 526, 2241 +43, 473, 903 +44, 418, 849 +45, 361, 803 +46, 303, 769 +47, 244, 713 +48, 184, 694 +49, 123, 644 +50, 62, 600 +51, 0, 595 +52, -62, 615 +53, -123, 644 +54, -184, 14886 +55, -244, 463 +56, -303, 444 +57, -361, 446 +58, -418, 516 +59, -473, 410 +60, -526, 400 +61, -578, 386 +62, -627, 370 +63, -674, 352 +64, -718, 372 +65, -759, 348 +66, -798, 342 +67, -834, 334 +68, -866, 326 +69, -895, 310 +70, -921, 274 +71, -943, 293 +72, -962, 285 +73, -977, 305 +74, -988, 1976 +75, -996, 230 +76, -1000, 236 +77, -1000, 236 +78, -996, 221 +79, -988, 215 +80, -977, 218 +81, -962, 216 +82, -943, 207 +83, -921, 192 +84, -895, 194 +85, -866, 194 +86, -834, 883 +87, -798, 176 +88, -759, 168 +89, -718, 163 +90, -674, 169 +91, -627, 169 +92, -578, 162 +93, -526, 150 +94, -473, 144 +95, -418, 150 +96, -361, 149 +97, -303, 142 +98, -244, 135 +99, -184, 133 +100, -123, 128 +101, -62, 128 +102, 0, 99 +103, 62, 118 +104, 123, 113 +105, 184, 105 +106, 244, 232 +107, 303, 106 +108, 361, 101 +109, 418, 96 +110, 473, 82 +111, 526, 87 +112, 578, 82 +113, 627, 80 +114, 674, 65 +115, 718, 56 +116, 759, 44 +117, 798, 77 +118, 834, 102926 +119, 866, 399 +120, 895, 207 +121, 921, 156 +122, 943, 131 +123, 962, 115 +124, 977, 106 +125, 988, 96 +126, 996, 88 +127, 1000, 78 +128, 1000, 92 +129, 996, 95 +130, 988, 84 +131, 977, 93 +132, 962, 88 +133, 943, 90 +134, 921, 68 +135, 895, 84 +136, 866, 86 +137, 834, 84 +138, 798, 230 +139, 759, 76 +140, 718, 76 +141, 674, 76 +142, 627, 72 +143, 578, 78 +144, 526, 74 +145, 473, 73 +146, 418, 75 +147, 361, 67 +148, 303, 73 +149, 244, 72 +150, 184, 355 +151, 123, 68 +152, 62, 65 +153, 0, 70 +154, -62, 77 +155, -123, 64 +156, -184, 61 +157, -244, 63 +158, -303, 59 +159, -361, 59 +160, -418, 63 +161, -473, 60 +162, -526, 61 +163, -578, 60 +164, -627, 61 +165, -674, 60 +166, -718, 47 +167, -759, 58 +168, -798, 56 +169, -834, 54 +170, -866, 105 +171, -895, 57 +172, -921, 55 +173, -943, 57 +174, -962, 55 +175, -977, 52 +176, -988, 53 +177, -996, 53 +178, -1000, 49 +179, -1000, 53 +180, -996, 54 +181, -988, 63 +182, -977, 2016 +183, -962, 45 +184, -943, 45 +185, -921, 45 +186, -895, 52 +187, -866, 48 +188, -834, 46 +189, -798, 50 +190, -759, 46 +191, -718, 43 +192, -674, 47 +193, -627, 50 +194, -578, 46 +195, -526, 45 +196, -473, 45 +197, -418, 44 +198, -361, 38 +199, -303, 43 +200, -244, 45 +201, -184, 49 +202, -123, 294 +203, -62, 40 +204, 0, 41 +205, 62, 41 +206, 123, 43 +207, 184, 42 +208, 244, 42 +209, 303, 43 +210, 361, 39 +211, 418, 42 +212, 473, 42 +213, 526, 42 +214, 578, 261 +215, 627, 44 +216, 674, 43 +217, 718, 44 +218, 759, 53 +219, 798, 43 +220, 834, 43 +221, 866, 43 +222, 895, 41 +223, 921, 46 +224, 943, 46 +225, 962, 45 +226, 977, 49 +227, 988, 45 +228, 996, 49 +229, 1000, 51 +230, 1000, 38 +231, 996, 59 +232, 988, 54 +233, 977, 58 +234, 962, 100 +235, 943, 61 +236, 921, 68 +237, 895, 73 +238, 866, 80 +239, 834, 88 +240, 798, 104 +241, 759, 125 +242, 718, 163 +243, 674, 232 +244, 627, 431 +245, 578, 1390 +246, 526, 646388 +247, 473, 688 +248, 418, 139 +249, 361, 48 +250, 303, 7 +251, 244, 21 +252, 184, 14 +253, 123, 13 +254, 62, 12 +255, 0, 16 +256, -62, 26 +257, -123, 26 +258, -184, 25 +259, -244, 26 +260, -303, 26 +261, -361, 26 +262, -418, 23 +263, -473, 28 +264, -526, 31 +265, -578, 45 +266, -627, 3940 +267, -674, 8 +268, -718, 12 +269, -759, 14 +270, -798, 13 +271, -834, 14 +272, -866, 17 +273, -895, 18 +274, -921, 18 +275, -943, 19 +276, -962, 18 +277, -977, 18 +278, -988, 117 +279, -996, 18 +280, -1000, 18 +281, -1000, 18 +282, -996, 23 +283, -988, 17 +284, -977, 18 +285, -962, 16 +286, -943, 17 +287, -921, 18 +288, -895, 18 +289, -866, 19 +290, -834, 18 +291, -798, 14 +292, -759, 19 +293, -718, 18 +294, -674, 13 +295, -627, 18 +296, -578, 18 +297, -526, 18 +298, -473, 34 +299, -418, 19 +300, -361, 19 +301, -303, 19 +302, -244, 17 +303, -184, 17 +304, -123, 18 +305, -62, 18 +306, 0, 17 +307, 62, 18 +308, 123, 19 +309, 184, 19 +310, 244, 760 +311, 303, 16 +312, 361, 15 +313, 418, 16 +314, 473, 17 +315, 526, 15 +316, 578, 16 +317, 627, 17 +318, 674, 17 +319, 718, 16 +320, 759, 17 +321, 798, 17 +322, 834, 16 +323, 866, 16 +324, 895, 16 +325, 921, 16 +326, 943, 14 +327, 962, 15 +328, 977, 15 +329, 988, 16 +330, 996, 98 +331, 1000, 7 +332, 1000, 19 +333, 996, 16 +334, 988, 17 +335, 977, 14 +336, 962, 16 +337, 943, 16 +338, 921, 16 +339, 895, 17 +340, 866, 16 +341, 834, 14 +342, 798, 93 +343, 759, 15 +344, 718, 14 +345, 674, 15 +346, 627, 16 +347, 578, 14 +348, 526, 13 +349, 473, 14 +350, 418, 12 +351, 361, 8 +352, 303, 20 +353, 244, 15 +354, 184, 15 +355, 123, 15 +356, 62, 15 +357, 0, 14 +358, -62, 11 +359, -123, 15 +360, -184, 13 +361, -244, 13 +362, -303, 25 +363, -361, 13 +364, -418, 13 +365, -473, 12 +366, -526, 12 +367, -578, 13 +368, -627, 11 +369, -674, 10 +370, -718, 11 +371, -759, 10 +372, -798, 4 +373, -834, 10 +374, -866, 18555 +375, -895, 70 +376, -921, 34 +377, -943, 27 +378, -962, 22 +379, -977, 21 +380, -988, 19 +381, -996, 20 +382, -1000, 16 +383, -1000, 15 +384, -996, 17 +385, -988, 16 +386, -977, 15 +387, -962, 17 +388, -943, 17 +389, -921, 18 +390, -895, 14 +391, -866, 16 +392, -834, 19 +393, -798, 18 +394, -759, 47 +395, -718, 17 +396, -674, 16 +397, -627, 16 +398, -578, 17 +399, -526, 18 +400, -473, 18 +401, -418, 17 +402, -361, 17 +403, -303, 17 +404, -244, 17 +405, -184, 18 +406, -123, 109 +407, -62, 17 +408, 0, 16 +409, 62, 17 +410, 123, 21 +411, 184, 17 +412, 244, 18 +413, 303, 17 +414, 361, 17 +415, 418, 18 +416, 473, 18 +417, 526, 17 +418, 578, 17 +419, 627, 18 +420, 674, 19 +421, 718, 20 +422, 759, 14 +423, 798, 18 +424, 834, 18 +425, 866, 18 +426, 895, 33 +427, 921, 18 +428, 943, 19 +429, 962, 18 +430, 977, 17 +431, 988, 17 +432, 996, 27 +433, 1000, 24 +434, 1000, 21 +435, 996, 22 +436, 988, 23 +437, 977, 25 +438, 962, 1032 +439, 943, 20 +440, 921, 20 +441, 895, 21 +442, 866, 27 +443, 834, 20 +444, 798, 21 +445, 759, 22 +446, 718, 21 +447, 674, 24 +448, 627, 24 +449, 578, 26 +450, 526, 24 +451, 473, 27 +452, 418, 32 +453, 361, 26 +454, 303, 23 +455, 244, 24 +456, 184, 27 +457, 123, 31 +458, 62, 164 +459, 0, 26 +460, -62, 28 +461, -123, 28 +462, -184, 29 +463, -244, 31 +464, -303, 33 +465, -361, 34 +466, -418, 35 +467, -473, 40 +468, -526, 38 +469, -578, 37 +470, -627, 329 +471, -674, 43 +472, -718, 40 +473, -759, 43 +474, -798, 58 +475, -834, 47 +476, -866, 50 +477, -895, 51 +478, -921, 58 +479, -943, 59 +480, -962, 65 +481, -977, 66 +482, -988, 74 +483, -996, 76 +484, -1000, 82 +485, -1000, 91 +486, -996, 80 +487, -988, 100 +488, -977, 117 +489, -962, 130 +490, -943, 186 +491, -921, 158 +492, -895, 180 +493, -866, 226 +494, -834, 254 +495, -798, 324 +496, -759, 410 +497, -718, 541 +498, -674, 796 +499, -627, 1308 +500, -578, 2740 +501, -526, 10251 +502, -473, 1262107 +503, -418, 6923 +504, -361, 1549 +505, -303, 571 +506, -244, 220 +507, -184, 132 +508, -123, 64 +509, -62, 28 +510, 0, 10 +511, 62, 4 +512, 123, 15 +513, 184, 15 +514, 244, 14 +515, 303, 15 +516, 361, 15 +517, 418, 14 +518, 473, 14 +519, 526, 19 +520, 578, 23 +521, 627, 40 +522, 674, 9215 +523, 718, 9 +524, 759, 6 +525, 798, 7 +526, 834, 8 +527, 866, 9 +528, 895, 9 +529, 921, 8 +530, 943, 9 +531, 962, 9 +532, 977, 8 +533, 988, 10 +534, 996, 69 +535, 1000, 9 +536, 1000, 9 +537, 996, 8 +538, 988, 11 +539, 977, 10 +540, 962, 10 +541, 943, 8 +542, 921, 9 +543, 895, 10 +544, 866, 10 +545, 834, 9 +546, 798, 9 +547, 759, 10 +548, 718, 10 +549, 674, 8 +550, 627, 9 +551, 578, 9 +552, 526, 9 +553, 473, 10 +554, 418, 19 +555, 361, 10 +556, 303, 11 +557, 244, 11 +558, 184, 11 +559, 123, 10 +560, 62, 10 +561, 0, 9 +562, -62, 10 +563, -123, 10 +564, -184, 10 +565, -244, 12 +566, -303, 473 +567, -361, 10 +568, -418, 9 +569, -473, 8 +570, -526, 12 +571, -578, 9 +572, -627, 12 +573, -674, 10 +574, -718, 9 +575, -759, 9 +576, -798, 10 +577, -834, 11 +578, -866, 9 +579, -895, 11 +580, -921, 11 +581, -943, 11 +582, -962, 8 +583, -977, 10 +584, -988, 10 +585, -996, 11 +586, -1000, 66 +587, -1000, 11 +588, -996, 10 +589, -988, 9 +590, -977, 11 +591, -962, 10 +592, -943, 15 +593, -921, 10 +594, -895, 8 +595, -866, 10 +596, -834, 10 +597, -798, 11 +598, -759, 65 +599, -718, 9 +600, -674, 10 +601, -627, 11 +602, -578, 11 +603, -526, 9 +604, -473, 10 +605, -418, 11 +606, -361, 9 +607, -303, 11 +608, -244, 10 +609, -184, 11 +610, -123, 10 +611, -62, 9 +612, 0, 10 +613, 62, 9 +614, 123, 8 +615, 184, 9 +616, 244, 9 +617, 303, 11 +618, 361, 19 +619, 418, 9 +620, 473, 9 +621, 526, 10 +622, 578, 8 +623, 627, 8 +624, 674, 8 +625, 718, 7 +626, 759, 7 +627, 798, 4 +628, 834, 4 +629, 866, 10 +630, 895, 14610 +631, 921, 55 +632, 943, 29 +633, 962, 20 +634, 977, 17 +635, 988, 14 +636, 996, 15 +637, 1000, 14 +638, 1000, 13 +639, 996, 11 +640, 988, 13 +641, 977, 13 +642, 962, 13 +643, 943, 13 +644, 921, 13 +645, 895, 13 +646, 866, 12 +647, 834, 12 +648, 798, 12 +649, 759, 14 +650, 718, 34 +651, 674, 9 +652, 627, 9 +653, 578, 21 +654, 526, 16 +655, 473, 14 +656, 418, 14 +657, 361, 14 +658, 303, 14 +659, 244, 14 +660, 184, 14 +661, 123, 15 +662, 62, 84 +663, 0, 15 +664, -62, 14 +665, -123, 15 +666, -184, 18 +667, -244, 14 +668, -303, 15 +669, -361, 16 +670, -418, 14 +671, -473, 15 +672, -526, 20 +673, -578, 7 +674, -627, 12 +675, -674, 14 +676, -718, 15 +677, -759, 14 +678, -798, 10 +679, -834, 13 +680, -866, 15 +681, -895, 15 +682, -921, 29 +683, -943, 15 +684, -962, 16 +685, -977, 16 +686, -988, 17 +687, -996, 16 +688, -1000, 17 +689, -1000, 16 +690, -996, 17 +691, -988, 18 +692, -977, 23 +693, -962, 12 +694, -943, 792 +695, -921, 15 +696, -895, 15 +697, -866, 15 +698, -834, 18 +699, -798, 15 +700, -759, 16 +701, -718, 17 +702, -674, 16 +703, -627, 17 +704, -578, 18 +705, -526, 17 +706, -473, 19 +707, -418, 18 +708, -361, 18 +709, -303, 19 +710, -244, 16 +711, -184, 17 +712, -123, 19 +713, -62, 23 +714, 0, 119 +715, 62, 18 +716, 123, 19 +717, 184, 19 +718, 244, 17 +719, 303, 21 +720, 361, 20 +721, 418, 20 +722, 473, 20 +723, 526, 22 +724, 578, 21 +725, 627, 23 +726, 674, 156 +727, 718, 20 +728, 759, 22 +729, 798, 23 +730, 834, 27 +731, 866, 23 +732, 895, 25 +733, 921, 18 +734, 943, 24 +735, 962, 27 +736, 977, 27 +737, 988, 27 +738, 996, 27 +739, 1000, 30 +740, 1000, 29 +741, 996, 31 +742, 988, 29 +743, 977, 31 +744, 962, 36 +745, 943, 38 +746, 921, 67 +747, 895, 42 +748, 866, 46 +749, 834, 52 +750, 798, 55 +751, 759, 64 +752, 718, 78 +753, 674, 98 +754, 627, 119 +755, 578, 189 +756, 526, 351 +757, 473, 1191 +758, 418, 569010 +759, 361, 631 +760, 303, 134 +761, 244, 47 +762, 184, 10 +763, 123, 17 +764, 62, 13 +765, 0, 11 +766, -62, 11 +767, -123, 13 +768, -184, 24 +769, -244, 22 +770, -303, 21 +771, -361, 24 +772, -418, 25 +773, -473, 27 +774, -526, 19 +775, -578, 26 +776, -627, 32 +777, -674, 48 +778, -718, 3595 +779, -759, 9 +780, -798, 13 +781, -834, 16 +782, -866, 16 +783, -895, 19 +784, -921, 19 +785, -943, 19 +786, -962, 19 +787, -977, 21 +788, -988, 22 +789, -996, 23 +790, -1000, 144 +791, -1000, 23 +792, -996, 22 +793, -988, 34 +794, -977, 21 +795, -962, 22 +796, -943, 22 +797, -921, 20 +798, -895, 23 +799, -866, 22 +800, -834, 24 +801, -798, 26 +802, -759, 24 +803, -718, 26 +804, -674, 25 +805, -627, 23 +806, -578, 21 +807, -526, 24 +808, -473, 25 +809, -418, 29 +810, -361, 40 +811, -303, 24 +812, -244, 26 +813, -184, 23 +814, -123, 26 +815, -62, 27 +816, 0, 28 +817, 62, 25 +818, 123, 31 +819, 184, 29 +820, 244, 31 +821, 303, 32 +822, 361, 1347 +823, 418, 26 +824, 473, 26 +825, 526, 25 +826, 578, 31 +827, 627, 26 +828, 674, 27 +829, 718, 25 +830, 759, 27 +831, 798, 29 +832, 834, 31 +833, 866, 29 +834, 895, 32 +835, 921, 32 +836, 943, 31 +837, 962, 34 +838, 977, 24 +839, 988, 33 +840, 996, 33 +841, 1000, 35 +842, 1000, 193 +843, 996, 30 +844, 988, 31 +845, 977, 34 +846, 962, 31 +847, 943, 33 +848, 921, 33 +849, 895, 34 +850, 866, 34 +851, 834, 35 +852, 798, 36 +853, 759, 35 +854, 718, 230 +855, 674, 32 +856, 627, 34 +857, 578, 35 +858, 526, 42 +859, 473, 36 +860, 418, 35 +861, 361, 36 +862, 303, 38 +863, 244, 38 +864, 184, 38 +865, 123, 39 +866, 62, 38 +867, 0, 36 +868, -62, 39 +869, -123, 38 +870, -184, 31 +871, -244, 40 +872, -303, 38 +873, -361, 38 +874, -418, 68 +875, -473, 41 +876, -526, 39 +877, -578, 35 +878, -627, 39 +879, -674, 35 +880, -718, 34 +881, -759, 31 +882, -798, 28 +883, -834, 23 +884, -866, 21 +885, -895, 49 +886, -921, 73240 +887, -943, 262 +888, -962, 134 +889, -977, 103 +890, -988, 85 +891, -996, 82 +892, -1000, 77 +893, -1000, 75 +894, -996, 60 +895, -988, 64 +896, -977, 67 +897, -962, 67 +898, -943, 67 +899, -921, 72 +900, -895, 71 +901, -866, 75 +902, -834, 58 +903, -798, 74 +904, -759, 78 +905, -718, 79 +906, -674, 184 +907, -627, 79 +908, -578, 81 +909, -526, 79 +910, -473, 81 +911, -418, 88 +912, -361, 90 +913, -303, 90 +914, -244, 80 +915, -184, 90 +916, -123, 91 +917, -62, 96 +918, 0, 649 +919, 62, 94 +920, 123, 95 +921, 184, 98 +922, 244, 111 +923, 303, 102 +924, 361, 104 +925, 418, 106 +926, 473, 103 +927, 526, 110 +928, 578, 116 +929, 627, 122 +930, 674, 119 +931, 718, 122 +932, 759, 128 +933, 798, 133 +934, 834, 85 +935, 866, 120 +936, 895, 129 +937, 921, 130 +938, 943, 217 +939, 962, 143 +940, 977, 149 +941, 988, 148 +942, 996, 155 +943, 1000, 161 +944, 1000, 166 +945, 996, 164 +946, 988, 171 +947, 977, 183 +948, 962, 195 +949, 943, 231 +950, 921, 10578 +951, 895, 173 +952, 866, 176 +953, 834, 184 +954, 798, 234 +955, 759, 190 +956, 718, 214 +957, 674, 211 +958, 627, 222 +959, 578, 226 +960, 526, 258 +961, 473, 266 +962, 418, 270 +963, 361, 281 +964, 303, 302 +965, 244, 309 +966, 184, 253 +967, 123, 315 +968, 62, 333 +969, 0, 362 +970, -62, 1713 +971, -123, 328 +972, -184, 369 +973, -244, 380 +974, -303, 376 +975, -361, 418 +976, -418, 447 +977, -473, 461 +978, -526, 490 +979, -578, 523 +980, -627, 544 +981, -674, 588 +982, -718, 5360 +983, -759, 626 +984, -798, 652 +985, -834, 692 +986, -866, 813 +987, -895, 780 +988, -921, 838 +989, -943, 879 +990, -962, 943 +991, -977, 1022 +992, -988, 1131 +993, -996, 1185 +994, -1000, 1370 +995, -1000, 1428 +996, -996, 1569 +997, -988, 1718 +998, -977, 1486 +999, -962, 2039 +1000, -943, 2341 +1001, -921, 2598 +1002, -895, 3162 +1003, -866, 3530 +1004, -834, 4107 +1005, -798, 4856 +1006, -759, 5951 +1007, -718, 7453 +1008, -674, 9738 +1009, -627, 13148 +1010, -578, 19483 +1011, -526, 32817 +1012, -473, 69678 +1013, -418, 264596 +1014, -361, 354996 +1015, -303, 185507 +1016, -244, 42489 +1017, -184, 16180 +1018, -123, 7601 +1019, -62, 3926 +1020, 0, 1991 +1021, 62, 933 +1022, 123, 435 +1023, 184, 117 diff --git a/FFT.py b/FFT.py index c2cef32..7efa27d 100644 --- a/FFT.py +++ b/FFT.py @@ -179,7 +179,7 @@ def FFT_arr_inplace(buf): -FP_acc = 1e3 +FP_acc = 2<<16 def FFT_arr_FP(buf): """ @@ -202,12 +202,19 @@ def FFT_arr_FP(buf): buf[i*2], buf[i*2+1], buf[j*2], buf[j*2+1] = buf[j*2], buf[j*2+1], buf[i*2], buf[i*2+1] # --- уровни бабочек --- + + c_delta_FP_arr = [cos(2.0 * pi / m) * FP_acc for m in range(2,N+2)] + s_delta_FP_arr = [-sin(2.0 * pi / m) * FP_acc for m in range(2,N+2)] m = 2 while m <= N: half = m // 2 # шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw) - c_delta_FP = cos(2.0 * pi / m) * FP_acc - s_delta_FP = -sin(2.0 * pi / m) * FP_acc + print("N, m:", N,m) + + c_delta_FP = c_delta_FP_arr[m -2] + s_delta_FP = s_delta_FP_arr[m -2] + #c_delta_FP = cos(2.0 * pi / m) * FP_acc + #s_delta_FP = -sin(2.0 * pi / m) * FP_acc for start in range(0, N, m): # w = e^{-j*0*Δ} = 1 + j0 @@ -238,7 +245,74 @@ def FFT_arr_FP(buf): buf = [val/ FP_acc for val in buf] return buf - +#global arrays. calculated once at compile time. They store sin and cos values as fixed-point + c_delta_FP_arr = [int(cos(2.0 * pi / m) * FP_acc) for m in range(2,inp_L//2 + 2)] + s_delta_FP_arr = [int(-sin(2.0 * pi / m) * FP_acc) for m in range(2, inp_L//2 + 2)] + +def FFT_arr_FP_for_C(inp, inp_L, buf): #version for translation to C directly + """ + In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im]. + Без массивов twiddle: твиддл на уровне обновляется рекуррентно. + """ + + #buf = [int(val*FP_acc) for val in buf] + for i in range(inp_L): + buf[i*2] = inp[i] #take data from input as real and store it to the buf as Re, Im pairs. + buf[i*2 + 1] = 0 + + N = inp_L//2 + # --- bit-reverse перестановка (чтобы бабочки шли последовательно) --- + j = 0 + for i in range(1, N): + bit = N >> 1 + while j & bit: + j ^= bit + bit >>= 1 + j |= bit + if i < j: + buf[i*2], buf[i*2+1], buf[j*2], buf[j*2+1] = buf[j*2], buf[j*2+1], buf[i*2], buf[i*2+1] + + # --- уровни бабочек --- + + c_delta_FP_arr = [cos(2.0 * pi / m) * FP_acc for m in range(2,N+2)] + s_delta_FP_arr = [-sin(2.0 * pi / m) * FP_acc for m in range(2,N+2)] + m = 2 + while m <= N: + half = m // 2 + # шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw) + #print("N, m:", N,m) + + c_delta_FP = c_delta_FP_arr[m -2] + s_delta_FP = s_delta_FP_arr[m -2] + #c_delta_FP = cos(2.0 * pi / m) * FP_acc + #s_delta_FP = -sin(2.0 * pi / m) * FP_acc + + for start in range(0, N, m): + # w = e^{-j*0*Δ} = 1 + j0 + wr, wi = 1.0 * FP_acc, 0.0 * FP_acc + for k in range(half): + u_re, u_im = buf[(start + k)*2], buf[(start + k)*2 +1] + v_re, v_im = buf[(start + k + half)*2], buf[(start + k + half)*2 +1] + + # t = w * v + #print("wr, v_re, wi, v_im:", wr, v_re, wi, v_im) + t_re = (wr * v_re - wi * v_im) / FP_acc + t_im = (wr * v_im + wi * v_re) / FP_acc + + # верх/низ + buf[(start + k)*2 +0] = u_re + t_re + buf[(start + k)*2 +1] = u_im + t_im + buf[(start + k + half)*2 + 0] = u_re - t_re + buf[(start + k + half)*2 + 1] = u_im - t_im + + + # w *= w_m (поворот на Δ с помощью рекуррентной формулы) + # (wr + j wi) * (cΔ + j sΔ) + wr, wi = (wr * c_delta_FP - wi * s_delta_FP)/ FP_acc, (wr * s_delta_FP + wi * c_delta_FP)/ FP_acc + + + m <<= 1 + diff --git a/__pycache__/FFT.cpython-312.pyc b/__pycache__/FFT.cpython-312.pyc index 4a1791b..d0fcf14 100644 Binary files a/__pycache__/FFT.cpython-312.pyc and b/__pycache__/FFT.cpython-312.pyc differ