Вот код
Точно. Off-by-one.
Формула одна и та же получается. M должно быть равно ASize.

vladisslav2011
Off-by-one.

Там написано:
float M = static_cast<float>(ntaps - 1);
Значит должно быть:
float M = static_cast<float>(ntaps); ?

Только тапками не кидайте, не понимаю я си.
А так да, с M := ASize всё верно.
Просто у меня одной функцией все окна считаются:

procedure SetAnalyzerWnd(Window, wSize: Longint);
var n, j : LongInt;
w : Double;
begin
j := wSize - 1;
for n:= 0 to j do
begin
w := n/wSize;
case Window of
0: AnWindowArr[n] := 1;
1: AnWindowArr[n] := 1-abs(n/(wSize/2)-1); {Triangular (Bartlett)}
2: AnWindowArr[n] := Sin(Pi*w); {Sinus}
3: AnWindowArr[n] := exp(-0.5*Sqr((n-wSize/2)/(0.5*(wSize/2)))); {Gaussian}
4: AnWindowArr[n] := FastSinc((2*n)/wSize-1); {Lanczos}
5: AnWindowArr[n] := 0.54-0.46*Cos(2*Pi*w); {Hamming}
6: AnWindowArr[n] := 0.5-0.5*Cos(2*Pi*w); {Hann}
7: AnWindowArr[n] := 0.42-0.5*cos(2*pi*w)+0.08*cos(4*pi*w); {Blackman}
8: AnWindowArr[n] := 0.355768-0.487396*Cos(2*Pi*w)+0.144232*Cos(4*Pi*w)-0.012604*Cos(6*Pi*w); {Nuttall}
9: AnWindowArr[n] := 1.0-1.93*Cos(2*Pi*w)+1.29*Cos(4*Pi*w)-0.388*Cos(6*Pi*w)+0.028*Cos(8*Pi*w); {Flat Top window}
end;
end;
end;

Programmist
Цикл паскаля for .. to выполняется для последнего значения, а на C написано так
for (int n = 0; n < ntaps; n++)
Для последнего значения не выполняется.
Если бы было
for (int n = 0; n <= ntaps; n++)
был бы цикл паскаля.

это не Бессель, проверки требует
Проверил (двойной расчет факториала в цикле оптимизировал и перевел на с с плюсов):

#include <stdio.h>
#include <math.h>

static double modified_bessel_first_kind_zero_order(double z, int max_k)
{
double ret = 0;
long fact = 1;
//максимум 12, т.к. факториал в квадрате иначе вылетает за long
if(max_k > 12)
max_k = 12;
for(int k = 0; k <= max_k; ++k)
{
if(k)
fact *= k;
long long divider = fact * fact;

double z_base = z * z * 0.25;
double top = pow(z_base, k);
ret = ret + top / (double)divider;
}
return ret;
};

#define IzeroEPSILON 1E-21 /* Max error acceptable in Izero */
static int iloops=0;
static double Izero(double x)
{
double sum, u, halfx, temp;
int n;

sum = u = n = 1;
halfx = x / 2.0;
do {
temp = halfx / (double)n;
n += 1;
temp *= temp;
u *= temp;
sum += u;
} while (u >= IzeroEPSILON * sum);
iloops=n;
return (sum);
}

int main()
{
for(int k=1;k<100;k++)
{
double b=k/10.0;
double c=modified_bessel_first_kind_zero_order(b,12);
double d=Izero(b);
printf("%lf: %lf %lf %lf %d\n",b,c,d,fabs(c-d)*200.0/(c+d),iloops);
}
return 0;
}




./bessel
0.100000: 1.002502 1.002502 0.000000 7
0.200000: 1.010025 1.010025 0.000000 8
0.300000: 1.022627 1.022627 0.000000 9
0.400000: 1.040402 1.040402 0.000000 10
0.500000: 1.063483 1.063483 0.000000 10
0.600000: 1.092045 1.092045 0.000000 11
0.700000: 1.126303 1.126303 0.000000 11
0.800000: 1.166515 1.166515 0.000000 11
0.900000: 1.212985 1.212985 0.000000 12
1.000000: 1.266066 1.266066 0.000000 12
1.100000: 1.326160 1.326160 0.000000 12
1.200000: 1.393726 1.393726 0.000000 13
1.300000: 1.469278 1.469278 0.000000 13
1.400000: 1.553395 1.553395 0.000000 13
1.500000: 1.646723 1.646723 0.000000 14
1.600000: 1.749981 1.749981 0.000000 14
1.700000: 1.863965 1.863965 0.000000 14
1.800000: 1.989559 1.989559 0.000000 14
1.900000: 2.127740 2.127740 0.000000 15
2.000000: 2.279585 2.279585 0.000000 15
2.100000: 2.446283 2.446283 0.000000 15
2.200000: 2.629143 2.629143 0.000000 15
2.300000: 2.829606 2.829606 0.000000 16
2.400000: 3.049257 3.049257 0.000000 16
2.500000: 3.289839 3.289839 0.000000 16
2.600000: 3.553269 3.553269 0.000000 16
2.700000: 3.841651 3.841651 0.000000 17
2.800000: 4.157298 4.157298 0.000000 17
2.900000: 4.502749 4.502749 0.000000 17
3.000000: 4.880793 4.880793 0.000000 17
3.100000: 5.294491 5.294491 0.000000 17
3.200000: 5.747207 5.747207 0.000000 18
3.300000: 6.242630 6.242630 0.000000 18
3.400000: 6.784813 6.784813 0.000000 18
3.500000: 7.378203 7.378203 0.000000 18
3.600000: 8.027685 8.027685 0.000000 18
3.700000: 8.738618 8.738618 0.000000 19
3.800000: 9.516888 9.516888 0.000000 19
3.900000: 10.368958 10.368958 0.000000 19
4.000000: 11.301922 11.301922 0.000000 19
4.100000: 12.323570 12.323570 0.000000 19
4.200000: 13.442456 13.442456 0.000000 19
4.300000: 14.667973 14.667973 0.000000 20
4.400000: 16.010436 16.010436 0.000000 20
4.500000: 17.481172 17.481172 0.000000 20
4.600000: 19.092623 19.092623 0.000000 20
4.700000: 20.858456 20.858456 0.000000 20
4.800000: 22.793678 22.793678 0.000000 20
4.900000: 24.914779 24.914779 0.000000 21
5.000000: 27.239872 27.239872 0.000000 21
5.100000: 29.788855 29.788855 0.000000 21
5.200000: 32.583593 32.583593 0.000000 21
5.300000: 35.648105 35.648105 0.000000 21
5.400000: 39.008788 39.008788 0.000000 21
5.500000: 42.694645 42.694645 0.000000 22
5.600000: 46.737551 46.737551 0.000000 22
5.700000: 51.172535 51.172536 0.000000 22
5.800000: 56.038097 56.038097 0.000000 22
5.900000: 61.376550 61.376550 0.000000 22
6.000000: 67.234407 67.234407 0.000000 22
6.100000: 73.662794 73.662794 0.000000 23
6.200000: 80.717913 80.717913 0.000000 23
6.300000: 88.461552 88.461553 0.000000 23
6.400000: 96.961639 96.961640 0.000000 23
6.500000: 106.292858 106.292858 0.000001 23
6.600000: 116.537324 116.537324 0.000001 23
6.700000: 127.785328 127.785330 0.000001 23
6.800000: 140.136158 140.136160 0.000001 24
6.900000: 153.698994 153.698996 0.000002 24
7.000000: 168.593905 168.593909 0.000002 24
7.100000: 184.952938 184.952944 0.000003 24
7.200000: 202.921322 202.921330 0.000004 24
7.300000: 222.658788 222.658800 0.000005 24
7.400000: 244.341026 244.341043 0.000007 24
7.500000: 268.161288 268.161312 0.000009 25
7.600000: 294.332151 294.332184 0.000011 25
7.700000: 323.087462 323.087508 0.000014 25
7.800000: 354.684471 354.684536 0.000018 25
7.900000: 389.406192 389.406283 0.000023 25
8.000000: 427.563989 427.564116 0.000030 25
8.100000: 469.500432 469.500607 0.000037 25
8.200000: 515.592436 515.592677 0.000047 26
8.300000: 566.254725 566.255056 0.000058 26
8.400000: 621.943634 621.944087 0.000073 26
8.500000: 683.161309 683.161927 0.000090 26
8.600000: 750.460320 750.461160 0.000112 26
8.700000: 824.448747 824.449884 0.000138 26
8.800000: 905.795781 905.797315 0.000169 26
8.900000: 995.237886 995.239948 0.000207 27
9.000000: 1093.585590 1093.588355 0.000253 27
9.100000: 1201.730963 1201.734657 0.000307 27
9.200000: 1320.655848 1320.660768 0.000373 27
9.300000: 1451.440932 1451.447466 0.000450 27
9.400000: 1595.275727 1595.284377 0.000542 27
9.500000: 1753.469570 1753.480991 0.000651 27
9.600000: 1927.463735 1927.478769 0.000780 27
9.700000: 2118.844767 2118.864504 0.000931 28
9.800000: 2329.359178 2329.385016 0.001109 28
9.900000: 2560.929617 2560.963353 0.001317 28

Разница в 0.001317% для 9.9. И не ясно, что точнее.

vladisslav2011
Цикл паскаля for .. to выполняется для последнего значения, а на C написано так
for (int n = 0; n < ntaps; n++)
Для последнего значения не выполняется.
Если бы было
for (int n = 0; n <= ntaps; n++)
был бы цикл паскаля.
Чёрт возьми, голову сломать можно.))
И так с этим Off-by-one всегда проблемы, а здесь запутали вообще в край.

Разница в 0.001317% для 9.9. И не ясно, что точнее.
Вот именно что не ясно. Потому мне моя функция и не нравилась, в виду своего "слишком много букв".
А до факториала я не добрался, устроило то что есть.
Вообще с Бесселем не всё так просто. Для всех его функций нужен десяток модулей, и каждый из них "до потолка".



Ещё бы неплохо с моей сравнить, но это я теперь сам смогу сделать.
Если будет разница более 0.01%, то это уже серьёзно.

vladisslav2011
Спасибо за сравнение, интересно. Я, как уже говорил, просто списал с экрана мат.формулу в код, не особо разбираясь - даже очевидные сейчас первым взглядом оптимизации не делал.

Programmist
Вот именно что не ясно.
Так а откуда ваша реализация появилась? Я где не посмотрю - везде мат.формула через факториал. Вот еще:
https://www.dsprelated.com/freebooks/sasp/Kaiser_Window.html - можете еще погуглить по Kaiser window
Математика ведь первичнее кода - откуда вы брали математику для своей реализации?

wazzoo
Математика ведь первичнее кода - откуда вы брали математику для своей реализации?
Не помню уже, выкроил откуда-то из математических библиотек на паскале.
Там настолько много всего и настолько серьёзно, что не было оснований им не доверять.
Переводить многоэтажные математические формулы на понятный язык я тоже не умею.
Попадёт одна ошибка и всё, приехали, а там их этих ошибок на каждой странице сотня.

можете еще погуглить по Kaiser window
А я тогда не Kaiser гуглил. Все окна были выдернуты откуда-то чуть ли не из википедии.
При переводе с математического языка коса нашла на Кайзер, начал гуглить Бессель, и пошло-поехало...

Programmist
И не ясно, что точнее.
Izero из GNU Radio это тот же факториал, только операции переставлены так, чтобы ошибка округления медленнее накапливалась... Так что по идее должно быть точнее для больших x, так как нет явного ограничения числа циклов. Для маленьких x разница мала, что мой тест и показывает.
Вообще с Бесселем не всё так просто.
Функции высших порядков можно получить из функции нулевого порядка через полином Чебышева. Вот вторая ссылка в гугле: https://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html

Но что-то мы отклонились от первоначальной темы про channelization в сторону оконных функций...

vladisslav2011
чтобы ошибка округления медленнее накапливалась...
Здесь ещё скорость работы надо смотреть.
Как, например, разница в функциях Sinc и FastSinc.
Это намеренно сделано, потому как без Sinc никакая интерполяция вообще не интерполяция, а там их о-го-го сколько.
Пока не сделал Бесселя, на FastSinc всё и крутилось.

Programmist
Вот примерно так будет выглядеть более развёрнутый ответ на Ваш вопрос, за исключением снижения потерь, конечно.
ага, да, хорошо... (ничего не понятно)
Я не программист, мне холивары не интересны, мне интересен результат. А он уже есть. :) Если я правильно понимаю upscale в 16 раз эквивалентен добавалению 4 бит в "разрядность" ЦАП. Что уже прямо хорошо.
Правда всё равно придётся разбираться с форматом "f32le" и писать скрипт для нормализации по амплитуде - общая ару умножает фэдинг. Я всё-таки с КВ сигналами играюсь.

sibrat
(ничего не понятно)
Не удивительно, upscale это и есть интерполяция, задача которой добавить в сигнал то, чего в нём нет, но в теории должно было быть.
Как это правильно сделать, не знает никто, поэтому и существует множество вариантов, от линейной интерполяции до применения различных экзотических функций.
Бессель же вообще был астрономом, и свои функции вывел исходя из многолетних наблюдений за движением небесных тел.
Вот я и подумал, почему бы не применить их в восстановлении недостающих отсчётов (upscale).
И оно сработало, но не в восстановлении сигналов из эфира, а в восстановлении формы звуковой волны, что и является конечной целью всех приёмников.

Я не программист
Я тоже не программист, и тем более не математик, поэтому и подошёл к изучению данного вопроса через "заднюю дверь", с чем мы здесь на форуме и разобрались.

Programmist
Не удивительно, upscale это и есть интерполяция,
upscale это повышение частоты выборок/разрядности. Она может быть с, а может быть без интерполяции. А ещё есть dithering...

P.S.
Cо скриптом я разобрался, результат гораздо лучше АРУ, но, боги, как же это медленно :) Два прохода многогигабайтного файла...

Вопрос к vladisslav2011 в формуле коэффициент считается исходя из максимальной мгновенной амплитуды, а пересчитываем как? перемножение на коэффициент каждой из компонент дало файл в котором максимальный сигнал имеет уровень примерно -15dB.

sibrat
может быть с, а может быть без интерполяции
Две совершенно разные вещи.
ещё есть dithering...
Третья разная вещь, с ней не так давно разбирались в клубе, и из этого вырос новый формат E-PCM.

Всё ещё будет, и Random AES 256 и E-PCM 32 Bit и Bessel Sound и Multi Channel, но потом.
Пока из всего выше перечисленного на SDR вы имеете радужный спектр.
И быстрое I/Q преобразование, у тех кто не пропустил данную разработку мимо ушей 12 лет назад.

Cо скриптом я разобрался, результат гораздо лучше АРУ, но, боги, как же это медленно
Для того и нужно программирование, чтобы было быстро, желательно на ассемблере.

sibrat
Если я правильно понимаю upscale в 16 раз эквивалентен добавалению 4 бит в "разрядность" ЦАП.

Не правильно понимаете, хотя в методичках, на которые я давал ссылку, все это есть.

YuriVR
в методичках, на которые я давал ссылку, все это есть.
У человека, который собирается начать учиться программированию, один взгляд на объём методичек отобьёт всё желание.

Programmist
Для этого они пронумерованы, базовая совсем небольшая и для американцев, т.е. наши студенты уже опустились до этого уровня, в самый раз будет.

YuriVR
для американцев
Тем более, представляете что они потом напишут с такой базовой.))
Единственный быстрый путь чему-то научиться, это взять рабочий исходный текст готового проекта и разобраться что к чему.
При помощи компилятора, естественно. А теорию сами нагуглят, в процессе так сказать проб и ошибок.

Сейчас наблюдаю за развитием одного проекта, людьми с образованием для американцев.
Смех сквозь слёзы. Люди говорят много, красиво, но совершенно не понимают о чём.
И почему их проект никогда не будет работать.

Читал тред. Много думал. Половину пока не понял. Всем спасибо за дискуссию, очень интересно.
Понял такое: мне хорошо подойдёт случайно озвученный мной же FFT-подход, потому что каналов я хочу слушать не один, а сразу штук 10 в видимой приёмником полосе. Тогда один раз сделанное дорогое FFT окупается при необходимости слушать несколько каналов.

Более того, возможна оптимизация: режим мониторинга: когда приёмник сначала просто "слушает" уровни в своих каналах, делая дешёвое FFT бинов на 400, где одна-две "палки" будут представлять всю активность в канале. Причём, это даст повышенное временнОе разрешения детектирования активности. При появлении активность хоть в одном канале, будет включаться режим увеличенного частотного разрешения, когда FFT будет резко детализироваться до 4К...40К палок.

1. Сделать FFT. FFT делается для I и Q независимо. Получается 4 массива чисел (Re + Im для I, Re + Im для Q). FFT хорош из-за аккуратной работы с фазой.
2. Умножить бины (результат FFT) нужных мне каналов на окна. На центр каждого нужного мне канала умножением насадить окно шириной с этот канал (какое окно - поэкспериментировать). Ну, совсем так, потому что каналы могут быть смежными и хвосты окон занулят бины соседних полезных каналов. А точнее так: скопировать (в новый спектр, шириной с нужный мне канал) бины каждого интересующего меня канала из "исходного спектра" с учётом веса натянутого на этот канал окна.

Пример:
1. есть I/Q-поток с приёмника, представляющий полосу 460-465 Мгц. Делаю FFT каждых N семплов (1 FFT для I, 1 FFT для Q) (с перекрытием N/2 и каким-то окном на эти N семплов). Получаю N бинов FFT (в виде 4 массивов). От N зависит разрешающая способность фильтрации по частоте, да. Допустим, N = 4000, тогда каждый бин - это полоса 1250 Гц. 10 бинов - канал 12500.
2. Предположим, каждый интересующий меня канал представлен 10 бинами, а всего меня интересуют 8 каналов.
3. Тогда я создаю 8 новых пустых FFT шириной 10 бинов.
4. В каждый такой новый FFT копирую 10 бинов из соответствующего канала. Лучше взять больше бинов и копировать, умножая на какое-то окно.
5. Делаем IFFT (обратное FFT) в каждом (4) и кладём готовый сигнал на выход. Но не совсем на выход, а суммируем текущее содержание выхода с результатом IFFT начиная с некой позиции выхода POS, которую сдвигам на N/2 после этого. Следующий акт вычисления начиная с (1) засуммирует выход этого IFFT уже с новой позицией выхода POS.

Как-то так?

Programmist
На самом деле, хорошая методичка, инженерная, без лишних формул, с картинками и примерами. Напомню, что это мы покупаем их элементную базу, и это их starlink летает выше уровня моря. И это у них имеются издательства типа Wiley. Нашу современную книжку, допустим, по OFDM можете назвать?

:) quadrature up down converter. Кажется же стандартная тема. Даже если и проспали-плохо давали-не проходили - гуглится за 5 минут. "Супергетеродин" в данном случае не очень умно, т.к., я так понял, I-Q последовательность уже есть. complex up-down - это по сути гетеродин на каждый канал отдельно. И Вы почти правильно это делаете в первых 4-х подпунктах первого пункта, единственно что умнее 3 и 4 объединить, очень неумно фильтровать на завышенной частоте.

Ну и с этого коммента всё началось. Было сказано, что это "по науке", а FFT - это дорого. Про объединение 3 и 4 вроде понял (и напрашивалось). Но потом как я понял LFSR был заборот в этом треде сторонниками FFT с аргументами:

1. FFT даст аккуратную фазу. (разве посоветованное LFSR ломало фазу - я не понял?)
2. Когда много каналов - проще 1 раз заплатить дорого за FFT, зато потом оптом легко надёргать дофига каналов.

Я пока так и не понял как точно реализуется quadrature down converter, только начинаю гуглить. Ну хотя почти понятно: https://en.wikipedia.org/wiki/Digital_down_converter

pavelkolodin
1. Сделать FFT. FFT делается для I и Q независимо. Получается 4 массива чисел (Re + Im для I, Re + Im для Q).
Получается 8 массивов, половина из которых "зеркальные", и обрабатывать надо все 8.
Поэтому лучше начать с реального сигнала.
Как-то так?
:) Как-то так.

Для начала надо подготовить инструмент.
1. Организация ввода вывода из файла.
2. Вывод на экран сигнала в частотном и временном домене.
Потом можно говорить о создании фильтров, для одного канала.
Без контроля всех операций, просто так взять и сделать не получится.
Далее можно добавлять каналы и переводить всю конструкцию на работу в реальном времени.

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

YuriVR
хорошая методичка, инженерная
Я скачал, и даже посмотрел кое-что, методичка хорошая, но очень далека от практической реализации.

pavelkolodin
делая дешёвое FFT бинов на 400
После этого дешевого FFT у вас сразу будет выхлоп на желаемой частоте дискретизации в 12.5kHz. Но он вас не устроит из-за перетекания энергии из соседних бинов. Так что придется увеличить частоту дискретизации до 25 или 50kHz применяя то же FFT на 400 бинов с перекрытием и дофильтровать те бины, в которых есть нужный вам сигнал, обычным FIR со сдвигом спектра и децимацией на 2 или 4, который будет достаточно низкого порядка и будет работать на низкой частоте в 12.5kHz.

FFT делается для I и Q независимо
Так ничего не получится. Читайте методички, которые тут выложили. Много думайте. Проведите несколько экспериментов с FFT, чтобы понять, что же там все-таки получается на выходе (подсказка: вам не потребуется обратное преобразование) и почему нельзя обрабатывать IQ по отдельности. В методичках есть красивые схемы вычисления маленького FFT на 8 бинов двумя способами. Просто внимательно посмотрите, что на что и в каком порядке там умножается, и станет ясно, что дает FFT на выходе.

2. Умножить бины (результат FFT) нужных мне каналов на окна.
Читайте методички. Умножение на окно делается до FFT. После FFT умножать на окно уже не нужно. После FFT у вас уже будет нарезка но узкие полосы.

Делаю FFT каждых N семплов (1 FFT для I, 1 FFT для Q)
IQ нельзя разделять. У вас зеркалки вылезут. Читайте методички.

5. Делаем IFFT (обратное FFT) в каждом (4) и кладём готовый сигнал на выход.
Это не потребуется, если прямое FFT будет сделано правильно. Читайте методички.

Вы как-то сильно все усложняете...

сильно все усложняете...
Не, нормально всё, просто описать весь процесс в 10 строк невозможно.

описать весь процесс в 10 строк невозможно.
Предположим каналов будет максимум 10
1. Выделяем буферы: входной кольцевой для наложения окна 2x(размер FFT), входной для FFT, выходной для FFT, входные кольцевые для FIR 2x(порядок FIR) 10 шт, выходные кольцевые для следующих этапов обработки 10 шт
2. Генерируем окно FFT, taps для 4х FIR (предположим, что бины FFT точно попадают на интересующие частоты), инциализируем фазовые аккумуляторы нулями
3. Заполняем входной буфер
4. Накладываем окно с копированием во входной буфер FFT
5. Выполняем FFT
6. Копируем из выходного буфера FFT те бины, в которые попали интересующие сигналы во входные буферы FIR, копируем сэмплы в начало и корректируем указатель, если требуется
7. Если накопилось 4 элемента во входных буферах FIR, выполняем FIR 10 раз по одному разу для каждого буфера, умножаем результаты на фазовые аккумуляторы,помещаем результаты в выходные буферы, фазовые аккумуляторы увеличиваем на соответствующие сдвиги фаз
8. Передвигаем указатель входного буфера на 0.25*(размер FFT) копируем сэмплы в начало и корректируем указатель, если требуется, дописываем на освободившееся место сэмплы, если есть и переходим к шагу 4
9. Если сэмплов больше нет, освобождаем буферы и выходим из thread function

9 шагов...

vladisslav2011
9 шагов...
Основных. И каждый из них разворачивается в дерево с картинками.
Иначе таких окон понакладывают, что на выходе в лучшем случае будет белый шум, а в обычном деление на нуль.)
И вообще задача не понятна, где фильтровать, зачем фильтровать, какой язык, какая операционная система, и т. д.

и почему нельзя обрабатывать IQ по отдельности
Можно. Можно обрабатывать и по отдельности. Просто нафиг не трогать FFT от слова совсем.

Но он вас не устроит из-за перетекания энергии из соседних бинов. Так что придется увеличить частоту дискретизации до 25 или 50kHz применяя то же FFT на 400 бинов с перекрытием и дофильтровать те бины...
Это бред. Топик-заводиле этого не нужно было совсем. Зачем ему проблемы на ровном месте а потом их решение причем неверное и костлявое!?

Что-то пока не увидел еще классический совет для данной задачи – CIC фильтр, первая попавшаяся ссылка - http://www.dsplib.ru/content/cic/cic.html

YuriVR
CIC фильтр
Это ещё зачем? Пока не ясен ответ на вопросы какой язык, какая операционная система?
Бумагу в процессор не вставишь.

pavelkolodin
И ещё хотелось бы знать, что уже имеется на этот счёт, начиная с операций ввода-вывода.

LFSR
трогать FFT от слова совсем
Можно обойтись пачкой FIR фильтров. FFT только позволит увеличить производительность в несколько раз.

их решение причем неверное и костлявое!?
Решение проверенное и рабочее, обрабатывает 64Msps в реальном времени на относительно скромном мобильном процессоре.
Предолжите свой вариант, который реализует сравнимую производительность.
5Msps IQ на входе, 10 полос по ~10kHz на 12.5ksps на выходе с тестами производительности.

YuriVR
CIC это частный случай FIR оптимизированный для целочисленных вычислений на FPGA/ASIC... На процессорах общего назначения... работать конечно тоже будет.

Programmist
какой язык, какая операционная система
А какая разница. Хоть фортран. Математика от этого не меняется.

начиная с операций ввода-вывода
Если читать файлы/раскидывать обработку на большое количество потоков, то просто увеличивать буферы до тех пор, пока производительность не перестанет расти или задержка станет слишком большой.

vladisslav2011
А какая разница. Хоть фортран. Математика от этого не меняется.
Огромная.
Если кому-то можно дать готовый проект для изучения, где уже всё есть, то для кого-то он окажется абсолютно бесполезен.
Включая бесконечное множество томов с математикой.

Если читать файлы/раскидывать обработку на большое количество потоков, то просто увеличивать буферы до тех пор, пока производительность не перестанет расти или задержка станет слишком большой.
Читать файлы тоже можно по-разному. Можно блоками с диска, а можно сразу в оперативку, потом там и крутить.
У меня второй вариант, что накладывает ограничение на размер файла, но значительно выше по скорости.
Благо для моих целей достаточно, а для записей с SDR уже не пойдёт.

sibrat
а пересчитываем как?
Все правильно. Сначала находим максимум, потом вычисляем коэффициент, на втором проходе умножаем каждый сэмпл на коэффициент (вектор на скаляр).
максимальный сигнал имеет уровень примерно -15dB.
Если сигналов в IQ много, то это вполне нормально. Еще возможно где-то в IQ есть широкополсные импульсы в 1..2 сэмпла - можно сначала найти и удалить их, но это еще +1 проход.
Скорость можно увеличить, читая файл большими блоками (16Mb хотя бы), используя SIMD и разбив обработку на несколько потоков, если все равно упираемся в процессор. Если в диск, то попробовать перегонять файл с одного физического диска на другой.