Чем сжать несколько сумм комплексных экспонент

[ilghiz]

Добрый день,

имеется незадачка, а именно имею на входе с десяток оцифрованных сигналов сумм комплексных экспонент, с амплитудами и фазами, которые медленно меняются относительно скорости оцифровки. Грубо говоря, оцифровываю на 25 мегасемплов в секунду, имею 10 каналов, в каждом - комплексная пара этой суммы:

$$a(t) = \sum_{n=1}^N exp(\alpha_n(t) t + i \beta_n(t) t), ~~ t=0, \dots, k t_d$$

$$\forall n: \alpha_n(t) t_d << 1, \beta_n(t) t_d << \pi,$$

$\alpha_n(t)$ и $\beta_n(t)$ медленно меняются на $t=0, \dots, k t_d$ - правда они могут и скачками меняться, но есть много протяженных участков с постоянной производной или самим значением.

$N$ - в военное время достигает нескольких сотен. Более того, во времени, бывает, что какая-то экспонента из суммы пропадает или появляется новая. Сигнал так устроен...

Если его в лоб оцифровать и не обрабатывать, то при 25 мегасемплах, 10 каналах, двух парах на канал и двух байтах на оцифрованный сигнал имеем 1ГБайт/с трафика - он у меня попадает на плиску (FPGA) и тут его хотелось бы пожать. Также есть и шум, который состоит или из очень слабых экспонент или случайного шума. Шум, понятно, не нужен. Сигнал/шум может в военное время быть меньше единицы :(

Между каналами обычно $\alpha_n(t)$ и $\beta_n(t)$ только незначительно отличаются, но нет четкой зависимости. Грубо говоря, эпсилон ранг 10 столбцов с оцифрованными каналами существенно меньше 10 на почти всех коротких промежудках, но даже если это использовать, то компрессия будет только в 2-3 раза.

Если смоделировать сигнал и взять только сами $\alpha_n(t)$ и $\beta_n(t)$, то каждая меняется на столько медленно, что их можно сохранять с объемом меньше килобайта информации на каждую функцию, то есть реальная информативность потока до мегабайта в секунду, то есть можно сжимать мой поток в 1000 раз без существенных потерь информации.

Хочу как-тот аппроксимировать этот поток, чтобы не потерять суть, но все-таки не тащить гигабайт в секунду.

На ум приходит идея:

оцифровываем блок например на 16К отсчетов, делаем Фурье, если есть четко выраженные пики (амплитуда и экспонента стоят на месте) их сохраняем. Вычитаем из этого Фурье, и делаем обратное преобразование.

Далее остаток сравниваем с несколькими аналогичными блоками и пытаемся аппроксимировать такие блоки друг за другом малым рангом.

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

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

Вдруг у кого-то есть идеи, как такое может называться или статьи или ссылки, пожалуйста, посоветуйте!

Спасибо!

[Евгений Машеров]

Авторегрессия порядка 2N, передавать коэффициенты и отклонения от авторегрессионного прогноза? Коэффициенты по блокам, или же менять на ходу и передавать только разности?

[=SSN=]

ilghiz в сообщении #1451004 писал(а):

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

MELP, CELP ?

[Евгений Машеров]

Считать авторегрессию можно методом Левинсона-Дурбина, сложность $O(n^2)$ по числу операций и $O(n)$ по памяти.
Можно покадрово, для кадра коэффициенты и отклонения от прогноза (тут надо как-то жать, имея в виду, что будет много нулей и малых величин, но иногда отклонения будут велики), можно начать с коэффициентов для первого кадра и адаптировать на лету, передавая только изменения коэффициентов (и отклонения от прогноза).
Комплексная экспонента (я правильно понимаю, что в наличии действительная часть, или снимаются именно комплексные отсчёты?) удовлетворяет разностному уравнению второго порядка, а сумма n экспонент - порядка $2n$

[=SSN=]

Евгений Машеров в сообщении #1451096 писал(а):

Считать авторегрессию можно методом Левинсона-Дурбина, ...

Там же используется вещественное деление.. Делить на FPGA то ещё удовольствие..
В FFT по крайней мере нет деления, и сложность $n\cdot\log(n)$.

[Евгений Машеров]

Только это разные n. Длина кадра в случае Фурье и порядок модели в случае авторегрессии... 16К и 20, соответственно.

[ilghiz]

Спасибо большое за интересные советы!

Евгений Машеров в сообщении #1451087 писал(а):

Авторегрессия порядка 2N, передавать коэффициенты и отклонения от авторегрессионного прогноза? Коэффициенты по блокам, или же менять на ходу и передавать только разности?

авторегрессию уже пробовал на симуллированных, забыл об этом написать в исходном сообщении. С ней есть проблема, много данных теряется. То есть если частота стоит на месте, авторегрессия ее вытягивает, но частота может медленно линейно меняться, и вот тут она не работает, то есть там возникает куча собственных векторов близких по значению и все идет юзом. То есть если так посмотреть, авторегрессия в принципе при коротком окне (100-500 отсчетов), если ее рассматривать абстрактно и не вытягивать частоты, позволяет в несколько раз сжать данные, так как блок в 100 отсчетов хорошо аппроксимируется несколькими предыдущими сдвинутыми блоками, но хочется сжимать сильнее, ну хотя бы в 100 раз.

-- 04.04.2020, 15:42 --

Евгений Машеров в сообщении #1451096 писал(а):

(я правильно понимаю, что в наличии действительная часть, или снимаются именно комплексные отсчёты?)

да у меня именно комплексная пара. Представьте - крутится у вас волчек, а вокруг него 2 пары дифференциальных измерителей со сдвигом 90 градусов, что даем мне $R(\cos(a),\sin(a))$. При необходимости, в FPGA я могу эту пару сразу в $R, a$ кордиком превращать, но пока не вижу смысла, если бы я сразу эксоненты растащил, то да, но они приходят ко мне в виде суммы.

[=SSN=]

ilghiz в сообщении #1451004 писал(а):

... оцифровываем блок например на 16К отсчетов, делаем Фурье, ...

А с чем связан выбор именно 16К отсчетов? Если частот несколько сот и не требуется большое разрешение по частоте, то можно же делать FFT на 4096 точек. Это в FPGA уже легко влезает и частота обработки будет выше. Если FFT считать на 500 МГц, то уже 10 каналов по 25 МГц посчитает без проблем.. Если, конечно, сильно не экономить на цене FPGA.

[Евгений Машеров]

Ну, авторегрессия и для комплексных данных работает. Можно попытаться подгонять модель
$y_t=a_1y_{t-1}+b_1y_{t-1}t+\cdots+a_ny_{t-n}+b_ny_{t-n}t$
где члены $yt$ соответствуют дрейфу коэффициентов, но я не знаю быстрого алгоритма для такой модели.
А на какой длине отрезка авторегрессия перестаёт работать?

[ilghiz]

Спасибо за ответы!

=SSN= в сообщении #1451239 писал(а):

ilghiz в сообщении #1451004 писал(а):

... оцифровываем блок например на 16К отсчетов, делаем Фурье, ...

А с чем связан выбор именно 16К отсчетов? Если частот несколько сот и не требуется большое разрешение по частоте, то можно же делать FFT на 4096 точек. Это в FPGA уже легко влезает и частота обработки будет выше. Если FFT считать на 500 МГц, то уже 10 каналов по 25 МГц посчитает без проблем.. Если, конечно, сильно не экономить на цене FPGA.

в принципе в ту плиску, что я пользую, можно запихнуть Фурье на 256К, но у меня несколько каналов, также надобно промежуточные копии хранить, поэтому 16К - это разумный максимум. 4К понятно тоже можно. Как я понимаю, идея с Фурье - не сильно хороша - так как на маленьком окне - маленькая точность, и получу я на нем 300 пиков, и что дальше, а на следующем окне - они слегка другие, и либо городить какой-то реальный огород по угадыванию как это потом доаппроксимировать, или мириться с тем фактом, что сжатие только в 3-5 раз.

-- 05.04.2020, 01:54 --

Спасибо за советы!!!

Евгений Машеров в сообщении #1451289 писал(а):

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

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

Евгений Машеров в сообщении #1451289 писал(а):

Можно попытаться подгонять модель
$y_t=a_1y_{t-1}+b_1y_{t-1}t+\cdots+a_ny_{t-n}+b_ny_{t-n}t$
где члены $yt$ соответствуют дрейфу коэффициентов, но я не знаю быстрого алгоритма для такой модели.

Красивая идея, спасибо!!! Возможно есть красивое решение, но я пока не вывел, хотя сегодня, после Вашего совета, только это и пытаюсь сделать.

Евгений Машеров в сообщении #1451289 писал(а):

А на какой длине отрезка авторегрессия перестаёт работать?

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

[Евгений Машеров]

А порядок модели авторегрессии какой? И в сравнении с числом экспонент.

[ilghiz]

Спасибо за ответ!

Евгений Машеров в сообщении #1451474 писал(а):

А порядок модели авторегрессии какой? И в сравнении с числом экспонент.

так до 300 или около того, пока все компоненты не вычерпаю. Там же проблема в том, что если это сумма экспонент, и порядок существенно меньше, то надо или третью размерность искать, или будет хорошо аппроксимироваться только несколько самых ярких и хорошо отстоящих от других экспонент, а остальные пойдут сильно замешанные.

EDIT: причем, как я говорил, если удлинять число отсчетов, по которой делать авторегрессию, то из-за того, что частоты иногда скачком могут измениться, эти частоты совсем не аппроксимируются. То есть короткое окно нужно для того, чтобы в этом окне не произошел скачек, но, в то же время, в нем мало данных, чтобы вытянуть все мои 300 экспонент в сумме.

[Евгений Машеров]

300 экспонент? Или авторегрессия порядка 300?

[ilghiz]

Евгений Машеров в сообщении #1451499 писал(а):

300 экспонент? Или авторегрессия порядка 300?

по физике там около 2000 экспонент, но из них ярких - около 300, остальные - фактически шум. Мне надо эти 300 вытащить. Я повидимому в терминологии запутался, давайте формулами.

Пусть мой оцифрованный сигнал находится в векторе $a \in R^{n+k}$, где $n$ - размер окна, а $k$ - число сдвигов авторегрессии.

Пусть

$$P = \left(\begin{tabular}{ccccc} 0 & 1 & 0 & $\dots$ \\ 0 & 0 & 1 & \\ $\vdots$ & & & $\vdots$ \\ $\dots$ & 0 & 0 & 1 \\ $\dots$ & 0 & 0 & 0 \\ \end{tabular}\right) \in R^{n+k \times n+k}, ~~ H = = \left( \begin{tabular}{ccccc} 1 & 0 & 0 & $\dots$ \\ 0 & 1 & 0 & \\ $\vdots$ & & & $\vdots$ \\ $\dots$ & 0 & 0 & 1 \\ \end{tabular}\right) \in R^{n+k \times n}$$

Строим $A=(HPa, HP^2a, \dots, HP^ka)$, и решаем $Ax=a$, а $x$ подставляем в коэффициенты полинома и ищем его корни.

Так вот в моем случае, $k$ примерно 300, а наиболее хорошо получаемая аппроксимация достигается при $n=500$.

[ilghiz]

Ой, выше $H$ должна быть вида

$$H = = \left( \begin{tabular}{ccccccc} 1 & 0 & 0 & 0 & 0& $\dots$ \\ 0 & 1 & 0 & & & \\ $\vdots$ & & & & $\vdots$ & $\vdots$ \\ $\dots$ & 0 & 0 & 1 & 0 & $\dots$\\ \end{tabular}\right) \in R^{n+k \times n}$$