ресурсоемкие софтовые фильтры без ПФ, и вопросы оптимизации ресурсов

[sewerin]

Предлагаю в этой ветке обсуждать примерно такой концепт. Для 4 полосного стерео нарезать полосы, только с помощью умножений на передаточные функции.
То, что преобразование Фурье взаимно однозначно переводит L2 на L2, пусть остаётся в учебниках по матану.))))
Смежные вопросы:
1. Как распределять ресурсы?
2. Как задействовать графический процессор?
3. Выбор типов и порядков фильтров.

[E.Sokol]

Я возможно чего-то не допонял, но "без ПФ" значит "без преобразования Фурье"? Но тогда остаются только IIR фильтры (с бесконечной импульсной характеристикой), а там вроде никакие особенные ресурсы не нужны.

[Alex]

но "без ПФ" значит "без преобразования Фурье"?

Видимо да.

Но тогда остаются только IIR фильтры

А FIR куда делисль?

[E.Sokol]

А FIR куда делисль?

Так FIR и реализуется через Фурье.
Делать его "в лоб" никаких ресурсов не хватит.

[Alex]

Так FIR и реализуется через Фурье.

Далеко не всегда, а еще точнее - в большинстве случаев как раз FIR делается без FFT а простой конволюцией.

Делать его "в лоб" никаких ресурсов не хватит.

Вполне хватит, если не увлекаться десяткам тысяч тапов.
Но похоже что в аудиотехнике, просто соревнование какое-то - кто больше тапов сделает....

[E.Sokol]

Далеко не всегда, а еще точнее - в большинстве случаев как раз FIR делается без FFT а простой конволюцией.

Ну, насколько я знаю, это "большинство случаев" относится к схемотехнике, ибо там это реализуется элементарно и работает моментально.

Вполне хватит, если не увлекаться десяткам тысяч тапов.
Но похоже что в аудиотехнике, просто соревнование какое-то - кто больше тапов сделает....

Ну, тут всё просто. Необходимое количество тапов зависит от требуемой крутизны фильтра, причём чем ниже частота среза, тем это критичнее.

[Alex]

Ну, тут всё просто. Необходимое количество тапов зависит от требуемой крутизны фильтра, причём чем ниже частота среза, тем это критичнее.

Ну это понятно, но делать десятки тысяч тапов даже на низких частотах большой надобности нет.
(а когда она есть, то это решается другими способами, например IIR-децимация-FIR-апсемлинг)

[E.Sokol]

Ну это понятно, но делать десятки тысяч тапов даже на низких частотах большой надобности нет.
(а когда она есть, то это решается другими способами, например IIR-децимация-FIR-апсемлинг)

Опять же, в железе так делается, потому что реализовать схемотехнически преобразование Фурье - намного более сложная задача. Но у IIR фильтров есть один фатальный недостаток - фазолинейность недостижима в принципе. И формы кривых среза ограничены. И чем больше порядок фильтра, тем больше сдвиги фаз.

[Alex]

Опять же, в железе так делается,

Не только.

Но у IIR фильтров есть один фатальный недостаток - фазолинейность недостижима в принципе.

В данном случае, это совершенно неважно, т.к. происходит вне полосы интереса.

[E.Sokol]

В данном случае, это совершенно неважно, т.к. происходит вне полосы интереса.

Не могу согласиться. Фаза не важна для отдельно взятой частоты, это да. Однако есть нюанс. В полосе раздела одни и те же частоты воспроизводятся одновременно двумя динамиками, которые потом акустически складываются. И суммарная амплитуда будет прямо зависеть от сдвига фаз между этими частотами, следовательно, будет нарушаться линейность АЧХ, а это уже нехорошо.

Вот картинка для демонстрации сути проблемы:

С левой стороны частоты (зелёным цветом) в фазе, с правой - в противофазе. Хорошо видно, как меняется их сумма (красным цветом) в зависимости от сдвига фаз.

[Alex]

Не могу согласиться. Фаза не важна для отдельно взятой частоты, это да.

Для всех частот, вне полосы интереса - тоже.

Хорошо видно, как меняется их сумма (красным цветом) в зависимости от сдвига фаз.

Хорошо видно, что ты просто не понял о чем я говорил

[E.Sokol]

Alex, предлагаю вести дискуссию более аргументировано.
О том, что что преобразование Фурье более эффективно для FIR фильтрации, я узнал из статьи википедии Overlap–add method
Кроме того, известный проект BruteFIR также использует преобразование Фурье, о чём прямым текстом написано в его описании, и в чём также можно убедиться, изучив его исходный код.
Приведите теперь Вы, пожалуйста, ссылку на проект, в котором используется

IIR-децимация-FIR-апсемлинг

[Alex]

предлагаю вести дискуссию

Я не собираюсь вести никаких дискуссий, каждый может делать все, что ему заблагорассудится, и как ему хочеться/можется.

P.S. Я это вполне успешно применял не один раз, правда, для других целей.

[E.Sokol]

Ну раз мы всё выяснили, можно перейти к теме обсуждения.

Свёртку нужно писать на ассемблере, тут без вариантов.
В 2 потока, если процессор двух-ядерный, и в 6 – если четырёх-ядерный.
Если нужна максимально возможная точность, то только на FPU, так как его регистры 80-битные. При этом сам фильтр желательно также сгенерировать с 80-битной точностью.
Если достаточно 64 бит точности, то использование SSE2, AVX или FMA даст 1.5-2.5 кратный прирост производительности.

Видеокарту для параллельных вычисления я не использовал, особенно сказать нечего. Это потому, что на тот момент они поддерживали только 32-битную точность, что для меня было недостаточно. Но сейчас вроде на новых картах можно делать вычисления с 64-битной точностью, этот вопрос нужно изучать отдельно.

[E.Sokol]

Для алгоритма свёртки вот в таком виде

    for(int i = 0;i<input_data.size();i++)
        for(int j = 0;j<fir_filter.size();j++)
            output_fir[i+j] += input_data[i] * fir_filter[j];

простейший вариант на ассемблере может быть таким:

Скрытый текст

conv_fpu PROC data:DWORD, data_size:DWORD, fir_filter:DWORD, fir_size:DWORD, output:DWORD
pusha; сохранение всех регистров

mov edi, data;=входной массив
mov ecx, data_size; длина входного массива
lea edx, [edi+ecx*8];=маркер конца входного массива

mov esi, fir_filter;=fir фильтр
mov ecx, fir_size;=длина фильтра
shl ecx, 3; *= 2 в степени 3 (это 8 = размер типа double в байтах)

mov ebx, output;=выходной массив

; предварительное обнуление выходного массива
fldz; загрузка нуля
mov eax, data_size; длина входного массива
lea ebp, [ebx+eax*8]
; примечание:
; так как значение ebp мы изменили,
; обращаться к параметрам функции больше нельзя
lea ebp, [ebp+ecx];
sub ebp, 8;=маркер конца выходного массива
mov eax, ebx;=выходной массив
loop_clear:
fst QWORD PTR [eax]; запись нуля в выходной массив
add eax, 8; += 8 байт
cmp eax, ebp; сравнить
jl loop_clear; переход, если eax < ebp

; примечание: этот цикл можно (и нужно) развернуть
; делается это так
; fst QWORD PTR [eax+0*8]
; fst QWORD PTR [eax+1*8]
; fst QWORD PTR [eax+2*8]
; fst QWORD PTR [eax+3*8]
; add eax, 4*8
; при этом маркер конца должен быть кратным четырём
; слишком сильно разворачивать тоже не имеет смысла,
; 8, максимум 16 итераций достаточно

fstp st(0); освободить регистр с нулём

; свёртка
loop_i:
fld QWORD PTR [edi]; загрузка input_data[I]
xor eax, eax; j=0

loop_j:
fld QWORD PTR [esi+eax]; загрузка input_filter[j]
fmul st(0), st(1); input_filter[j] * input_data[I]
fadd QWORD PTR [ebx+eax]; input_filter[j] * input_data[I] + output_fir[i+j]
fstp QWORD PTR [ebx+eax]; сохранение output_fir[i+j]
add eax, 8; j++
cmp eax, ecx
jl loop_j
; примечание: и этот цикл можно развернуть, но уже более сложным образом, с учётом кэш линий

fstp st(0); освободить регистр
add ebx, 8
add edi, 8; i++
cmp edi, edx
jl loop_i

popa; восстановление всех регистров
ret
conv_fpu ENDP

[свернуть]

Это не тот вариант с максимально возможной точностью, так как накопление суммы умножений происходит в памяти. В c/c++ нет поддержки типа типа extended double (80 бит) по ряду причин, поэтому чтобы делать накопление в 80-битном FPU-регистре, нужен другой, чуть более сложный алгоритм.

Нужно также помнить, что если свёртка производится не над всем массивом входных, а блоками, то при таком алгоритме нужно складывать хвост одного блока с началом другого.

[E.Sokol]

В реализации многопоточной обработки тоже есть свою особенности.

Во первых, это организация данных. Так как у каждого ядра своя кэш-память, необходимо исключить конкурентность за разделяемую память. То есть вариант, при котором одно ядро обрабатывает только чётные элементы массива, а другое ядро – только нечётные – приведёт только к ухудшению производительности.

Во вторых, это синхронизация потоков, что является не самой простой задачей. Ошибки в синхронизации приводят к взаимоблокировкам, которые приводят к зависанию приложения (навсегда). И кстати, когда Windows 95 висла – то не потому, что она глючная, а вследствие именно взаимоблокировок.
Кроме того, использование объектов синхронизации тоже несколько накладно по ресурсам. Чем реже потоки синхронизируются между собой, тем лучше.

[E.Sokol]

Дополнительная сложность в том, что по многопоточному программированию много теории, но мало конкретных, Copy-Paste совместимых примеров.
Для себя я придумал такое универсальное решение:

1. Определяется базовый класс с реализацией нужного количества дополнительных потоков, вызываемых через виртуальные функции (здесь их два):

IThreadLoop2

#pragma once
#include <windows.h>

DWORD WINAPI IThreadLoop2Func1(LPVOID lpThreadParameter);
DWORD WINAPI IThreadLoop2Func2(LPVOID lpThreadParameter);

class IThreadLoop2
{
private:
    HANDLE hThreads[2];
    DWORD dwThreadIDs[2];    

public: 
    bool _started;
    bool thread_must_exit;

    IThreadLoop2::IThreadLoop2()
    {
        hThreads[0] = NULL;
        hThreads[1] = NULL;
        _started = false;
    }

    IThreadLoop2::~IThreadLoop2()
    {
        if(_started) stop(true);
    }

    bool IThreadLoop2::start()
    {
        IThreadLoop2::stop(true);
        thread_must_exit = false;
        //
        hThreads[0] = CreateThread(
            NULL, 
            0,
            IThreadLoop2Func1, 
            (LPVOID) this, 
            CREATE_SUSPENDED, 
            &dwThreadIDs[0]);
        //
        hThreads[1] = CreateThread(
            NULL, 
            0,
            IThreadLoop2Func2, 
            (LPVOID) this, 
            CREATE_SUSPENDED, 
            &dwThreadIDs[1]);
        //
        ResumeThread(hThreads[0]);        
        ResumeThread(hThreads[1]);        
        
        _started = true;
        return true;
    }
    
    bool IThreadLoop2::stop(bool wait)
    {
        thread_must_exit = true;
        if(wait)
        {
            if(hThreads[0] != NULL)
            {
                WaitForSingleObject(hThreads[0],-1);
                WaitForSingleObject(hThreads[1],-1);
                CloseHandle(hThreads[0]);
                CloseHandle(hThreads[1]);
                hThreads[0] = NULL;        
                hThreads[1] = NULL;    
            }
            _started = false;
        }
        return true;
    }
    
    virtual void IThreadLoop2::thread_loop1()
    {
        
    }    
    
    virtual void IThreadLoop2::thread_loop2()
    {
        
    }    
};

DWORD WINAPI IThreadLoop2Func1(LPVOID lpThreadParameter)
{
    IThreadLoop2* wo = reinterpret_cast<IThreadLoop2*>(lpThreadParameter);
    do
    {
        wo->thread_loop1();
    }
    while(!wo->thread_must_exit);
    wo->_started = false;
    return 0;
}

DWORD WINAPI IThreadLoop2Func2(LPVOID lpThreadParameter)
{
    IThreadLoop2* wo = reinterpret_cast<IThreadLoop2*>(lpThreadParameter);
    do
    {
        wo->thread_loop2();
    }
    while(!wo->thread_must_exit);
    wo->_started = false;
    return 0;
}

[свернуть]

2. Для синхронизации - гибридный объект BarrierLock. Если представить, что потоки - это машины, которые едут с разной скоростью, то BarrierLock - это шлагбаум, который синхронизирует начало их движения и остановку.

BarrierLock

#pragma once
#include <windows.h>

class BarrierLock
{
private:
    long threads_count;
    long threads_started;
    long threads_finished;

    HANDLE evt_wait_for_thread_start;
    HANDLE evt_wait_for_thread_finish;

    bool _exit;

public: 
    BarrierLock::BarrierLock(long threads_count=2)
    {
        this->threads_count=threads_count;
        threads_started = 0;
        threads_finished = 0;

        evt_wait_for_thread_start = CreateEvent(NULL,true,false,NULL);
        evt_wait_for_thread_finish = CreateEvent(NULL,true,false,NULL);
        
        _exit=false;
    }

    BarrierLock::~BarrierLock()
    {
        CloseHandle(evt_wait_for_thread_start);        
        CloseHandle(evt_wait_for_thread_finish);        
    }

    void BarrierLock::init(long threads_count)
    {
        this->threads_count=threads_count;
        _exit=false;

        InterlockedExchange(&threads_started,  0);
        InterlockedExchange(&threads_finished, 0);

        SetEvent(evt_wait_for_thread_finish);
        ResetEvent(evt_wait_for_thread_start);
    }

    void BarrierLock::run()
    {
        if(_exit) return;
        
        if(InterlockedIncrement(&threads_started)==threads_count)
        {
            InterlockedExchange(&threads_started,  0);
            InterlockedExchange(&threads_finished, 0);

            ResetEvent(evt_wait_for_thread_finish);
            SetEvent(evt_wait_for_thread_start);
        }
        else
            WaitForSingleObject(evt_wait_for_thread_start, INFINITE);
    }

    void BarrierLock::finish()
    {
        if(_exit) return;
        
        if(InterlockedIncrement(&threads_finished)==threads_count)
        {
            InterlockedExchange(&threads_started,  0);
            InterlockedExchange(&threads_finished, 0);

            ResetEvent(evt_wait_for_thread_start);
            SetEvent(evt_wait_for_thread_finish);
        }
        else
            WaitForSingleObject(evt_wait_for_thread_finish, INFINITE);
    }

    void BarrierLock::sync()
    {
        finish();
        run();
    }

    void BarrierLock::exit()
    {
        _exit=true;
        SetEvent(evt_wait_for_thread_finish);
        SetEvent(evt_wait_for_thread_start);    
    }
};

[свернуть]

3. Далее мы просто наследуем наш класс от IThreadLoop2, в конструкторе/деструкторе запускаем/останавливаем потоки, а также инициализируем BarrierLock для заданного количества потоков.

M_TripleStereo (сокращённый вариант)

#include "IThreadLoop2.h"
#include "BarrierLock.h"
class M_TripleStereo: public IThreadLoop2
{
private:
    BarrierLock barrier;
    double** _buffers;
    FFT* _fft;

public:    
    M_TripleStereo::M_TripleStereo()
    {    
        barrier.init(3);
        IThreadLoop2::start();
    }

    M_TripleStereo::~M_TripleStereo()
    {
        IThreadLoop2::stop(false); // запрос на остановку всех потоков без ожидания
        barrier.exit();
        IThreadLoop2::stop(true); // ожидание остановки всех потоков
    }
    
    // главный метод, выполняется в основном потоке
    void M_TripleStereo::Process(double** buffers, FFT* fft)
    {
        _buffers = buffers;
        _fft =fft;
         double* dataL = buffers[0];
         double* dataR = buffers[1];
         double* dataC = buffers[2];
        //        
        barrier.run();
        fft->transform(dataL);        
        barrier.finish();
        //
        Split(dataL, dataR, dataC, fft->size);
        // 
        barrier.run();
        fft->back_transform(dataL);
        barrier.finish();
    }

    // вспомогательный поток 1 из IThreadLoop2
    void M_TripleStereo::thread_loop1() 
    {
        barrier.run();        
        if(IThreadLoop2::thread_must_exit) return;
        double** buffers = _buffers;
        FFT* fft = _fft;
        double* dataR = buffers[1];
        fft->transform(dataR);
        barrier.finish();
        //
        barrier.run();
        fft->back_transform(dataR);
        barrier.finish();
    }

    // вспомогательный поток 2 из IThreadLoop2
    void M_TripleStereo::thread_loop2() 
    {
        barrier.run();        
        if(IThreadLoop2::thread_must_exit) return;
        double** buffers = _buffers;
        FFT* fft = _fft;
        double* dataC = buffers[2];
        barrier.finish();
        //
        barrier.run();
        fft->back_transform(dataC);
        barrier.finish();
    }
};

[свернуть]

[E.Sokol]

Что касается фильтров, то тут у меня есть некоторые соображения, несколько отличающиеся от общепринятых.

Так как я программист, а не электронщик, то не вижу логики переносить в программы решения из схемотехники, такие как фильтры Баттерворта, Чебышева и прочих. Если кривые для FIR фильтра так и так приходится рисовать при помощи формул, тогда почему бы сразу не вывести формулу, которая наиболее оптимально решает конкретную задачу?

По моей логике, если задача стоит в декомпозиции сигнала на частотные полосы, это значит:

1. Декомпозиция должна быть обратимой, т.е. сложив полосы обратно, мы должны получить исходный сигнал,
2. Кривая спада АЧХ должна быть зеркально симметричной для НЧ и ВЧ фильтров в логарифмическом масштабе.

Под эти критерии подходит фильтр Линквица-Райли. Но s_valentin пояснил мне, что это – не самый лучший выбор, так как у него пологий спад характеристики, а на высоких порядках он «звенит» из-за слишком резкого перегиба. Гораздо более лучшим, по его мнению, является инверсный фильтр Чебышёва, который затухает сначала плавно, а потом резко. А то, что у него есть пульсации в полосе подавления, так это особо не принципиально, потому как их можно увести за -60 децибел, где их никто и не услышит.
Так-то это может быть и так, но добровольно рисовать пульсации, когда их можно НЕ рисовать… Короче, это добавляет третье условие:

3. Кривая спада АЧХ должна максимально гладко переходить из полосы пропускания в полосу подавления (что бы это не значило).

[Alex]

решения из схемотехники, такие как фильтры Баттерворта, Чебышева и прочих.

"фильтры Баттерворта, Чебышева и прочих" - это не "решения из схемотехники" а аппроксимационные полиномы.

[sewerin]

Александр, напишите, пожалуйста, статью-инструкцию. Расширенное руководство пользователя к xVM.
Подробно про то, как организавать аудиофильскую цифрофильтрацию. Какой софт, какие настройки. В Винде и в Линуксе.