chatlanin/harmonica
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

initial commit

Browse Source
This commit is contained in:
Андрей Зимин
2026-02-16 19:47:03 +03:00
commit 1e06bb4dc7
27 changed files with 1786 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored Normal file
View File

@@ -0,0 +1,5 @@
build
.cache
subprojects/*
!subprojects/*.wrap
src/noincl.hpp

0
README.md Normal file
View File

23
bin/main.cpp Normal file
View File

@@ -0,0 +1,23 @@
#include <hack/logger/logger.hpp>
#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<hr::plugins::magnitude>(setup);
hack::log()("size:", r.m_data.size());
if (!r.empty())
{
std::vector<float> s;
for (auto p : r.m_data) s.push_back(p.m_value[0]);
// hack::log()(s);
}
}

8
bin/meson.build Normal file
View File

@@ -0,0 +1,8 @@
executable(
meson.project_name(),
'main.cpp',
dependencies : deps,
cpp_args: args,
include_directories : inc
)

39
meson.build Normal file
View File

@@ -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')

29
run.sh Executable file
View File

@@ -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

75
src/adapter/adapter.hpp Normal file
View File

@@ -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<typename Plugin>
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);
}
};
}

95
src/harmonica.hpp Normal file
View File

@@ -0,0 +1,95 @@
#pragma once
#include <sndfile.h>
#include <hack/logger/logger.hpp>
#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<typename Plugin>
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<base_t>(setup.m_channels);
}
// Подготовка к следующей итерации
read = hack::math::min(length, static_cast<uint_t>(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();
}
}

46
src/meson.build Normal file
View File

@@ -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

28
src/plugins/magnitude/magnitude.cpp Normal file
View File

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

21
src/plugins/magnitude/magnitude.hpp Normal file
View File

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

21
src/utils/cvec/cvec.cpp Normal file
View File

@@ -0,0 +1,21 @@
#include "cvec.hpp"
#include <math.h>
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);
}
}

18
src/utils/cvec/cvec.hpp Normal file
View File

@@ -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);
};
}

116
src/utils/fft/fft.hpp Normal file
View File

@@ -0,0 +1,116 @@
#pragma once
#include "utils/using.hpp"
#include "utils/cvec/cvec.hpp"
#include "rdft.hpp"
#include <hack/logger/logger.hpp>
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;
}
};
}

924
src/utils/fft/rdft.hpp Normal file
View File

@@ -0,0 +1,924 @@
#include "utils/using.hpp"
#include <math.h>
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);
}
}
}
}

60
src/utils/fvec/fvec.cpp Normal file
View File

@@ -0,0 +1,60 @@
#include "fvec.hpp"
#include <math.h>
#include <algorithm>
#include <hack/math/math.hpp>
#include <hack/logger/logger.hpp>
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]);
}
}

44
src/utils/fvec/fvec.hpp Normal file
View File

@@ -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<base_t>& 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<base_t> m_data;
};
}

18
src/utils/real_time/real_time.cpp Normal file
View File

@@ -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);
}
}

48
src/utils/real_time/real_time.hpp Normal file
View File

@@ -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; // Приводим к положительному значению
}
}
};
}

14
src/utils/using.hpp Normal file
View File

@@ -0,0 +1,14 @@
#pragma once
#include <vector>
namespace hr
{
using base_t = float;
using uint_t = unsigned int;
// HERE
// убрать это чтоб не было желяния превраить это по аналогии с fvec_t
using ivec_t = std::vector<int>;
}

9
src/utils/var.hpp Normal file
View File

@@ -0,0 +1,9 @@
#pragma once
#include "noincl.hpp" // IWYU pragma: keep
namespace hr::var
{
}

13
src/utils/windows/hann/hann.cpp Normal file
View File

@@ -0,0 +1,13 @@
#include "hann.hpp"
#include <numbers>
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));
}
}

33
src/utils/windows/hann/hann.hpp Normal file
View File

@@ -0,0 +1,33 @@
#pragma once
#include <math.h>
#include <hack/math/math.hpp>
#include "utils/fvec/fvec.hpp"
namespace hr::math
{
class hann
{
public:
template <typename T>
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 <typename T>
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;
};
}

28
src/utils/workers/plugin.hpp Normal file
View File

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

39
src/utils/workers/result.hpp Normal file
View File

@@ -0,0 +1,39 @@
#pragma once
#include <vector>
#include "utils/real_time/real_time.hpp"
#include "utils/fvec/fvec.hpp"
#include <hack/logger/logger.hpp>
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<bit> m_data;
};
}

25
src/utils/workers/setup.hpp Normal file
View File

@@ -0,0 +1,25 @@
#pragma once
#include <filesystem>
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;
};
}

7
subprojects/hack.wrap Normal file
View File

@@ -0,0 +1,7 @@
[wrap-git]
url = chatlanin@gitcast.ru:chatlanin/hack.git
revision = master
[provide]
hack = hack_dep