From 1e06bb4dc734b927c9c8d374bcd82f0a00c322fb Mon Sep 17 00:00:00 2001 From: chatlanin Date: Mon, 16 Feb 2026 19:47:03 +0300 Subject: [PATCH] initial commit --- .gitignore | 5 + README.md | 0 bin/main.cpp | 23 + bin/meson.build | 8 + meson.build | 39 ++ run.sh | 29 + src/adapter/adapter.hpp | 75 +++ src/harmonica.hpp | 95 +++ src/meson.build | 46 ++ src/plugins/magnitude/magnitude.cpp | 28 + src/plugins/magnitude/magnitude.hpp | 21 + src/utils/cvec/cvec.cpp | 21 + src/utils/cvec/cvec.hpp | 18 + src/utils/fft/fft.hpp | 116 ++++ src/utils/fft/rdft.hpp | 924 ++++++++++++++++++++++++++++ src/utils/fvec/fvec.cpp | 60 ++ src/utils/fvec/fvec.hpp | 44 ++ src/utils/real_time/real_time.cpp | 18 + src/utils/real_time/real_time.hpp | 48 ++ src/utils/using.hpp | 14 + src/utils/var.hpp | 9 + src/utils/windows/hann/hann.cpp | 13 + src/utils/windows/hann/hann.hpp | 33 + src/utils/workers/plugin.hpp | 28 + src/utils/workers/result.hpp | 39 ++ src/utils/workers/setup.hpp | 25 + subprojects/hack.wrap | 7 + 27 files changed, 1786 insertions(+) create mode 100644 .gitignore create mode 100644 README.md create mode 100644 bin/main.cpp create mode 100644 bin/meson.build create mode 100644 meson.build create mode 100755 run.sh create mode 100644 src/adapter/adapter.hpp create mode 100644 src/harmonica.hpp create mode 100644 src/meson.build create mode 100644 src/plugins/magnitude/magnitude.cpp create mode 100644 src/plugins/magnitude/magnitude.hpp create mode 100644 src/utils/cvec/cvec.cpp create mode 100644 src/utils/cvec/cvec.hpp create mode 100644 src/utils/fft/fft.hpp create mode 100644 src/utils/fft/rdft.hpp create mode 100644 src/utils/fvec/fvec.cpp create mode 100644 src/utils/fvec/fvec.hpp create mode 100644 src/utils/real_time/real_time.cpp create mode 100644 src/utils/real_time/real_time.hpp create mode 100644 src/utils/using.hpp create mode 100644 src/utils/var.hpp create mode 100644 src/utils/windows/hann/hann.cpp create mode 100644 src/utils/windows/hann/hann.hpp create mode 100644 src/utils/workers/plugin.hpp create mode 100644 src/utils/workers/result.hpp create mode 100644 src/utils/workers/setup.hpp create mode 100644 subprojects/hack.wrap diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5b2a04f --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +build +.cache +subprojects/* +!subprojects/*.wrap +src/noincl.hpp diff --git a/README.md b/README.md new file mode 100644 index 0000000..e69de29 diff --git a/bin/main.cpp b/bin/main.cpp new file mode 100644 index 0000000..ed63dbb --- /dev/null +++ b/bin/main.cpp @@ -0,0 +1,23 @@ +#include +#include "harmonica.hpp" + +auto main() -> int +{ + // setup создается для каждого файла свой + // т.к. при чтении данных из файла уже должен быть определен размер блока + // данных для чтения m_block_size; см. установки по умолчанию. + // Передается по ссылке и заполняется необходимыми данными + hr::setup setup; + setup.m_domain = hr::DOMAIN_PLUGIN::FREQUENSY; + setup.m_file = hr::var::SOUND; + + auto r = hr::run(setup); + hack::log()("size:", r.m_data.size()); + + if (!r.empty()) + { + std::vector s; + for (auto p : r.m_data) s.push_back(p.m_value[0]); + // hack::log()(s); + } +} diff --git a/bin/meson.build b/bin/meson.build new file mode 100644 index 0000000..f6ed8ee --- /dev/null +++ b/bin/meson.build @@ -0,0 +1,8 @@ +executable( + meson.project_name(), + 'main.cpp', + dependencies : deps, + cpp_args: args, + include_directories : inc +) + diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..8c9f340 --- /dev/null +++ b/meson.build @@ -0,0 +1,39 @@ +project( + meson.current_source_dir().split('/').get(-1), + 'cpp', + version : run_command('git', 'rev-parse', '--short', 'HEAD', check: false).stdout().strip(), + default_options : [ + 'warning_level=1', + 'optimization=3', + 'cpp_std=c++20', +]) + +add_project_arguments ( + '-Wpedantic', + '-Wno-shadow', + '-Wno-unused-but-set-variable', + '-Wno-comment', + '-Wno-unused-parameter', + '-Wno-unused-value', + '-Wno-missing-field-initializers', + '-Wno-narrowing', + '-Wno-deprecated-enum-enum-conversion', + '-Wno-volatile', + '-Wno-format-security', + '-Wno-switch', + '-Wno-ignored-attributes', + '-Wno-unused-variable', + '-Wno-deprecated-enum-enum-conversion', + language: 'cpp' +) + +args = [] +inc = [] +deps = [ + dependency('sndfile'), + dependency('cblas'), + subproject('hack').get_variable('hack_dep') +] + +subdir('src') +subdir('bin') diff --git a/run.sh b/run.sh new file mode 100755 index 0000000..dd4589e --- /dev/null +++ b/run.sh @@ -0,0 +1,29 @@ +#!/bin/zsh + +PROJECT_NAME=$(basename $PWD) + +run() { + command meson compile -C build + cd build + ./bin/$PROJECT_NAME + cd .. +} + +# run test [name_test] +# example: run test pattrens +if [[ "$1" == "test" ]]; then + echo "" + meson test $2 -C build + echo "" + awk '/^-------------------------------------------------------------------------------/{flag=1} /===============================================================================/{flag=0} flag' ./build/meson-logs/testlog.txt +elif [[ "$1" == "tests" ]]; then + echo "" + meson test -C build + echo "" + # awk '/^-------------------------------------------------------------------------------/{flag=1} /===============================================================================/{flag=0} flag' ./build/meson-logs/testlog.txt +elif [[ -d "build" ]]; then + run +else + command meson setup build + run +fi diff --git a/src/adapter/adapter.hpp b/src/adapter/adapter.hpp new file mode 100644 index 0000000..a113be6 --- /dev/null +++ b/src/adapter/adapter.hpp @@ -0,0 +1,75 @@ +#pragma once + +#include "utils/workers/result.hpp" +#include "utils/workers/setup.hpp" +#include "utils/fft/fft.hpp" +#include "utils/fvec/fvec.hpp" +#include "utils/windows/hann/hann.hpp" + +namespace hr +{ + template + class adapter + { + public: + adapter(Plugin& p) : m_plugin { p } + { + m_end = m_plugin.m_setup.m_block_size - m_plugin.m_setup.m_step_size; + m_data.resize(m_plugin.m_setup.m_block_size, 0.0); + m_data_old.resize(m_end, 0.0); + m_fft.init(m_plugin.m_setup.m_block_size); + m_hann.creaate(m_plugin.m_setup.m_block_size); + } + + virtual ~adapter() { } + + protected: + Plugin& m_plugin; + real_time m_timestamp; + fft m_fft; + math::hann m_hann; + fvec_t m_data; + fvec_t m_data_old; + size_t m_end; + + public: + void process(fvec_t& in, real_time timestamp) + { + m_timestamp = timestamp; + + switch(m_plugin.m_setup.m_domain) + { + case DOMAIN_PLUGIN::TIME: + { + m_plugin.process(in, m_timestamp); + break; + } + case DOMAIN_PLUGIN::FREQUENSY: + { + frequensy(in); + break; + } + } + } + + result get_result() { return m_plugin.get_result(); } + + private: + void swap_buffer(fvec_t& in) + { + size_t i = 0; + for (i = 0; i < m_end; ++i) m_data[i] = m_data_old[i]; + for (i = 0; i < m_plugin.m_setup.m_step_size; ++i) m_data[m_end + i] = in[i]; + for (i = 0; i < m_end; ++i) m_data_old[i] = m_data[i + m_plugin.m_setup.m_step_size]; + } + + void frequensy(fvec_t& in) + { + swap_buffer(in); + m_hann.apply(m_data); + m_data.shift(); + auto r = m_fft.process(m_data); + m_plugin.process(r, in, m_timestamp); + } + }; +} diff --git a/src/harmonica.hpp b/src/harmonica.hpp new file mode 100644 index 0000000..022b287 --- /dev/null +++ b/src/harmonica.hpp @@ -0,0 +1,95 @@ +#pragma once + +#include +#include + +#include "utils/var.hpp" // IWYU pragma: keep +#include "utils/using.hpp" +#include "utils/workers/setup.hpp" +#include "utils/workers/result.hpp" +#include "adapter/adapter.hpp" + +#include "plugins/magnitude/magnitude.hpp" // IWYU pragma: keep + +namespace hr +{ + /** + * @brief Запуск на чтение аудиофайла и обработка его через плагин + * @tparam Plugin Тип плагина для обработки аудиоданных + * @param setup Настройки обработки (путь к файлу, параметры и т.д.) + * @return Результат обработки аудио + */ + template + inline result run(setup& setup) + { + // Инициализация структуры для libsndfile и открытие файла + SF_INFO sf_info; + SNDFILE* file = sf_open(setup.m_file.c_str(), SFM_READ, &sf_info); + if (!file) + { + // Обработка ошибки открытия файла + hack::exception ex; + hack::log().on_file(); + hack::log().on_func(); + hack::log().on_row(); + ex.title("Error of open file"); + ex.description(sf_strerror(file)); + ex.set("file", setup.m_file); + hack::error()(ex); + throw ex; + } + + // Сохранение информации о файле в настройки + setup.m_sample_rate = sf_info.samplerate; + setup.m_frames = sf_info.frames; + setup.m_channels = sf_info.channels; + + if (setup.m_channels == 0) throw std::runtime_error("Нет каналов в аудио файле"); + + // Инициализация переменных для чтения + std::size_t read = 0; // Количество обработанных кадров + fvec_t read_data(setup.m_channels * setup.m_step_size, .0); // Буфер для чтения (интерливированные данные) + fvec_t in(setup.m_step_size, .0); // Буфер для моно-данных + + // Создание плагина и адаптера для обработки + Plugin pl { setup }; + adapter ad { pl }; + + do + { + // Определение длины читаемого блока (защита от выхода за границы) + auto length = hack::math::min(setup.m_step_size, in.size()); + auto read_samples = sf_read_float(file, read_data.data(), read_data.size()); + + uint_t read_length = read_samples / setup.m_channels; // Перевод в кадры + read_length = hack::math::min(length, read_length); // Ограничение длиной буфера + + // Де-интерливирование и down-mixing (преобразование многоканального в моно) + for (std::size_t i = 0; i < read_length; ++i) + { + in[i] = 0.0; + // Суммирование всех каналов + for (int c = 0; c < setup.m_channels; ++c) + in[i] += read_data[setup.m_channels * i + c]; + // Усреднение для получения моно-сигнала + in[i] /= static_cast(setup.m_channels); + } + + // Подготовка к следующей итерации + read = hack::math::min(length, static_cast(floorf(read_length + .5))); + + // Дополнение буфера нулями если считано неполный блок (конец файла) + if (in.size() > read) std::fill(in.begin() + read, in.end(), 0.0); + + // Вычисление временной метки и обработка данных через адаптер + real_time timestamp = real_time::frame2rt(read, sf_info.samplerate); + ad.process(in, timestamp); + } + while (read == setup.m_step_size); // Продолжать пока читаются полные блоки + + // Закрытие файла и возврат результата + sf_close(file); + + return ad.get_result(); + } +} diff --git a/src/meson.build b/src/meson.build new file mode 100644 index 0000000..80bcacb --- /dev/null +++ b/src/meson.build @@ -0,0 +1,46 @@ +inc += include_directories('.') + +headers = [ + 'adapter/adapter.hpp', + + 'utils/fft/fft.hpp', + 'utils/fft/rdft.hpp', + 'utils/cvec/cvec.hpp', + 'utils/fvec/fvec.hpp', + 'utils/real_time/real_time.hpp', + 'utils/workers/plugin.hpp', + 'utils/workers/result.hpp', + 'utils/workers/setup.hpp', + 'utils/windows/hann/hann.hpp', + + # plugins + 'plugins/magnitude/magnitude.hpp', + + 'harmonica.hpp' +] + +sources = [ + 'utils/real_time/real_time.cpp', + 'utils/cvec/cvec.cpp', + 'utils/fvec/fvec.cpp', + + 'utils/windows/hann/hann.cpp', + + # plugins + 'plugins/magnitude/magnitude.cpp', +] + +lib = library( + meson.project_name(), + include_directories : inc, + sources: [headers, sources], + dependencies : deps, + cpp_args: args +) + +harmonica_sdk_dep = declare_dependency( + include_directories: inc, + link_with: lib +) + +deps += harmonica_sdk_dep diff --git a/src/plugins/magnitude/magnitude.cpp b/src/plugins/magnitude/magnitude.cpp new file mode 100644 index 0000000..5e86ea8 --- /dev/null +++ b/src/plugins/magnitude/magnitude.cpp @@ -0,0 +1,28 @@ +#include "magnitude.hpp" + +namespace hr::plugins +{ + magnitude::magnitude(const setup& st) : plugin{ st } + { + } + + void magnitude::process(cvec_t& fft, fvec_t& base, real_time timestamp) + { + std::size_t frames = fft.size(); + + base_t m = 0.f; + for (std::size_t i = 0; i < frames; ++i) m += fft.m_norm[i]; + m /= frames; + + result::bit b; + b.m_value.push_back(m); + b.m_duration = timestamp; + + m_result.set_bit(b); + } + + result magnitude::get_result() + { + return m_result; + } +} diff --git a/src/plugins/magnitude/magnitude.hpp b/src/plugins/magnitude/magnitude.hpp new file mode 100644 index 0000000..56ba2f6 --- /dev/null +++ b/src/plugins/magnitude/magnitude.hpp @@ -0,0 +1,21 @@ +#pragma once + +#include "utils/workers/plugin.hpp" + +namespace hr::plugins +{ + class magnitude : public plugin + { + public: + magnitude(const setup& st); + virtual ~magnitude() = default; + + private: + result m_result; + + public: + void process(fvec_t& base, real_time timestamp) override {}; + void process(cvec_t& fft, fvec_t& base, real_time timestamp) override; + result get_result() override; + }; +} diff --git a/src/utils/cvec/cvec.cpp b/src/utils/cvec/cvec.cpp new file mode 100644 index 0000000..a0dfb86 --- /dev/null +++ b/src/utils/cvec/cvec.cpp @@ -0,0 +1,21 @@ +#include "cvec.hpp" +#include + +namespace hr +{ + cvec_t::cvec_t(std::size_t s, base_t v) + { + resize(s, v); + } + + size_t cvec_t::size() const noexcept + { + return m_norm.size(); + } + + void cvec_t::resize(std::size_t s, base_t v) + { + m_norm.resize(s, v); + m_phas.resize(s, v); + } +} diff --git a/src/utils/cvec/cvec.hpp b/src/utils/cvec/cvec.hpp new file mode 100644 index 0000000..560329a --- /dev/null +++ b/src/utils/cvec/cvec.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include "utils/fvec/fvec.hpp" + +namespace hr +{ + struct cvec_t + { + cvec_t() =default; + cvec_t(std::size_t s, base_t v = 0.0); + fvec_t m_norm; + fvec_t m_phas; + + size_t size() const noexcept; + void resize(std::size_t s, base_t v = 0.0); + }; +} + diff --git a/src/utils/fft/fft.hpp b/src/utils/fft/fft.hpp new file mode 100644 index 0000000..139b854 --- /dev/null +++ b/src/utils/fft/fft.hpp @@ -0,0 +1,116 @@ +#pragma once + +#include "utils/using.hpp" +#include "utils/cvec/cvec.hpp" +#include "rdft.hpp" + +#include + +namespace hr +{ + struct fft + { + uint_t m_win_size; + uint_t m_fft_size; + fvec_t m_w; + ivec_t m_ip; + fvec_t m_complex_spectr; + + void init(size_t size) noexcept + { + m_win_size = size; + m_fft_size = m_win_size / 2 + 1; + m_complex_spectr.resize(m_win_size, 0.0); + m_ip.resize(m_fft_size, 0); + m_w.resize(m_fft_size, 0.0); + } + + cvec_t process(fvec_t& in) + { + m_complex_spectr = get_complex_spectr(in); + + // Подготавливаем результат: векторы для фазы и амплитуды + cvec_t res; + res.m_phas.resize(m_fft_size, 0.0); // Фаза в радианах для каждой частотной компоненты + res.m_norm.resize(m_fft_size, 0.0); // Амплитуда (магнитуда) для каждой частотной компоненты + + // =========================================================================== + // ВЫЧИСЛЕНИЕ ФАЗЫ: + // Для постоянной составляющей (0 Гц): фаза может быть 0 или π в зависимости от знака + if (m_complex_spectr[0] < 0.0) + res.m_phas[0] = std::numbers::pi; // Отрицательное значение → фаза 180° + else + res.m_phas[0] = 0.0; // Положительное значение → фаза 0° + + // Для основных частот (от 1 до N/2-1): вычисляем фазу через арктангенс + for (size_t i = 1; i < m_fft_size - 1; ++i) + { + // atan2(Im, Re) дает правильную фазу с учетом квадрантов + // m_complex_spectr[m_win_size - i] - мнимая часть + // m_complex_spectr[i] - действительная часть + res.m_phas[i] = std::atan2f(m_complex_spectr[m_win_size - i], m_complex_spectr[i]); + } + + // Для частоты Найквиста (N/2): аналогично постоянной составляющей + if (m_complex_spectr[m_fft_size / 2] < 0) + // res.m_phas[m_fft_size - 1] = std::numbers::pi; // Отрицательное значение → фаза 180° + res.m_phas[m_fft_size - 1] = 3.14159265358979323846; + else + res.m_phas[m_fft_size - 1] = 0.0; // Положительное значение → фаза 0° + + // =========================================================================== + // ВЫЧИСЛЕНИЕ АМПЛИТУДЫ (МАГНИТУДЫ): + // Для постоянной составляющей (0 Гц): просто абсолютное значение + res.m_norm[0] = std::fabsf(m_complex_spectr[0]); + + // Для основных частот (от 1 до N/2-1): вычисляем магнитуду как sqrt(Re² + Im²) + for (size_t i = 1; i < m_fft_size - 1; ++i) + // Теорема Пифагора для комплексных чисел + // закоментированные ниже не дают точность по сравнению с аубио + // res.m_norm[i] = std::hypotf(m_complex_spectr[i], m_complex_spectr[m_complex_spectr.size() - i]); + // res.m_norm[i] = std::sqrtf(std::pow(m_complex_spectr[i], 2)+ std::pow(m_complex_spectr[m_complex_spectr.size() - i], 2)); + res.m_norm[i] = std::sqrtf(m_complex_spectr[i] * m_complex_spectr[i] + m_complex_spectr[m_complex_spectr.size() - i] * m_complex_spectr[m_complex_spectr.size() - i]); + + + // Для частоты Найквиста (N/2): просто абсолютное значение + res.m_norm[res.m_norm.size() - 1] = std::fabs(m_complex_spectr[m_complex_spectr.size() / 2]); + + return res; + } + + fvec_t get_complex_spectr(fvec_t& in) + { + fvec_t complex_spectr(m_win_size, 0.0); + + // Выполняем прямое вещественное БПФ + // Результат сохраняется в том же массиве 'in' в специальном формате + rdft(m_win_size, in.data(), m_ip.data(), m_w.data()); + + // ПРЕОБРАЗОВАНИЕ: конвертируем упакованный формат БПФ в полноценный комплексный спектр + // Формат после rdft(): + // in[0] = Re[0] // Постоянная составляющая (0 Гц) + // in[1] = Re[N/2] // Частота Найквиста (f_sample/2) + // in[2] = Re[1] // Действительная часть 1-й частоты + // in[3] = Im[1] // Мнимая часть 1-й частоты + // in[4] = Re[2] // Действительная часть 2-й частоты + // in[5] = Im[2] // Мнимая часть 2-й частоты + // in[6] = Re[3] // Действительная часть 3-й частоты + // in[7] = Im[3] // Мнимая часть 3-й частоты + // ... + // in[N-2] = Re[N/2-1] // Действительная часть (N/2-1)-й частоты + // in[N-1] = Im[N/2-1] // Мнимая часть (N/2-1)-й частоты + + complex_spectr[0] = in[0]; // Постоянная составляющая (0 Гц) - только действительная часть + complex_spectr[m_win_size / 2] = in[1]; // Частота Найквиста (f_sample/2) - только действительная часть + + // Заполняем комплексный спектр для частот от 1 до N/2-1 + for (size_t i = 1; i < m_fft_size - 1; ++i) + { + complex_spectr[i] = in[2 * i]; // Действительная часть i-й частоты + complex_spectr[m_win_size - i] = -in[2 * i + 1]; // Мнимая часть i-й частоты (знак изменен для правильной фазы) + } + + return complex_spectr; + } + }; +} diff --git a/src/utils/fft/rdft.hpp b/src/utils/fft/rdft.hpp new file mode 100644 index 0000000..2124145 --- /dev/null +++ b/src/utils/fft/rdft.hpp @@ -0,0 +1,924 @@ +#include "utils/using.hpp" +#include + +namespace hr +{ + namespace { + void makewt(int nw, int *ip, base_t *w); + void makect(int nc, int *ip, base_t *c); + void bitrv2(int n, int *ip, base_t *a); + void cftfsub(int n, base_t *a, base_t *w); + void cftbsub(int n, base_t *a, base_t *w); + void cft1st(int n, base_t *a, base_t *w); + void cftmdl(int n, int l, base_t *a, base_t *w); + void rftfsub(int n, base_t *a, int nc, base_t *c); + void rftbsub(int n, base_t *a, int nc, base_t *c); + + void makewt(int nw, int* ip, base_t* w) + { + int j, nwh; + base_t delta, x, y; + + ip[0] = nw; + ip[1] = 1; + + if (nw > 2) + { + nwh = nw >> 1; + delta = std::atan(1.0) / nwh; + w[0] = 1; + w[1] = 0; + w[nwh] = std::cos(delta * nwh); + w[nwh + 1] = w[nwh]; + if (nwh > 2) + { + for (j = 2; j < nwh; j += 2) + { + x = std::cos(delta * j); + y = std::sin(delta * j); + w[j] = x; + w[j + 1] = y; + w[nw - j] = y; + w[nw - j + 1] = x; + } + for (j = nwh - 2; j >= 2; j -= 2) + { + x = w[2 * j]; + y = w[2 * j + 1]; + w[nwh + j] = x; + w[nwh + j + 1] = y; + } + + + bitrv2(nw, ip + 2, w); + } + } + } + + void makect(int nc, int* ip, base_t* c) + { + int j, nch; + base_t delta; + + ip[1] = nc; + if (nc > 1) + { + nch = nc >> 1; + delta = std::atan(1.0) / nch; + c[0] = std::cos(delta * nch); + c[nch] = (base_t)0.5 * c[0]; + for (j = 1; j < nch; j++) + { + c[j] = (base_t)0.5 * std::cos(delta * j); + c[nc - j] = (base_t)0.5 * std::sin(delta * j); + } + } + } + + void bitrv2(int n, int* ip, base_t* a) + { + int j, j1, k, k1, l, m, m2; + base_t xr, xi, yr, yi; + + ip[0] = 0; + l = n; + m = 1; + while ((m << 3) < l) + { + l >>= 1; + for (j = 0; j < m; j++) + ip[m + j] = ip[j] + l; + m <<= 1; + } + m2 = 2 * m; + if ((m << 3) == l) + { + for (k = 0; k < m; k++) + { + for (j = 0; j < k; j++) + { + j1 = 2 * j + ip[k]; + k1 = 2 * k + ip[j]; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + j1 += m2; + k1 += 2 * m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + j1 += m2; + k1 -= m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + j1 += m2; + k1 += 2 * m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + } + j1 = 2 * k + m2 + ip[k]; + k1 = j1 + m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + } + } + else + { + for (k = 1; k < m; k++) + { + for (j = 0; j < k; j++) + { + j1 = 2 * j + ip[k]; + k1 = 2 * k + ip[j]; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + j1 += m2; + k1 += m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + } + } + } + } + + void cftfsub(int n, base_t* a, base_t* w) + { + int j, j1, j2, j3, l; + base_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; + + l = 2; + if (n >= 16) + { + cft1st(n, a, w); + l = 16; + while ((l << 3) <= n) + { + cftmdl(n, l, a, w); + l <<= 3; + } + } + if ((l << 1) < n) + { + for (j = 0; j < l; j += 2) + { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + a[j] = x0r + x2r; + a[j + 1] = x0i + x2i; + a[j2] = x0r - x2r; + a[j2 + 1] = x0i - x2i; + a[j1] = x1r - x3i; + a[j1 + 1] = x1i + x3r; + a[j3] = x1r + x3i; + a[j3 + 1] = x1i - x3r; + } + } + else if ((l << 1) == n) + { + for (j = 0; j < l; j += 2) + { + j1 = j + l; + x0r = a[j] - a[j1]; + x0i = a[j + 1] - a[j1 + 1]; + a[j] += a[j1]; + a[j + 1] += a[j1 + 1]; + a[j1] = x0r; + a[j1 + 1] = x0i; + } + } + } + + + void cftbsub(int n, base_t* a, base_t* w) + { + int j, j1, j2, j3, j4, j5, j6, j7, l; + base_t wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i; + + l = 2; + if (n > 16) + { + cft1st(n, a, w); + l = 16; + while ((l << 3) < n) + { + cftmdl(n, l, a, w); + l <<= 3; + } + } + if ((l << 2) < n) + { + wn4r = w[2]; + for (j = 0; j < l; j += 2) + { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + j4 = j3 + l; + j5 = j4 + l; + j6 = j5 + l; + j7 = j6 + l; + x0r = a[j] + a[j1]; + x0i = -a[j + 1] - a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = -a[j + 1] + a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + y0r = x0r + x2r; + y0i = x0i - x2i; + y2r = x0r - x2r; + y2i = x0i + x2i; + y1r = x1r - x3i; + y1i = x1i - x3r; + y3r = x1r + x3i; + y3i = x1i + x3r; + x0r = a[j4] + a[j5]; + x0i = a[j4 + 1] + a[j5 + 1]; + x1r = a[j4] - a[j5]; + x1i = a[j4 + 1] - a[j5 + 1]; + x2r = a[j6] + a[j7]; + x2i = a[j6 + 1] + a[j7 + 1]; + x3r = a[j6] - a[j7]; + x3i = a[j6 + 1] - a[j7 + 1]; + y4r = x0r + x2r; + y4i = x0i + x2i; + y6r = x0r - x2r; + y6i = x0i - x2i; + x0r = x1r - x3i; + x0i = x1i + x3r; + x2r = x1r + x3i; + x2i = x1i - x3r; + y5r = wn4r * (x0r - x0i); + y5i = wn4r * (x0r + x0i); + y7r = wn4r * (x2r - x2i); + y7i = wn4r * (x2r + x2i); + a[j1] = y1r + y5r; + a[j1 + 1] = y1i - y5i; + a[j5] = y1r - y5r; + a[j5 + 1] = y1i + y5i; + a[j3] = y3r - y7i; + a[j3 + 1] = y3i - y7r; + a[j7] = y3r + y7i; + a[j7 + 1] = y3i + y7r; + a[j] = y0r + y4r; + a[j + 1] = y0i - y4i; + a[j4] = y0r - y4r; + a[j4 + 1] = y0i + y4i; + a[j2] = y2r - y6i; + a[j2 + 1] = y2i - y6r; + a[j6] = y2r + y6i; + a[j6 + 1] = y2i + y6r; + } + } + else if ((l << 2) == n) + { + for (j = 0; j < l; j += 2) + { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + x0r = a[j] + a[j1]; + x0i = -a[j + 1] - a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = -a[j + 1] + a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + a[j] = x0r + x2r; + a[j + 1] = x0i - x2i; + a[j2] = x0r - x2r; + a[j2 + 1] = x0i + x2i; + a[j1] = x1r - x3i; + a[j1 + 1] = x1i - x3r; + a[j3] = x1r + x3i; + a[j3 + 1] = x1i + x3r; + } + } + else + { + for (j = 0; j < l; j += 2) + { + j1 = j + l; + x0r = a[j] - a[j1]; + x0i = -a[j + 1] + a[j1 + 1]; + a[j] += a[j1]; + a[j + 1] = -a[j + 1] - a[j1 + 1]; + a[j1] = x0r; + a[j1 + 1] = x0i; + } + } + } + + void cft1st(int n, base_t* a, base_t* w) + { + int j, k1; + base_t wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i; + base_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i; + + wn4r = w[2]; + x0r = a[0] + a[2]; + x0i = a[1] + a[3]; + x1r = a[0] - a[2]; + x1i = a[1] - a[3]; + x2r = a[4] + a[6]; + x2i = a[5] + a[7]; + x3r = a[4] - a[6]; + x3i = a[5] - a[7]; + y0r = x0r + x2r; + y0i = x0i + x2i; + y2r = x0r - x2r; + y2i = x0i - x2i; + y1r = x1r - x3i; + y1i = x1i + x3r; + y3r = x1r + x3i; + y3i = x1i - x3r; + x0r = a[8] + a[10]; + x0i = a[9] + a[11]; + x1r = a[8] - a[10]; + x1i = a[9] - a[11]; + x2r = a[12] + a[14]; + x2i = a[13] + a[15]; + x3r = a[12] - a[14]; + x3i = a[13] - a[15]; + y4r = x0r + x2r; + y4i = x0i + x2i; + y6r = x0r - x2r; + y6i = x0i - x2i; + x0r = x1r - x3i; + x0i = x1i + x3r; + x2r = x1r + x3i; + x2i = x1i - x3r; + y5r = wn4r * (x0r - x0i); + y5i = wn4r * (x0r + x0i); + y7r = wn4r * (x2r - x2i); + y7i = wn4r * (x2r + x2i); + a[2] = y1r + y5r; + a[3] = y1i + y5i; + a[10] = y1r - y5r; + a[11] = y1i - y5i; + a[6] = y3r - y7i; + a[7] = y3i + y7r; + a[14] = y3r + y7i; + a[15] = y3i - y7r; + a[0] = y0r + y4r; + a[1] = y0i + y4i; + a[8] = y0r - y4r; + a[9] = y0i - y4i; + a[4] = y2r - y6i; + a[5] = y2i + y6r; + a[12] = y2r + y6i; + a[13] = y2i - y6r; + if (n > 16) + { + wk1r = w[4]; + wk1i = w[5]; + x0r = a[16] + a[18]; + x0i = a[17] + a[19]; + x1r = a[16] - a[18]; + x1i = a[17] - a[19]; + x2r = a[20] + a[22]; + x2i = a[21] + a[23]; + x3r = a[20] - a[22]; + x3i = a[21] - a[23]; + y0r = x0r + x2r; + y0i = x0i + x2i; + y2r = x0r - x2r; + y2i = x0i - x2i; + y1r = x1r - x3i; + y1i = x1i + x3r; + y3r = x1r + x3i; + y3i = x1i - x3r; + x0r = a[24] + a[26]; + x0i = a[25] + a[27]; + x1r = a[24] - a[26]; + x1i = a[25] - a[27]; + x2r = a[28] + a[30]; + x2i = a[29] + a[31]; + x3r = a[28] - a[30]; + x3i = a[29] - a[31]; + y4r = x0r + x2r; + y4i = x0i + x2i; + y6r = x0r - x2r; + y6i = x0i - x2i; + x0r = x1r - x3i; + x0i = x1i + x3r; + x2r = x1r + x3i; + x2i = x3r - x1i; + y5r = wk1i * x0r - wk1r * x0i; + y5i = wk1i * x0i + wk1r * x0r; + y7r = wk1r * x2r + wk1i * x2i; + y7i = wk1r * x2i - wk1i * x2r; + x0r = wk1r * y1r - wk1i * y1i; + x0i = wk1r * y1i + wk1i * y1r; + a[18] = x0r + y5r; + a[19] = x0i + y5i; + a[26] = y5i - x0i; + a[27] = x0r - y5r; + x0r = wk1i * y3r - wk1r * y3i; + x0i = wk1i * y3i + wk1r * y3r; + a[22] = x0r - y7r; + a[23] = x0i + y7i; + a[30] = y7i - x0i; + a[31] = x0r + y7r; + a[16] = y0r + y4r; + a[17] = y0i + y4i; + a[24] = y4i - y0i; + a[25] = y0r - y4r; + x0r = y2r - y6i; + x0i = y2i + y6r; + a[20] = wn4r * (x0r - x0i); + a[21] = wn4r * (x0i + x0r); + x0r = y6r - y2i; + x0i = y2r + y6i; + a[28] = wn4r * (x0r - x0i); + a[29] = wn4r * (x0i + x0r); + k1 = 4; + for (j = 32; j < n; j += 16) + { + k1 += 4; + wk1r = w[k1]; + wk1i = w[k1 + 1]; + wk2r = w[k1 + 2]; + wk2i = w[k1 + 3]; + wtmp = 2 * wk2i; + wk3r = wk1r - wtmp * wk1i; + wk3i = wtmp * wk1r - wk1i; + wk4r = 1 - wtmp * wk2i; + wk4i = wtmp * wk2r; + wtmp = 2 * wk4i; + wk5r = wk3r - wtmp * wk1i; + wk5i = wtmp * wk1r - wk3i; + wk6r = wk2r - wtmp * wk2i; + wk6i = wtmp * wk2r - wk2i; + wk7r = wk1r - wtmp * wk3i; + wk7i = wtmp * wk3r - wk1i; + x0r = a[j] + a[j + 2]; + x0i = a[j + 1] + a[j + 3]; + x1r = a[j] - a[j + 2]; + x1i = a[j + 1] - a[j + 3]; + x2r = a[j + 4] + a[j + 6]; + x2i = a[j + 5] + a[j + 7]; + x3r = a[j + 4] - a[j + 6]; + x3i = a[j + 5] - a[j + 7]; + y0r = x0r + x2r; + y0i = x0i + x2i; + y2r = x0r - x2r; + y2i = x0i - x2i; + y1r = x1r - x3i; + y1i = x1i + x3r; + y3r = x1r + x3i; + y3i = x1i - x3r; + x0r = a[j + 8] + a[j + 10]; + x0i = a[j + 9] + a[j + 11]; + x1r = a[j + 8] - a[j + 10]; + x1i = a[j + 9] - a[j + 11]; + x2r = a[j + 12] + a[j + 14]; + x2i = a[j + 13] + a[j + 15]; + x3r = a[j + 12] - a[j + 14]; + x3i = a[j + 13] - a[j + 15]; + y4r = x0r + x2r; + y4i = x0i + x2i; + y6r = x0r - x2r; + y6i = x0i - x2i; + x0r = x1r - x3i; + x0i = x1i + x3r; + x2r = x1r + x3i; + x2i = x1i - x3r; + y5r = wn4r * (x0r - x0i); + y5i = wn4r * (x0r + x0i); + y7r = wn4r * (x2r - x2i); + y7i = wn4r * (x2r + x2i); + x0r = y1r + y5r; + x0i = y1i + y5i; + a[j + 2] = wk1r * x0r - wk1i * x0i; + a[j + 3] = wk1r * x0i + wk1i * x0r; + x0r = y1r - y5r; + x0i = y1i - y5i; + a[j + 10] = wk5r * x0r - wk5i * x0i; + a[j + 11] = wk5r * x0i + wk5i * x0r; + x0r = y3r - y7i; + x0i = y3i + y7r; + a[j + 6] = wk3r * x0r - wk3i * x0i; + a[j + 7] = wk3r * x0i + wk3i * x0r; + x0r = y3r + y7i; + x0i = y3i - y7r; + a[j + 14] = wk7r * x0r - wk7i * x0i; + a[j + 15] = wk7r * x0i + wk7i * x0r; + a[j] = y0r + y4r; + a[j + 1] = y0i + y4i; + x0r = y0r - y4r; + x0i = y0i - y4i; + a[j + 8] = wk4r * x0r - wk4i * x0i; + a[j + 9] = wk4r * x0i + wk4i * x0r; + x0r = y2r - y6i; + x0i = y2i + y6r; + a[j + 4] = wk2r * x0r - wk2i * x0i; + a[j + 5] = wk2r * x0i + wk2i * x0r; + x0r = y2r + y6i; + x0i = y2i - y6r; + a[j + 12] = wk6r * x0r - wk6i * x0i; + a[j + 13] = wk6r * x0i + wk6i * x0r; + } + } + } + + void cftmdl(int n, int l, base_t* a, base_t* w) + { + int j, j1, j2, j3, j4, j5, j6, j7, k, k1, m; + base_t wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i; + base_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i; + + m = l << 3; + wn4r = w[2]; + for (j = 0; j < l; j += 2) + { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + j4 = j3 + l; + j5 = j4 + l; + j6 = j5 + l; + j7 = j6 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + y0r = x0r + x2r; + y0i = x0i + x2i; + y2r = x0r - x2r; + y2i = x0i - x2i; + y1r = x1r - x3i; + y1i = x1i + x3r; + y3r = x1r + x3i; + y3i = x1i - x3r; + x0r = a[j4] + a[j5]; + x0i = a[j4 + 1] + a[j5 + 1]; + x1r = a[j4] - a[j5]; + x1i = a[j4 + 1] - a[j5 + 1]; + x2r = a[j6] + a[j7]; + x2i = a[j6 + 1] + a[j7 + 1]; + x3r = a[j6] - a[j7]; + x3i = a[j6 + 1] - a[j7 + 1]; + y4r = x0r + x2r; + y4i = x0i + x2i; + y6r = x0r - x2r; + y6i = x0i - x2i; + x0r = x1r - x3i; + x0i = x1i + x3r; + x2r = x1r + x3i; + x2i = x1i - x3r; + y5r = wn4r * (x0r - x0i); + y5i = wn4r * (x0r + x0i); + y7r = wn4r * (x2r - x2i); + y7i = wn4r * (x2r + x2i); + a[j1] = y1r + y5r; + a[j1 + 1] = y1i + y5i; + a[j5] = y1r - y5r; + a[j5 + 1] = y1i - y5i; + a[j3] = y3r - y7i; + a[j3 + 1] = y3i + y7r; + a[j7] = y3r + y7i; + a[j7 + 1] = y3i - y7r; + a[j] = y0r + y4r; + a[j + 1] = y0i + y4i; + a[j4] = y0r - y4r; + a[j4 + 1] = y0i - y4i; + a[j2] = y2r - y6i; + a[j2 + 1] = y2i + y6r; + a[j6] = y2r + y6i; + a[j6 + 1] = y2i - y6r; + } + if (m < n) + { + wk1r = w[4]; + wk1i = w[5]; + for (j = m; j < l + m; j += 2) + { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + j4 = j3 + l; + j5 = j4 + l; + j6 = j5 + l; + j7 = j6 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + y0r = x0r + x2r; + y0i = x0i + x2i; + y2r = x0r - x2r; + y2i = x0i - x2i; + y1r = x1r - x3i; + y1i = x1i + x3r; + y3r = x1r + x3i; + y3i = x1i - x3r; + x0r = a[j4] + a[j5]; + x0i = a[j4 + 1] + a[j5 + 1]; + x1r = a[j4] - a[j5]; + x1i = a[j4 + 1] - a[j5 + 1]; + x2r = a[j6] + a[j7]; + x2i = a[j6 + 1] + a[j7 + 1]; + x3r = a[j6] - a[j7]; + x3i = a[j6 + 1] - a[j7 + 1]; + y4r = x0r + x2r; + y4i = x0i + x2i; + y6r = x0r - x2r; + y6i = x0i - x2i; + x0r = x1r - x3i; + x0i = x1i + x3r; + x2r = x1r + x3i; + x2i = x3r - x1i; + y5r = wk1i * x0r - wk1r * x0i; + y5i = wk1i * x0i + wk1r * x0r; + y7r = wk1r * x2r + wk1i * x2i; + y7i = wk1r * x2i - wk1i * x2r; + x0r = wk1r * y1r - wk1i * y1i; + x0i = wk1r * y1i + wk1i * y1r; + a[j1] = x0r + y5r; + a[j1 + 1] = x0i + y5i; + a[j5] = y5i - x0i; + a[j5 + 1] = x0r - y5r; + x0r = wk1i * y3r - wk1r * y3i; + x0i = wk1i * y3i + wk1r * y3r; + a[j3] = x0r - y7r; + a[j3 + 1] = x0i + y7i; + a[j7] = y7i - x0i; + a[j7 + 1] = x0r + y7r; + a[j] = y0r + y4r; + a[j + 1] = y0i + y4i; + a[j4] = y4i - y0i; + a[j4 + 1] = y0r - y4r; + x0r = y2r - y6i; + x0i = y2i + y6r; + a[j2] = wn4r * (x0r - x0i); + a[j2 + 1] = wn4r * (x0i + x0r); + x0r = y6r - y2i; + x0i = y2r + y6i; + a[j6] = wn4r * (x0r - x0i); + a[j6 + 1] = wn4r * (x0i + x0r); + } + k1 = 4; + for (k = 2 * m; k < n; k += m) + { + k1 += 4; + wk1r = w[k1]; + wk1i = w[k1 + 1]; + wk2r = w[k1 + 2]; + wk2i = w[k1 + 3]; + wtmp = 2 * wk2i; + wk3r = wk1r - wtmp * wk1i; + wk3i = wtmp * wk1r - wk1i; + wk4r = 1 - wtmp * wk2i; + wk4i = wtmp * wk2r; + wtmp = 2 * wk4i; + wk5r = wk3r - wtmp * wk1i; + wk5i = wtmp * wk1r - wk3i; + wk6r = wk2r - wtmp * wk2i; + wk6i = wtmp * wk2r - wk2i; + wk7r = wk1r - wtmp * wk3i; + wk7i = wtmp * wk3r - wk1i; + for (j = k; j < l + k; j += 2) + { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + j4 = j3 + l; + j5 = j4 + l; + j6 = j5 + l; + j7 = j6 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + y0r = x0r + x2r; + y0i = x0i + x2i; + y2r = x0r - x2r; + y2i = x0i - x2i; + y1r = x1r - x3i; + y1i = x1i + x3r; + y3r = x1r + x3i; + y3i = x1i - x3r; + x0r = a[j4] + a[j5]; + x0i = a[j4 + 1] + a[j5 + 1]; + x1r = a[j4] - a[j5]; + x1i = a[j4 + 1] - a[j5 + 1]; + x2r = a[j6] + a[j7]; + x2i = a[j6 + 1] + a[j7 + 1]; + x3r = a[j6] - a[j7]; + x3i = a[j6 + 1] - a[j7 + 1]; + y4r = x0r + x2r; + y4i = x0i + x2i; + y6r = x0r - x2r; + y6i = x0i - x2i; + x0r = x1r - x3i; + x0i = x1i + x3r; + x2r = x1r + x3i; + x2i = x1i - x3r; + y5r = wn4r * (x0r - x0i); + y5i = wn4r * (x0r + x0i); + y7r = wn4r * (x2r - x2i); + y7i = wn4r * (x2r + x2i); + x0r = y1r + y5r; + x0i = y1i + y5i; + a[j1] = wk1r * x0r - wk1i * x0i; + a[j1 + 1] = wk1r * x0i + wk1i * x0r; + x0r = y1r - y5r; + x0i = y1i - y5i; + a[j5] = wk5r * x0r - wk5i * x0i; + a[j5 + 1] = wk5r * x0i + wk5i * x0r; + x0r = y3r - y7i; + x0i = y3i + y7r; + a[j3] = wk3r * x0r - wk3i * x0i; + a[j3 + 1] = wk3r * x0i + wk3i * x0r; + x0r = y3r + y7i; + x0i = y3i - y7r; + a[j7] = wk7r * x0r - wk7i * x0i; + a[j7 + 1] = wk7r * x0i + wk7i * x0r; + a[j] = y0r + y4r; + a[j + 1] = y0i + y4i; + x0r = y0r - y4r; + x0i = y0i - y4i; + a[j4] = wk4r * x0r - wk4i * x0i; + a[j4 + 1] = wk4r * x0i + wk4i * x0r; + x0r = y2r - y6i; + x0i = y2i + y6r; + a[j2] = wk2r * x0r - wk2i * x0i; + a[j2 + 1] = wk2r * x0i + wk2i * x0r; + x0r = y2r + y6i; + x0i = y2i - y6r; + a[j6] = wk6r * x0r - wk6i * x0i; + a[j6 + 1] = wk6r * x0i + wk6i * x0r; + } + } + } + } + + void rftfsub(int n, base_t* a, int nc, base_t* c) + { + int j, k, kk, ks, m; + base_t wkr, wki, xr, xi, yr, yi; + + m = n >> 1; + ks = 2 * nc / m; + kk = 0; + for (j = 2; j < m; j += 2) + { + k = n - j; + kk += ks; + wkr = (base_t)0.5 - c[nc - kk]; + wki = c[kk]; + xr = a[j] - a[k]; + xi = a[j + 1] + a[k + 1]; + yr = wkr * xr - wki * xi; + yi = wkr * xi + wki * xr; + a[j] -= yr; + a[j + 1] -= yi; + a[k] += yr; + a[k + 1] -= yi; + } + } + + void rftbsub(int n, base_t* a, int nc, base_t* c) + { + int j, k, kk, ks, m; + base_t wkr, wki, xr, xi, yr, yi; + + a[1] = -a[1]; + m = n >> 1; + ks = 2 * nc / m; + kk = 0; + for (j = 2; j < m; j += 2) + { + k = n - j; + kk += ks; + wkr = (base_t)0.5 - c[nc - kk]; + wki = c[kk]; + xr = a[j] - a[k]; + xi = a[j + 1] + a[k + 1]; + yr = wkr * xr + wki * xi; + yi = wkr * xi - wki * xr; + a[j] -= yr; + a[j + 1] = yi - a[j + 1]; + a[k] += yr; + a[k + 1] = yi - a[k + 1]; + } + a[m + 1] = -a[m + 1]; + } + } + + // isgn = 1 - пямое преобразованиея время -> частота + // isgn = -1 - обратное преобразованиея частота -> время + inline void rdft(int n, base_t* a, int* ip, base_t* w, int isgn = 1) + { + int nw, nc; + base_t xi; + + nw = ip[0]; + if (n > (nw << 2)) + { + nw = n >> 2; + makewt(nw, ip, w); + } + nc = ip[1]; + if (n > (nc << 2)) + { + nc = n >> 2; + makect(nc, ip, w + nw); + } + if (isgn >= 0) + { + if (n > 4) + { + bitrv2(n, ip + 2, a); + cftfsub(n, a, w); + rftfsub(n, a, nc, w + nw); + } + else if (n == 4) + { + cftfsub(n, a, w); + } + xi = a[0] - a[1]; + a[0] += a[1]; + a[1] = xi; + } + else + { + a[1] = (base_t)0.5 * (a[0] - a[1]); + a[0] -= a[1]; + if (n > 4) + { + rftbsub(n, a, nc, w + nw); + bitrv2(n, ip + 2, a); + cftbsub(n, a, w); + } + else if (n == 4) + { + cftfsub(n, a, w); + } + } + } +} diff --git a/src/utils/fvec/fvec.cpp b/src/utils/fvec/fvec.cpp new file mode 100644 index 0000000..784770b --- /dev/null +++ b/src/utils/fvec/fvec.cpp @@ -0,0 +1,60 @@ +#include "fvec.hpp" +#include +#include +#include +#include + +namespace hr +{ + // Конструктор: создает вектор заданного размера, заполненный значением v + fvec_t::fvec_t(std::size_t size, base_t v) : m_data(size, v) + { + } + + // Возвращает текущий размер вектора + std::size_t fvec_t::size() const + { + return m_data.size(); + } + + // Проверяет, пуст ли вектор + bool fvec_t::empty() const + { + return m_data.empty(); + } + + // Добавляет элемент в конец вектора + void fvec_t::push_back(const base_t& v) + { + m_data.push_back(v); + } + + // Изменяет размер вектора, новые элементы заполняются значением el + void fvec_t::resize(std::size_t new_size, base_t el) + { + m_data.resize(new_size, el); + } + + // Циклический сдвиг вектора: первая половина меняется местами со второй + // Пример: [1,2,3,4,5,6] -> [4,5,6,1,2,3] + // Для нечетных размеров: [1,2,3,4,5,6,7] -> [5,6,7,1,2,3,4] + // Алгоритм: + // 1. Основной сдвиг: меняем местами первую и вторую половины + // 2. Для нечетных размеров: корректируем средний элемент + void fvec_t::shift() + { + std::size_t half = size() / 2; + std::size_t start = half; + + // Для нечетных размеров start будет на 1 больше half + if (2 * half < size()) ++start; + + // Основной сдвиг: меняем первую половину со второй + for (size_t i = 0; i < half; ++i) std::swap(m_data[i], m_data[i + start]); + + // Корректировка для нечетных размеров: сдвигаем средний элемент + if (half != start) + for (std::size_t i = 0; i < half; ++i) + std::swap(m_data[i + start - 1], m_data[i + start]); + } +} diff --git a/src/utils/fvec/fvec.hpp b/src/utils/fvec/fvec.hpp new file mode 100644 index 0000000..1807e47 --- /dev/null +++ b/src/utils/fvec/fvec.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include "utils/using.hpp" + +namespace hr +{ + class fvec_t + { + public: + fvec_t() = default; + fvec_t(std::size_t size, base_t v); + fvec_t(const fvec_t&) = default; + fvec_t& operator=(const fvec_t&) = default; + fvec_t(fvec_t&&) noexcept = default; + fvec_t& operator=(fvec_t&&) noexcept = default; + ~fvec_t() = default; + + public: + // FOR READING + const base_t& operator[](std::size_t index) const { return m_data[index]; } + const base_t* data() const noexcept { return m_data.data(); } + auto begin() const noexcept { return m_data.begin(); } + auto end() const noexcept { return m_data.end(); } + + public: + // FOR CHANGING + base_t* data() noexcept { return m_data.data(); } + std::vector& src() noexcept { return m_data; } + base_t& operator[](std::size_t index) { return m_data[index]; } + auto begin() noexcept { return m_data.begin(); } + auto end() noexcept { return m_data.end(); } + + public: + std::size_t size() const; + bool empty() const; + void push_back(const base_t& v); + void resize(std::size_t new_size, const base_t el); + void shift(); + + private: + std::vector m_data; + }; +} + diff --git a/src/utils/real_time/real_time.cpp b/src/utils/real_time/real_time.cpp new file mode 100644 index 0000000..fc1e9b0 --- /dev/null +++ b/src/utils/real_time/real_time.cpp @@ -0,0 +1,18 @@ +#include "real_time.hpp" + +namespace hr +{ + real_time::real_time(int s, int n) : m_sec{ s }, m_nsec{ n } + { + normalize(); + } + + real_time real_time::frame2rt(long frame, unsigned int sample_rate) + { + if (frame < 0) return -frame2rt(-frame, sample_rate); + int sec = int(frame / long(sample_rate)); + frame -= sec * long(sample_rate); + int nsec = (int)((double(frame) / double(sample_rate)) * ONE_BILLION + 0.5); + return real_time(sec, nsec); + } +} diff --git a/src/utils/real_time/real_time.hpp b/src/utils/real_time/real_time.hpp new file mode 100644 index 0000000..55f9610 --- /dev/null +++ b/src/utils/real_time/real_time.hpp @@ -0,0 +1,48 @@ +#pragma once + +namespace hr +{ + class real_time + { + public: + real_time(): m_sec{ 0 }, m_nsec{ 0 } {} + real_time(int s, int n); + real_time(const real_time &r) : m_sec{r.m_sec }, m_nsec{ r.m_nsec } { } + + public: + real_time operator+(const real_time &r) const { return real_time(m_sec + r.m_sec, m_nsec + r.m_nsec); } + real_time operator-(const real_time &r) const { return real_time(m_sec - r.m_sec, m_nsec - r.m_nsec); } + real_time operator-() const { return real_time(-m_sec, -m_nsec); } + bool operator <(const real_time &r) const { if (m_sec == r.m_sec) return m_nsec < r.m_nsec; else return m_sec < r.m_sec; } + + public: + static real_time frame2rt(long frame, unsigned int sample_rate); + + public: + int sec() const { return m_sec; } + int msec() const { return m_nsec / 1'000'000; } + int nsec() const { return m_nsec; } + + public: + static const int ONE_BILLION { 1'000'000'000 }; + static const real_time zero; + + private: + int m_sec; + int m_nsec; + + private: + void normalize() { + if (m_nsec >= ONE_BILLION) + { + m_sec += m_nsec / ONE_BILLION; + m_nsec %= ONE_BILLION; + } + else if (m_nsec < 0) + { + m_sec += (m_nsec / ONE_BILLION) - 1; // Учитываем округление + m_nsec = (m_nsec % ONE_BILLION + ONE_BILLION) % ONE_BILLION; // Приводим к положительному значению + } + } + }; +} diff --git a/src/utils/using.hpp b/src/utils/using.hpp new file mode 100644 index 0000000..9dc7a62 --- /dev/null +++ b/src/utils/using.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include + +namespace hr +{ + using base_t = float; + using uint_t = unsigned int; + + // HERE + // убрать это чтоб не было желяния превраить это по аналогии с fvec_t + using ivec_t = std::vector; +} + diff --git a/src/utils/var.hpp b/src/utils/var.hpp new file mode 100644 index 0000000..cea4be7 --- /dev/null +++ b/src/utils/var.hpp @@ -0,0 +1,9 @@ +#pragma once + +#include "noincl.hpp" // IWYU pragma: keep + +namespace hr::var +{ +} + + diff --git a/src/utils/windows/hann/hann.cpp b/src/utils/windows/hann/hann.cpp new file mode 100644 index 0000000..f92986d --- /dev/null +++ b/src/utils/windows/hann/hann.cpp @@ -0,0 +1,13 @@ +#include "hann.hpp" +#include + +namespace hr::math +{ + void hann::creaate(size_t size) + { + m_size = size; + m_win.resize(m_size, 0.0); + for (size_t i = 0; i < m_size; ++i) + m_win[i] = 0.5 * (1.0 - std::cosf((2.0 * std::numbers::pi * i) / m_size)); + } +} diff --git a/src/utils/windows/hann/hann.hpp b/src/utils/windows/hann/hann.hpp new file mode 100644 index 0000000..1b0441d --- /dev/null +++ b/src/utils/windows/hann/hann.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include +#include +#include "utils/fvec/fvec.hpp" + +namespace hr::math +{ + class hann + { + public: + template + void apply(T& target) const + { + size_t length = hack::math::min(target.size(), m_size); + for (size_t i = 0; i < length; ++i) target[i] *= m_win[i]; + } + + // это применяет окно сомвестно с еще одним вектором + template + void apply(T& target, T& base) const + { + uint_t length = hack::math::min(target.size(), hack::math::min(base.size(), m_size)); + for (size_t i = 0; i < length; ++i) target[i] = base[i] * m_win[i]; + } + + void creaate(size_t size); + + public: + size_t m_size; + fvec_t m_win; + }; +} diff --git a/src/utils/workers/plugin.hpp b/src/utils/workers/plugin.hpp new file mode 100644 index 0000000..670f825 --- /dev/null +++ b/src/utils/workers/plugin.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include "utils/cvec/cvec.hpp" +#include "utils/real_time/real_time.hpp" +#include "setup.hpp" +#include "result.hpp" + +namespace hr +{ + class plugin + { + public: + plugin(const setup& st) : m_setup { st } { } + + public: + virtual ~plugin() { } + + public: + // in - результат БПФ + // base - базовый сигнал + virtual void process(fvec_t& base, real_time timestamp){}; + virtual void process(cvec_t& fft, fvec_t& base, real_time timestamp){}; + virtual result get_result() = 0; + + public: + setup m_setup; + }; +} diff --git a/src/utils/workers/result.hpp b/src/utils/workers/result.hpp new file mode 100644 index 0000000..c7679be --- /dev/null +++ b/src/utils/workers/result.hpp @@ -0,0 +1,39 @@ +#pragma once + +#include +#include "utils/real_time/real_time.hpp" +#include "utils/fvec/fvec.hpp" +#include + +namespace hr +{ + struct result + { + struct bit + { + real_time m_duration; + fvec_t m_value; + }; + + void set_bit(bit& b) + { + m_data.push_back(b); + } + + bool empty() const + { + bool res = true; + try + { + if (!m_data.empty()) res = m_data.at(0).m_value.empty(); + } + catch(std::exception& e) + { + hack::error()(e.what()); + } + return res; + } + + std::vector m_data; + }; +} diff --git a/src/utils/workers/setup.hpp b/src/utils/workers/setup.hpp new file mode 100644 index 0000000..3e69afe --- /dev/null +++ b/src/utils/workers/setup.hpp @@ -0,0 +1,25 @@ +#pragma once + +#include + +namespace hr +{ + enum class DOMAIN_PLUGIN + { + TIME, + FREQUENSY + }; + + struct setup + { + // Эти данные заполняются из прочитанного файла (sndfile) + int m_sample_rate; + int m_frames; + int m_channels; + + std::filesystem::path m_file; + std::size_t m_block_size = 1'024; + std::size_t m_step_size = 512; + DOMAIN_PLUGIN m_domain = DOMAIN_PLUGIN::FREQUENSY; + }; +} diff --git a/subprojects/hack.wrap b/subprojects/hack.wrap new file mode 100644 index 0000000..2c6aecd --- /dev/null +++ b/subprojects/hack.wrap @@ -0,0 +1,7 @@ +[wrap-git] +url = chatlanin@gitcast.ru:chatlanin/hack.git +revision = master + +[provide] +hack = hack_dep +