Вычисление PDF формы сигнала по его образцам

27

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

введите описание изображения здесь

Наивный способ сделать это:

  1. Разделите аудиофайл на один фрагмент на горизонтальный пиксель в выходном изображении
  2. Рассчитайте гистограмму амплитуд выборки для каждого куска
  3. Постройте гистограмму по яркости в виде столбца пикселей

Это производит что-то вроде этого: введите описание изображения здесь

Это прекрасно работает, если на блоке много выборок, а частота сигнала не связана с частотой дискретизации, но не иначе. Например, если частота сигнала является точной кратной частоте дискретизации, выборки всегда будут происходить с одинаковыми амплитудами в каждом цикле, а гистограмма будет состоять только из нескольких точек, даже если между этими точками существует реальный восстановленный сигнал. Этот синусоидальный импульс должен быть таким же плавным, как левый выше, но это не так, потому что он ровно 1 кГц, и отсчеты всегда происходят вокруг одних и тех же точек:

введите описание изображения здесь

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

Так что мне бы очень хотелось, чтобы был способ рассчитать истинный PDF (вероятность и амплитуда) непрерывного восстановленного сигнала из его цифровых выборок (амплитуда в зависимости от времени). Я не знаю, какой алгоритм использовать для этого. В общем случае PDF функции является производной ее обратной функции .

PDF греха (х): ddxarcsinx=11x2

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

Этот «PDF интерполированных данных» также применим к попытке, которую я сделал для оценки плотности ядра трека GPS. Он должен был иметь форму кольца, но поскольку он смотрел только на образцы и не учитывал интерполированные точки между образцами, KDE больше походил на горб, чем на кольцо. Если образцы все, что мы знаем, то это лучшее, что мы можем сделать. Но образцы не все, что мы знаем. Мы также знаем, что между образцами есть путь. Для GPS не существует идеальной реконструкции Найквиста, как для аудио с ограниченной полосой пропускания, но основная идея все еще применима, с некоторыми догадками в функции интерполяции.

эндолиты
источник
У вас есть пример многозначной функции, которая вас интересует? Вам, вероятно, придется оценивать его по разрезу веток, который наиболее подходит для ваших физических данных.
Lorem Ipsum
Вас больше интересуют способы рисования такого графика, или этот сюжет является лишь мотивом для вопроса о расчете PDF?
обработке данных
@yoda: Ну, функция выше для синусоидальной волны получается путем взятия только половины цикла, инвертирования и получения производной, потому что каждый полупериод имеет тот же PDF, что и следующий. Но чтобы получить значение для всего произвольного аудиосигнала, вы не можете сделать это предположение. Я думаю, что вам нужно разделить его на «ветки», взять PDF каждого из них, и сложить их все вместе?
эндолит
@datageist: Хм. Я заинтересован в том, чтобы нарисовать такой сюжет, но такой сюжет - PDF. Ярлык, который дает такой же или очень похожий результат, в порядке.
эндолит
@endolith, О да, я понимаю. Просто вопрос об акценте (то есть, какие сочетания клавиш являются разумными).
обработке данных

Ответы:

7

Интерполируйте в несколько раз исходную скорость (например, 8-кратную передискретизацию). Это позволяет принять кусочно-линейный сигнал. Этот сигнал будет иметь очень небольшую погрешность по сравнению с бесконечным разрешением, непрерывной синусоидальной (x) / x-интерполяцией формы сигнала.

Предположим, что каждая пара значений передискретизации имеет непрерывную линию от одного значения к другому. Используйте все значения между. Это дает вам один тонкий горизонтальный срез от y1 до y2 для накопления в PDF произвольного разрешения. Каждый прямоугольный срез вероятности должен быть масштабирован до области 1 / nsamples.

Использование линии между сэмплами, а не самих сэмплов, предотвращает «скачкообразный» PDF, даже в том случае, когда существует фундаментальная связь между периодом выборки и осциллограммой.

Марк Боргердинг
источник
Я написал функцию для линейно-интерполированной гистограммы, но она хитрая. Вы знаете существующий код для этого?
эндолит
Линейная интерполяция имеет огромное значение для большинства сигналов даже без передискретизации. Синус 1 кГц теперь выглядит в основном как синус 997 Гц. Вместо просто горизонтальных линий на значениях образца, теперь это горизонтальные полосы цвета между ними. С передискретизацией полосы тоже сглаживаются. Благодаря FFT-ресэмплингу и некоторому перекрытию со смежными блоками я смогу добиться того, чтобы он достиг реальных пиков между сэмплами. Мне нужно ускорить мой интерполированный код гистограммы ...
endolith
Я полностью переписал мой сценарий для этого, и я думаю , что я получил гистограмму и сглаживание права на этот раз: gist.github.com/endolith/652d3ba1a68b629ed328
эндолиты
Последняя версия находится на github.com/endolith/scopeplot
endolith
7

По сути, я бы использовал «случайный пересэмплер» Джейсона Р., который, в свою очередь, является реализацией стохастической выборки yoda на основе предварительно отобранных сигналов.

Я использовал простую кубическую интерполяцию для одной случайной точки между каждыми двумя выборками. Для примитивного звучания синтезатора (затухание от насыщенного, не ограниченного по площади прямоугольного сигнала + даже гармоник до синуса) это выглядит так:

Произвольно-сэмплированный синтезатор PDF

Давайте сравним это с версией с более высокой выборкой,

введите описание изображения здесь

и странный с той же частотой дискретизации, но без интерполяции.

введите описание изображения здесь

Примечательным артефактом этого метода является превышение в квадратичной области, но на самом деле это то, как PDF-сигнал от sinc-фильтра (как я уже говорил, мой сигнал не имеет полосы пропускания) также выглядел бы и представлял воспринимаемую громкость намного лучше чем пики, если бы это был аудиосигнал.

Код (Haskell):

cubInterpolate vll vl v vr vrr vrrr x
    = v*lSpline x + vr*rSpline x
      + ((vr-vl) - (vrr-vll)/4)*ldSpline x
      + ((vrr-v) - (vrrr-vl)/4)*rdSpline x
     where lSpline x = rSpline (1-x)
           rSpline x = x*x * (3-2*x)
           ldSpline x = x * (1 + x*(x-2))
           rdSpline x = -ldSpline (1-x)

                   --  rand list   IN samples  OUT samples
stochasticAntiAlias :: [Double] -> [Double] -> [Double]
stochasticAntiAlias rs (lsll:lsl:lsc:lsr:lsrr:[]) = []
stochasticAntiAlias (r:rLst) (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)
    = ( cubInterpolate lsll lsl lsc lsr lsrr lsrrr r )
          : stochasticAntiAlias rLst (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)

rand list список случайных величин в диапазоне [0,1].

leftaroundabout
источник
1
Выглядит потрясающе +1 для кода на Haskell.
Datageist
Да, он должен превышать значения выборки. На самом деле я планировал также иметь пиковое значение для каждого столбца пикселей, возможно, нарисованное по-разному, основываясь на максимальных пиках между выборками, а не только на максимальных выборках. Такие формы волны, как flic.kr/p/7QAScX, показывают, почему это необходимо.
эндолит
Под «версией с более высокой частотой дискретизации» вы подразумеваете, что она с повышенной частотой дискретизации, но все еще с одинаковой частотой дискретизации? И это синие точки?
эндолит
1
@endolith Во-первых, это просто исходная форма волны, рассчитанная с более высокой частотой дискретизации. По сути, синие точки представляют звук, дискретизированный на частоте 192 кГц, а самые нижние желтые точки представляют наивно выполненную понижающую частоту до 24 кГц. Верхние желтые точки stochasticAntiAliasэтого. Но версия с более высокой частотой действительно одинакова в обоих случаях.
оставлено около
5

Хотя ваш подход является теоретически правильным (и его необходимо слегка изменить для немонотонных функций), чрезвычайно трудно вычислить обратное значение обобщенной функции. Как вы говорите, вам придется иметь дело с точками ветвления и срезами ветвей, что выполнимо, но серьезно вы бы этого не хотели.

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

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

f(x)=sin(20πx)+sin(100πx)fs=1000fN=1001000 выборок (равномерное распределение) в секунду (я здесь не использую Гц, потому что это подразумевает другое значение) в течение 30 секунд дает график справа (тот же биннинг).

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

введите описание изображения здесь

Лорем Ипсум
источник
1
«вычислить обратное для обобщенной функции крайне сложно». Ну, это не столько функция, сколько последовательность выборок, поэтому для нахождения инверсии нужно просто поменять координаты x и y выборок, а затем выполнить повторную выборку для соответствия новая система координат. Я не могу изменить выборку в любом случае. Мы говорим о ранее существующих данных, созданных с использованием единой выборки.
эндолит
4

Оценка плотности ядра

Одним из способов оценки PDF формы волны является использование оценки плотности ядра .

x(n)K(x)δ(xx(n))P^

P^(x)=n=0NK(xx(n))

Обновление: интересная дополнительная информация.

x(n)n=0,1,...,N1X(k)

X(k)=n=0N1x(n)eȷ2πnk/N

X(k)eȷ2πnk/N

x(n)=1Nk=0N1X(k)eȷ2πnk/N

Итак, угадайте, что вы можете сделать после того, как объедините все PDF-файлы каждого компонента Фурье вместе:

|X(k)|11x2

X(k)x(n)

Однако требуется больше мыслей!

Питер К.
источник
Я думал об этом, но оценка плотности используется для оценки неизвестной функции плотности вероятности. Из-за теоремы отсчетов Найквиста вся форма волны точно известна, и точная функция плотности вероятности тоже должна быть известна. Я согласен с оценкой, если это компромисс между скоростью и точностью, но должен быть способ извлечь из него реальный PDF. Например, восстановленная форма волны может быть получена путем помещения функции sinc в каждую выборку и суммирования их вместе. Может ли PDF быть создан с использованием PDF функции sinc в качестве ядра? Я не думаю, что это работает так.
эндолит
Мол, я не думаю, что это решает проблему, когда выборки сигнала являются кратными частоте дискретизации. Это не учитывает восстановленную форму волны между образцами, не так ли? Это просто размывает каждую точку в PDF, чтобы попытаться заполнить пробелы. У меня была похожая проблема с попыткой оценки плотности ядра трассировки GPS, потому что она не учитывает значения между выборками.
эндолит
4

Как вы указали в одном из ваших комментариев, было бы привлекательным иметь возможность вычислять гистограмму восстановленного сигнала, используя только выборки и PDF-функцию sinc, которая интерполирует сигналы с ограниченной полосой частот. К сожалению, я не думаю, что это возможно, потому что гистограмма sinc не имеет всей информации, которую имеет сам сигнал; вся информация о позициях во временной области, где встречается каждое значение, теряется. Это делает невозможным моделирование того, как масштабированные и задержанные по времени версии sinc будут суммироваться вместе, что и нужно для расчета гистограммы «непрерывной» или дискретизированной версии сигнала без фактического выполнения повышающая дискретизация.

Я думаю, что вы оставили интерполяцию как лучший вариант. Вы указали пару проблем, которые помешали вам сделать это, и я думаю, что они могут быть решены:

  • Затраты на вычисления: это, конечно, всегда относительная проблема, в зависимости от конкретного приложения, для которого вы хотите его использовать. На основании ссылки, которую вы разместили в галерее визуализаций, которые вы собрали, я предполагаю, что вы хотите сделать это для визуализации аудиосигналов. Если вы заинтересованы в этом для приложений в режиме реального времени или в автономном режиме, я бы посоветовал вам создать прототип эффективного интерполятора и посмотреть, действительно ли он слишком дорог. Многофазная повторная выборка - хороший способ сделать это гибким (вы можете использовать любой рациональный фактор).

  • π

Джейсон Р
источник
Но что если сигнал находится на частоте 44,1 / π кГц? :) Это хороший совет, хотя. Есть ли такая вещь, как случайная повторная выборка? Или, на самом деле, я думаю, что идеально сработало бы повторное сэмплирование, чтобы новые сэмплы идеально помещались в ячейки в измерении y, а не были равномерно распределены в измерении x. Не уверен, есть ли способ сделать это
эндолит
2
Вы можете легко реализовать «случайный» ресэмплер, используя структуру Фарроу. Это схема, которая допускает произвольную задержку дробной выборки путем интерполяции с использованием полиномов (часто кубических). Вы могли бы поддерживать фазовый накопитель между выборками , подобный тому, который используется в NCO , который увеличивается на псевдослучайные доли интервала выборки для каждой выходной (повторно дискретизированной) выборки. Значение аккумулятора используется в качестве входа для интерполятора Фарроу, определяющего величину дробной задержки для каждого выхода.
Джейсон Р
Хм, чтобы уточнить, Фарроу - просто оптимизированная для процессора / памяти версия обычной старой полиномиальной интерполяции?
эндолит
1
Да. Это просто эффективная структура для реализации произвольной дробной задержки на основе полиномов.
Джейсон Р.
Тем не менее, кубическая интерполяция является лишь приближением. Я хочу знать истинные пики между выборками, и они не очень хорошо работают на экстремальных пиках: stackoverflow.com/questions/1851384/… На самом деле, кажется, что бесконечный ряд с разрывом типа [..., -1, 1, -1, 1, 1, -1, 1, -1, ...] будет производить бесконечный пик между выборками, поэтому я не уверен, насколько это будет иметь значение на практике.
эндолит
0

Вам необходимо сгладить гистограмму (это даст результаты, аналогичные использованию метода ядра). Как именно должно выполняться сглаживание, требует экспериментов. Может быть, это также можно сделать путем интерполяции. В дополнение к сглаживанию, я полагаю, вы также получите улучшенные результаты, если увеличите частоту дискретизации таким образом, чтобы частота дискретизации была «значительно выше», чем самая высокая частота на входе. Это должно помочь в «сложном» случае, когда синусоида связана с частотой дискретизации таким образом, что в гистограмме заполняются всего несколько элементов. Если довести до крайности, достаточно высокая частота дискретизации должна давать хорошие графики без сглаживания. Таким образом, повышение частоты дискретизации в сочетании с некоторым типом сглаживания должно дать лучшие графики.

Вы приводите пример тона 1 кГц, где график не такой, как вы ожидаете. Вот мое предложение (код Matlab / Octave)

pixels_vertical = 100;
% This needs to be tuned to your configuration and acceptance
upsampling_factor = 16*(pixels_vertical/100); 
fs_original = 48000;
fsine = 1000; % in Hz
fs_up = upsampling_factor*fs_original;
duration = 1; % in seconds
x = sin(2*pi*fsine*[0:duration*fs_up]/fs_up);
period_in_samples = fs_up/fsine;
hist_points = linspace(-1,1,pixels_vertical);
istart = 1;
iend   = period_in_samples;
pixel_values = hist(x(istart:iend), hist_points);
% smooth pixel values
[b,a] = butter(2,0.2);
pixel_values_smooth = filtfilt(b,a,pixel_values);
figure;hold on;
plot(hist_points, pixel_values);
plot(hist_points, pixel_values_smooth,'r');

Для вашего тона 1000 Гц вы получаете это введите описание изображения здесь

Что вам нужно сделать, это настроить выражение upsampling_factor в соответствии с вашими предпочтениями.

Все еще не уверен на 100%, каковы ваши требования. Но используя вышеупомянутый принцип повышения частоты дискретизации и сглаживания, вы получите это для тона 1 кГц (сделанный с Matlab). Обратите внимание, что в исходной гистограмме есть много бинов с нулевыми попаданиями.

введите описание изображения здесь

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