Забегая вперёд: слышал, что фильтровать I/Q - не мудрено, просто рассматриваешь каждый I и Q как отдельный сигнал и независимо фильтруешь. Вопрос не в этом.

Есть SDR, откуда валится I/Q полосы 5 мгц. Это 2 канала по 2 байта 5М раз в секунду. В полосе есть несколько каналов, скажем по 12500 Гц. Хочется их выделить, желательно попутно пересемплировав в человеческий низкий битрейт из оригинального 5-мегасемплового.

Цель - из исходного I/Q - потока получить несколько I/Q - потоков с ЧД = 12500, несущих сигнал отдельных интересных нам подполос исходного спектра.

Верным ли будет такой подход:

1. Задизайнить цифровой полосопропускающий фильтр (FIR или IIR - не важно; я слышал IIR - это выгоднее, меньше операций сложения/умножения во время фильтрации, а артефакты не критичны) с полосой пропускания 0...12500.
2. Сдвинуть весь 5-мегагерцовый спектр влево (вниз по частоте) так, чтобы нижняя частота интересующей нас полосы упёрлась в 0. Делается это умножением на какую-то частоту, видимо на частоту равную начальной частоте интересующей нас полосы - как в супергетеродинах. Правда, у меня пробел в знаниях насчёт отрицательных частот, но это наживное)
3. Получившийся на этапе (2) сигнал подвергнуть фильтрации фильтром (1).
4. Результат пересемплировать.

Шаги, начиная с (2) - проделываются для каждой интересующей нас полосы.

Пример: нужна полоса 2000000 - 2012500 Гц
1. Фильтр есть
2. Умножаем сигнал на синус 2000000 Гц.
3. Результат прогоняем через (2)
4. Пересемплируем на ЧД = 12500 (двойную "по Котельникову" - не надо, у нас же I/Q).
5. PROFIT!!!

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

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

pavelkolodin
Подскажите, что в этой схеме плохо
Медленно, криво, и портит фазу.

как сделать вычислительно эффективнее или правильнее
FFT, ничего эффективнее нет.
После прямого преобразования границы требуемой полосы умножаются на половинку окна, во всех полученных четырёх массивах.
Тип и размер окна выбирается в зависимости от требуемых параметров фильтра.
I/Q преобразования и обратно можно делать в рамках того же цикла, но это уже вопросы оптимизации.

FFT, ничего эффективнее нет.
Не упоминайте его в суе.
Автор не может работать с FFT поскольку у него временнОе представление. Пропускать через FFT обнуляя ненужные бины по частоте - неумно. Очень много ненужных операций. FT для области - будет просто DFT.

Плюс, FFT - оконный алгоритм. А up-down - непрерывный.

LFSR
Пропускать через FFT обнуляя ненужные бины по частоте - неумно
Я не говорил обнулять, говорил обнулять плавно, умножением на половинку окна.
FIR, IIR и прочие свёрточные фильтры свёрткой же и решаются, если есть лишнее процессорное время.

FFT - оконный алгоритм. А up-down - непрерывный
Лёгким движением руки брюки превращаются в элегантные шорты.)
Ели нет возможности сделать FFT, то нет и разговора об эффективности.

Уже сколько лет прошло, когда Sinus сделал I/Q преобразование на FFT? А воз и ныне там...

pavelkolodin
Цель - из исходного I/Q - потока получить несколько I/Q - потоков с ЧД = 12500, несущих сигнал отдельных интересных нам подполос исходного спектра.

Все гораздо проще. Качаете старую версию программы SDR Console V2 и там есть программа SDRDataFileDecimator. Выбираете IQ и вырезаете нужный фрагмент на частоте.

melom
Все гораздо проще.
Однозначно! За последние годы написано столько мусора, что лучше даже не начинать.))

pavelkolodin
рассматриваешь каждый I и Q как отдельный сигнал и независимо фильтруешь.
Неверно.

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

FIR с комплексными коэффициентами вам в помощь.

2. Сдвинуть весь 5-мегагерцовый спектр влево (вниз по частоте) так, чтобы нижняя частота интересующей нас полосы упёрлась в 0
Этого можно добиться просто умножив коэффициенты FIR на синусоиду и сдвинув результат умножением на синусоиду уже на низкой частоте после прореживания...

4. Результат пересемплировать.
Этого можно добиться вычисляя FIR c шагом в n сэмплов...

Возьмите GNU Radio и поэкспериментируйте, если не хочется много писать. Можно прямо в GNU Radio companion flowgraph набросать и подобрать параметры для синтезатора фильтров сразу наблюдая спектр во всех интересующих точках.
Если хочется посмотреть реализацию того, что вам нужно: https://github.com/gnuradio/gnuradio/blob/main/gr-filter/lib/freq_xlating_fir_filter_impl.cc

Programmist
FFT, ничего эффективнее нет.
По результатам моих тестов в некоторых случаях FIR оказывается быстрее, да и FFT не позволяет сдвинуть центральную частоту произвольным и независимым от передискретизации образом. И дробное ресэмплирование с помощью FFT не сделать. Все равно придется ресэмплировать /доворачивать банком FIR фильтров на выходе FFT. Да и для получения хорошего внеполосного подавления приходится выполнять FFT с 4-х кратным оверсэмплингом.

vladisslav2011
По результатам моих тестов в некоторых случаях FIR оказывается быстрее
FIR какого порядка? Хотя, не надо порядка, знаю что и 10 и 15 будет быстрей, но это не фильтр, а лишь ослабление сигнала до требуемой величины.
С ужасными фазовыми искажениями.
Умножайте результат прямого преобразования на окно, как я говорил - получите любую характеристику, но без "хвостов" и искажений.

FFT не позволяет сдвинуть центральную частоту произвольным и независимым от передискретизации образом
Как это не позволяет? FFT всё позволяет, даже быструю свёртку делать в частотном домене.

дробное ресэмплирование с помощью FFT не сделать.
И не надо, ресэмплирование одна из самых сложных операций, надо делать один раз, до рабочей частоты ЦАП звука.
Чего обычно никто не делает, а слушают сигналы с низкой частотой дискретизации, доверяя линейной интерполяции звукового устройства. Очень зря.

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

pavelkolodin
https://disk.yandex.ru/d/79YG6CTO3KmyMZ
скачать и прочитать в порядке сортировки папок, ну или хотя бы самую первую

pavelkolodin
Вам советуют ерунду, делайте так как указал я. Это по науке.

Неверно.
??? Это еще почему? Ортогональность = независимость. Фильтруй наздоровье.

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

Спасибо. Позже попробую осознать, что вы написали, и поэкспериментировать.

Programmist
FIR какого порядка?

Ну например возьмем ddc из Gqrx в случае HackRF на 8Msps в моем форке:


python
from gnuradio import filter
len(filter.firdes.low_pass(1,8000,120,400-2*120))-1
120


120 порядок получается.
Для 1..2 каналов работает быстрее, чем FFT. 3 канала и больше - выигрывает FFT.

В стоке:


python
from gnuradio import filter
len(filter.firdes.low_pass(1,8000,120,1000-2*120))-1
24


24 порядок - тут по любому быстрее.

Может быть, конечно, я как-то неправильно FFT готовлю (я использую только наложение окна+прямое преобразование)...
https://github.com/vladisslav2011/gqrx/blob/fft_channelizer/src/dsp/rx_fft.cpp#L560-L582

??? Это еще почему? Ортогональность = независимость. Фильтруй наздоровье.
Как например отфильтровать сигнал на частоте +1MHz с полосой 10kHz с частотой сэмплирования 5MHz обрабатывая I и Q независимо?
Никак. При обработке по отдельности останется неподавленный зеркальный канал на частоте -1MHz.
Можно конечно сначала сдвинуть частоту в центр, тогда симметричный относительно центра фильтр отработает нормально, но если мы уменьшаем частоту в 400 раз, наверное целесообразнее домножать на синусоиду уже после FIR - в 400 раз меньше операций умножения на канал (не считая FIR, порядок которого при переносе центральной частоты не меняется).

vladisslav2011
Может быть, конечно, я как-то неправильно FFT готовлю
Может быть. FFT.
DoubleToRGB там не нужно, это для отрисовки радужного спектра.

Programmist
Всегда хватало двух.
А подавление идеальное, то есть в нуль.
Пользуясь случаем... не подскажете, как рассчитать окно с хорошим подавлением для перекрытия 2. У меня Kaiser с beta=7.5 дает удовлетворительное подавление только с 4-х кратным перекрытием.

+
Посмотрел по ссылке - там delphi, наверное 2.0 с i387 inline assembler. Планы libfftw с AXV будут по любому эффективнее...

vladisslav2011
как рассчитать окно с хорошим подавлением для перекрытия 2
Для перекрытия есть всего одно окно - окно Ханна, никакие другие не годятся.

AXV будут по любому эффективнее...
Не знаю, до сегодняшнего дня конкурентов по скорости работы в готовых приложениях не было.)

там delphi, наверное 2.0
7.0, всё что выше, набор бесполезных функций с глюками.

Я тестировал все окна, которые генерирует gr::fft::widow::build - наилучший результат получился с Kaiser с beta=7.2 ... 7.7 Все остальные варианты не дали желаемого подавления, Hann в том числе. Видимо у меня слишком высокие требования (60db подавления)...
+ Еще хочется минимум пульсций в полосе пропускания - чтобы не было завала по краям. Этого никак добиться пока не удалось. Видимо придется все-таки вешать дополнительный FIR на выходе.

vladisslav2011
Видимо у меня слишком высокие требования (60db подавления)...
Какие 60, там должно быть все 120, и ограничено оно исключительно точностью данных с плавающей точкой, которые в Delphi называются Single.
32 бита, сколько там дБ? Можно сделать 64, но будет работать медленней, и вообще не имеет никакого физического смысла.
Да и код переписывать с ума сойдёшь.

чтобы не было завала по краям
Не реально, не получится без артефактов обратно собрать блоки. Поэтому и нужно не обнуление, а обнуление плавное.
И размер блока чем больше, тем лучше, при 65536 вполне достаточно окна 128.

Видимо придется все-таки вешать дополнительный FIR на выходе.
Вот не надо самодеятельности.
Похоже, Вам придётся испытать ещё не мало приятных ощущений, при правильной реализации FFT, даже завидую белой завистью.)

Programmist
В вашем файле нет реализации LPF/BPF с использованием FFT. Там только реализация FFT на i387 FPU для ряда фиксированных длин вектора. Это как бы неинтересно. Реализации FFT уже писали все, кому не лень. Есть популяные KISSFFT и LibFFTW + коммерческие реализации. LibFFTW меня вполне устраивает.
Мне больше интересно, как устроен фильтр, так как
Не реально, не получится без артефактов обратно собрать блоки.
Вот этот пункт я как-то не понял. Фильтр работает по классической схеме: наложение окна, FFT, умножение на FFT от импульсной характеристики, iFFT, сложение с предыдущим результатом с перекрытием?
Если так, то мой способ не дает искажений и работает минимум в 2 раза быстрее получается, так как ничего обратно собирать не нужно - результат получается сразу после прямого преобразования...

vladisslav2011
С какой операционной системой работаете? Windows есть?

Ubuntu 18.04 x86-64 на данный момент. Есть Ubuntu 16.04 x86, Ubuntu 20.04 x86-64, Ubuntu 22.04 x86-64 в chroot.
Windows нет и не будет (оно мне помножило на 0 всю радость от покупки первого компьютера). MacOS тоже не будет. Оно на пару с brew мне и сегодня неплохо крови попило, да и в прошлом со своим 512k стеком в процессе реализации поддержки новых спутников BeiDou в GNSS-SDR.

vladisslav2011
Windows нет и не будет
В таком случае, мне сложно будет Вам помочь.
FFT фильтр реализован вот так, на основе представленного выше модуля.
RejectWindowingArr и есть то самое маленькое окно размером FadeSize, на которое умножаются края полосы, перед обнулением.
ComplexArray1 и ComplexArray2 - массивы для двойного перекрытия под окном Ханна.
AudioAReal - дополнительный массив (65536) для сборки сигнала обратно во временной домен.
Лишнее надо выбросить.
Исходный код всего проекта, к сожалению, дать не могу.

//
Давно делал, уже не помню точно что там к чему, но проще не получится.
Все результаты надо выводить на спектр, иначе беда...

Programmist
Все нормально. Я паскаль со школы еще пока совсем не забыл. Полный проект не интересует. Тип RejectWindowingArr интересует, но это наверное Array of Single, так?
Тут фильтр с real taps: умножения на комплексные и сигнал не комплексный, так в строке 143 выполняется восстановление зеркалок.
iFFT выполняется в строках 152...153.
decimation не наблюдаю.
Так что это то, о чем я сразу подумал.
Такой фильтр я пробовал - не устроил ни производительностью ни качеством результата. Polyphase filterbank channelizer не устроил производительностью и гибкостью - нужен внешний роутер для подключения демодуляторов к активным каналам и затыкание NULL sink неиспользуемых каналов, что плодит кучу потоков впустую...

vladisslav2011
Тип RejectWindowingArr интересует, но это наверное Array of Single, так?
Совершенно верно, задаётся вот этим.

decimation не наблюдаю.
Её там нет. Мне не понятно, зачем при фильтрации делать децимацию. Не хватает процессора на всю ширину спектра?
не устроил ни производительностью ни качеством результата
Ну, тогда я не знаю, что Вы там делаете, если меня устраивает, при полном неприятии всего остального.

Programmist
зачем при фильтрации делать децимацию
Мы же нарезаем baseband на каналы, а не фильтруем аудио.
Даже если процессора хватает, зачем тратить его время впустую, если можно снизить частоту дискретизации и слушать демодулированный сигнал без шума вентиляторов в фоне и гораздо дольше при работе от батареи.
Если нужна полоса например максимум 240kHz при частоте дискретизации на входе в 8Msps, то логично будет перейти например на 400ksps, в которые 240kHz влезает с запасом, а нагрузка на процессор на последующих этапах обработки (да и на этапе фильтрации тоже) упадет в 20 раз. Что, вероятно, ТС и хочет.
Выбор промежуточной частоты дискретизации это отдельный интересный вопрос. Если выбирать максимально близко к ширине фильтра, то из-за крутых скатов сильно растет порядок фильтра, и нагрузка на процессор. Если наоборот брать побольше, сильно растет нагрузка на последующих этапах обработки.

vladisslav2011
У меня другое устройство. Максимальная частота дискретизации 160 kHz, при разрядности в 32 бита.
Разрядность фейковая, достигается логарифмическим кодированием, так что можете считать 16 бит, но ДД-таки больше 120 дБ.
И это всё, что удаётся пропихнуть без потерь по сети.
Ну и можно пропихнуть 10 таких потоков по разным портам, на большее и не рассчитывал.
Весь набор фильтров имеется, но больше всего ресурсов жрёт децимация, вот это сделано качественно, и на понимание того, как сделать, ушло 8 лет.


P. S.
FFT не позволяет сдвинуть центральную частоту произвольным и независимым от передискретизации образом
Вот, пожалуйста, сдвиг спектра на произвольную частоту, вместе с фильтром в одном цикле FFT.
Для MaxFFTSize = FFTSize = 65536. Фильтр можно выбросить, будет быстрей.
Надо заметить, без малейших фазовых искажений.
// Всё это куски кода с локальной оптимизацией, поэтому требуют новой оптимизации.

В идеале, всё должно делаться на аппаратном уровне, а на компьютер гнать 48 кГц PCM и картинку по UDP.
Для любителей сохранять весь эфир, диск туда же, в аппаратную часть.
И никаких I/Q, место в 2 раза сэкономит, и вопросов будет меньше, какую часть куда фильтровать.
Мне с устройством проще, такой ширины потоков нет, поэтому пишет всё на флэш, как резервный носитель, и на сервер, как основной.

Programmist
пожалуйста, сдвиг спектра на произвольную частоту
Нет. Не на произвольную. На целое число бинов и без коррекции фазы, так что будут искажения...
В строке 10 произвольная частота уничтожается операцией округления.
И опять операции над real samples, которые на самом деле можно обрабатывать раза в полтора быстрее вообще не переходя к комплексным числам и память сэкономить...

И никаких I/Q, место в 2 раза сэкономит, и вопросов будет меньше, какую часть куда фильтровать.
Нет. Котельникова не обманишь. Или хранить 2 сэмпла real на ПЧ, или домножить на ПЧ/4, отфильтровать и хранить 1 сэмпл IQ (в 2 раза больше размером) на вдвое более низкой частоте. Вот если в IQ лежат на самом деле real samples - Q=0, тогда переход к IQ даст выигрыш в 2 раза. У меня была запись в IQ на центральной частоте 0 - перенес на 5MHz FIR фильтром с decimation 2, получил уменьшение занимаемого места в 2 раза.

Для любителей сохранять весь эфир, диск туда же, в аппаратную часть.
Это что, если FFT 4M, которое я хочу видеть на экране, не влезает в плисину, плисину перепаивать и менять 1.25G SFP на 10G SFP, и цепи питания менять, так как жирная плисина жрет больше? Неее. Я лучше программно все сделаю, а железка пусть гонит IQ с минимальной обработкой.

Программа мониторит/пишет на диск несколько каналов по активности?

vladisslav2011
Не на произвольную. На целое число бинов
Более высокая точность достигается увеличением размера блока и соответственно большим временем задержки.
Здесь без вариантов.
без коррекции фазы, так что будут искажения...
Не будет искажений, FFT единственный инструмент, который корректно работает с фазой, других я не знаю.
Фаза двигается иначе, вот так (из той же серии).

Программа мониторит/пишет на диск несколько каналов по активности?
Устройство пишет на флэшку (аварийный резерв при потере связи) и гонит потоком несколько каналов на сервер, сервер пишет всё круглосуточно.
Один сервер поддерживает 10 устройств. Одна мощная серверная машина 10 серверов.
На скриншоте клиент мониторинга сервера, включая архив, максимальное число каналов - 800.
Максимальное воспринимаемое человеком - 8. Через 8-ми канальное ASIO устройство и 8 динамиков.
Как-то так.

Вариантов оптимизации, как сами понимаете, здесь великое множество, но я уже остановился на достигнутом и ничего больше делать не буду.)
Ребята подключали asm модуль FFT к Unix-подобной системе. Говорят, всё работает, но я им не верю.
Модуль X86 работать может, всё - никогда.

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

Programmist

величением размера блока
Это делается с задержкой в размер буфера (1 сэмпл минимум) умножением на синусоиду заданной частоты без FFT. Точность ограничена точностью float. Надо сдвинуть на миллигерц - будет миллигерц.

клиент мониторинга
Ясно. У меня тоже нечно подобное есть. [Видео YouTube - нажмите для просмотра]
Специальных железок не надо. Подходит все, что поддерживается gr-osmosdr/soapysdr.

Через 8-ми канальное ASIO устройство и 8 динамиков.
Хм. У меня до 256 каналов приема (можно поменять 1 define и будет например 1024) сделано в Gqrx через 2 динамика с расстановкой в стереопространстве. 4..5 одновременно говорящих нормально воспринимается. Запись можно включить по открытию squelch. Имя файлы формируется из частоты-даты-времени.
В общем каждая жаба хвалит свое болото...

vladisslav2011
каждая жаба хвалит свое болото...
Да, только болото одно, и проблемы с ним одинаковые.