есть измеренные данные про которые известно, что по ним прошлись скользяшим средним с неким известным размером D.
каким образом можно восстановить исходные данные которые были до усреднения?
из простого на ум приходит только пройтись по этим уже усреднённым данным еще раз таким же скользящим средним и получить оценку ошибки вызванную усреднением, которую потом добавить обратно к первоначальным данным. Оно вроде как работает, но понятно что это хоть и очень простой, но не совсем честный способ.
Как это делать математически правильно?
понятно что можно сделать Фурье, поделить спектр на характеристику скользящего среднего (по сути КИХ фильтра с одинаковыми единичными коэффициентами) и преобразовать обратно. но при делении на нули в том месте где характеристика имеет провалы пожалуй будут проблемы, оно и понятно так как частоты кратные размеру фильтра давятся в ноль и обратному восстановлению не подлежат.
А если без преобразования в частотную область, с какой функцией надо сделать свёртку чтобы получить фильтр обратный скользящему среднему? или даже в общем случае, каким образом преобразовать коэффициенты КИХ фильтра, чтобы произведение исходного фильтра и пробразованного давало 1.

Для решения задачи нужно прежде всего почитать про z-преобразование.

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

и конечно же не спутать с z-преобразованием, которое z-score, или, которое, преобразование фишера, а можно вообще без этого преобразования на пальцах написать, что скользящее среднее - есть умножение матрицы на исходный вектор, и просто найти обратную к этой матрице. Если есть необходимые знания о теплицевости-циркулянтности - можно из этого через Фурье и то самое зет-преобразование самому вывести, а можно в лоб как нижнюю треугольную решить. Тоже кстати достаточно устойчиво и программируется в 6 строчек без сложных слов типа БПФ и комплексных чисел.

Z-преобразование это дискретное преобразование Лапласса.
Умножение матрицы на вектор? Расскажите подробнее пожалуйста. Как с использованием z-преобразования знаю, а про такой подход не слышал. Фактически скользящее среднее есть свертка. По теореме о свертке ее можно свести к умножению проделав z-преобразование или дискретное преобразование фурье (которое в данном случае совпадет с z-преобразованием). Каким боком применить к этому делу матрицу для мне неизвестно. Или вы имеете ввиду алгоритмы быстрого преобразования Фурье?

Мне известен только один способ восстановления последовательности, обработаной сверткой. Он называется оптимальная гомоморфная фильтрация и детально изложен в книге Рабинера и Голда по цифровой обработке сигналов. Книга старая 70-х годов.
Ключевое слово для поиска в гугле кепстр (измененное слово спектр).
http://en.wikipedia.org/wiki/Cepstrum

да, всё оказалось просто, надо посчитать производную.
скользящее среднее:
Y[i] = (X[i-N] + X[i-N+1] + ... +X[i+N]) / 2N
или если выразить через предыдущий отсчёт:
Y[i] = Y[i-1] + (X[i+N] - X[i-N-1]) / 2N
соответственно
X[i+N] = X[i-N-1] + (Y[i] - Y[i-1]) * 2N.

upd: только не 2N, a 2N+1 конечно же.

Возможно прямоугольная форма импульсной характеристики фильтра производящего скользящее среднее позволит что-нибудь упростить.




Поздравляю Как все просто оказалось.

Пусть $x \in R^n$ исходный вектор до обработки скользящим средним. Тогда если скользящее среднее имеет длину K, то результат его применения можно записать как умножение матрицы
(1 0 ... 0)
(1/2 1/2 0 ... 0)
(... {1/K ... 1/K}(K раз) 0 ... 0)
(...)
(0 ... 0 1/K 1/K ... 1/K)
(...)
(0... 0 1)

- надеюсь, понятно, без ЛаТеХа не могу красиво нарисовать.

Теперь внимательно смотрим, что получилось. Если отбросить первые и последние K-1 строчек - получается Теплицева матрица, если перейти к бесконечным, то в пределе по спектру она должна совпадать с соответствующим циркулянтом, спектр которого получается их БПФа оного. Как я понимаю, формулы должны с зет-преобразованием совпасть, но кропотливо не проверял.

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

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

Лучший метод, конечно же через Фурье, так как есть контроль спектра, но с умением дела можно и без оного обойтись.

А какого рода данные там были? Если строго гармонические, то какие частоты сверху и снизу.

Никак.

В общем случае исходные данные теряются навсегда. Если лень изучить теорию обратных операций (хотя интегрирование это, конечно, функционал, а не оператор), то хотя бы поразмышляйте над средним арифметическим (x+y)/2 = a. Отсюда имеем линейное уравнение x+y = 2a. а известно. Сколько решений можно получить для двух переменных из одного уравнения? Правильно.

я решение выше вроде как привёл.
пока спектр сигнала меньше чем первый провал до нуля на характеристике скользящего среднего (1 / длину фильтра) работает нормально.
ну а после умножения на 0, дейтствительно, восстанавить обратно непросто.
Нажмите для просмотра прикрепленного файла
синенький - исходный сигнал, красненький - скользящее среднее по 15 точкам, зелёненький - восстановлен обратно из красненького.

Это не решение, а хотение. Вы _не_ знаете величину X[i-N-1] в общем случае. А то, что у вас в ряду данных первые отсчёты - нули, это случайность. Разность соседних средних = разности _двух_ исходных отсчётов. Опять одно уравнение для _двух_ переменных. Т.е. ваша формула на самом деле

X'[i+N] = ??? + (Y[i] - Y[i-1]) * (2N+1), штрих обозначает оценку X[i+N], а модуль |???| не ограничен сверху, т.е. ошибка может быть произвольно большой.

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

Скользящее среднее - это аппроксимация определённого интеграла в непрерывном случае. Восстановить подинтегральную функцию по значению интеграла - невозможно.

В общем случае - никак.
Для приближенного восстановления требуются оговорки.

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

Методов синтеза много. Если в качестве инструмента Матлаб, то посмотрите fdesign.arbmagnphase

z-плоскость дала бы вам идею того, с чем надо бороться.

Я может что-то не понимаю, но в чем загвоздка? Взятие среднего - пропуск сигнала через фильтр (для третьего порядка) (z^2+z+1)/(3*z^3). Пропустите полученный сигнал через обратный фильтр (3*z^3)/(z^2+z+1) и будет исходный сигнал.

Посмотрите, где лежат полюсы вашего "обратного фильтра".

Топикстартер задал тривиальный вопрос, на который интуитивно нашел правильный ответ.

Скользящее среднее - однородный фильтр
Передаточную ф-цию такого фильтра можно записать как



H(z) = 1/N*(1-z^-N)/(1-Z^-1)

разностное уравнение y(n)=( x(n) - x(n-N) )/N + y(n-1)

Восстанавливающий фильтр будет с пф 1/H(z) = N*(1-z^-1)/(1-z^-N)

разностное уравнение z(n)=( y(n) - y(n-1) )*N + z(n-N) что и получил топикстартер.

При чем тут определенные интегралы и др дивергенции роторов - совершенно непонятно...

Вот скрипт для особо неверующих

Верующие могут отправлять свои естественные культовые надобности в ином месте. Интуиция - хороша с девочками - согласится/несогласится. А формулы - вещь, напротив, аналитическая. В них смотреть надо.
А тут - вроде смотрим, а видим...

Ещё раз, последний. Есть, известны только y(i). _Все_ z(i) давно пропали. Начинаем по вашей "формуле" восстанавливать z(n). y(n), y(n-1) известны, а z(n-N) - нет. Все, процесс закончен, не успев начаться.

Поясняю, при чём тут интеграл. Формула скользящего среднего z = (y(1) + ... + y(n))/n - это вычисление определённого интеграла аппроксимацией прямоугольниками, т.е. нулевого порядка.

Учите численные методы для 1 курса.

Особо неверующим будет вдвойне интересно, если в канале передачи присутствует шум, или если мы попытаемся восстанавливать данные не с начала передачи.

Топикстартер писал где-то про канал с шумом? Или про несначало? Зачем фантазировать?



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

Откуда взялся предыдущий отсчёт?
В той же формуле кол-во элементов X, если это не ошибка, = 2N+2.

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

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

Да. Для непрерывной обработки длинной последовательности / со случайного момента / в условиях шумов этот способ конечно же не годится. Но у топикстартера про все эти частности нет ни слова.



Доброхоты навыдумывали себе кучу условий, назадавали вопросов и принялись бодро на них отвечать.
Отвечать-то надо по существу, а не захламлять топики.

По существу же сказано, если имеются только данные с выхода фильтра, то для вселенского счастья не хватит одного неотфильтрованного элемента

Не согласен. Топикстартер имеет некую отфильтрованную последовательность данных, и хочет из нее получить исходные. С чего начинались исходные данные, не известно. Явно не с нуля. Включается некий прибор, выходит на рабочий режим, запускается измерение сигнала, обрабатывается, сохраняется последняя часть. Вот эта часть и передается. "Восстановлению не подлежит."

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

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

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

С рекурсивным алгоритмом в принципе "нечестно". Исследовать точность/достоверность (приближения) необходимо без рекурсии. В общем случае стартовый неотфильтрованный элемент неизвестен. А проход всей последовательности рекурсивным фильтром автоматом расширяет ширину фильтра до всей длины последовательности.

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

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

При нулевых начальных условиях и наличии прописи старта фильтра - все вычисляется. Но и если нет стартапа, то тоже можно кое-что сделать. Как пример - если в прописи есть участок, на котором в течении N тактов значение суммы не меняется (т.е. текущее значение равно N-му) - мы точно определяем все элементы суммы и дальше можем все точно вычислять.

Отвечу сам себе. Для задачи топикстартера - плохо. Никакой дополнительной точности не даст. Наоборот, вроде как шума добавит.


Подать сигнал калибратора на вход фильтра на некоторое время.

Из чего это следует?
Если на выходе N элементов подряд имеют одинаковое значение, то они являются функцией N+N-1 входных элементов.

Дык, если задаться вопросом "что можно сделать если неизвестны фаза/начальные условия/длина итд?" то ответ тут однозначен - ничего. Именно это и имели ввиду местные телепаты. Объясняется этот факт наличием нулей в частотной характеристике однородного фильтра. Скомпенсировать их можно только полюсами, что в свою очередь требует полной когерентности последовательностей. Так что все фантазии на эту тему останутся фантазиями (реверанс в сторону AndrewN) и не более того.

Скользящее среднее нелинейно искажает высокие частоты. Через фурье или свёртку недостоверно будет. В частоту Fs/2D попадут искажения (обрезки) частот, некратных Fs/2D. Если свернуть выход фильтра с обратной скользящему фильтру функцией, то частота Fs/2D будет плавать по фазе. А она ещё должна будет усилиться в 1.41 раза для восстановления амплитуды. И увеличения амплитуды "мусора".

На единичной окружности.

интегрирование - линейная операция, а скользяцее среднее - вообще КИХ фильтр с одинаковыми коэффициентами, на 0 некоторые частоты помножит, но ничего оно нелинейно не искажает.

Выше дописал. D задаёт интервал, на котором будут искажаться некратные интервалу частоты, причём и выше 1 и ниже. У меня еденичный импульс вылез в области Fs/2D. (PS на коротком интервале) Хотя вопрос встанет - как считать. Если при восстановлении сигнала требуется таки визуальное приближение к оригиналу как достоверность оценки, то плохо дело. А если оценкой будет результат вроде БПФ, то это другой разговор. Спектр по БПФ очень интересная вещь, к тому же двумерная. Может как-то на относительно коротких последовательностях (2..10 раз выше D) можно попытаться минимизировать искажение фазы. Но если сигнал в дальнейшем как-то узкополосно фильтруется, то смысла нет.

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

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

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

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

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

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

PS
Ошибки округления скользящего среднего "на обратном пути" размножаются очень сильно по нарастающей в области ВЧ. А в стартовых значениях даже добавление шума (которого не было в исхонике) 10-20% от размаха сигнала не сильно искажает восстановление. ИМХО можно улучшать разрешение сигнала даже из округлённых значений скользящего среднего. Но с фильтрацией при восстановлении, вероятно, аналогичной предворительному поиску стартовых значений. Проблемы только будут с частотами кратными Fs/(2N+1). Но их фильтр давит без вариантов.

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

Там для счастья не хватает 2N+1 неизвестный элемент исходника. И алгоритм скользящего без потерь, что легко делается.

Пока есть идея взять для стартового поиска элементы с выхода скользящего, как реально приближённое значение к оригиналу. По формуле получается, что известны зависимости между крайними элементами фильтра. Т.о. можно итеративно "раздвигать" их в нужную сторону, одновременно соблюдая постоянную составляющую (т.к. точно известно среднее арифметическое этой и соседних к ней точек). Мелкими шагами, к примеру по 1/20 от отклонения. Есть вероятность, что результат будет сходиться к оригиналу.

Т.е. получается как с алиасом. Среднее по N точкам будет нормально передавать частоты с периодом >= 2N, далее - спады до нуля на частотах с периодом N, N/2, N/3 и т.д. Т.е. что-то правдоподобное можно выжать до первого нуля (период - N), т.е. до удвоенной частоты, а дальше будет каша с полной неопределенностью.

А, с другой стороны, ежели однократно узнать все N членов, то можно восстановить исходные отсчеты точно, т.е. все частоты до периода в 2 такта. Интересная ситуация.

Если считать, что сигнал был негармонический, то определённость более чем достаточная. Эффект алиасов имеет смысл именно к гармоническим сигналам. У негармонических должна быть мера оценки только приближение одномерного результата к исходнику. Что такое алиас применительно к одномерному негармоническому сигналу? Или если сигнал широкополосный и далее умножается/сворачивается с широкополосным сигналом совсем с другими характеристиками. Широкополосные аналоговые фильтры тоже не заточены на ограничения в гармоническом сигнале по сетке дискретизации. Но лучше сразу "договориться", что если источник информации о сигнале является одномерным, то вне зависимости от того, гармонический сигнал или нет, единственный достоверный критерий оценки (и восстановленного тоже) сигналов - их похожесть в одномерном виде. Причём даже спектральное сравнение корректно восстановленного сигнала будет совпадать с оригиналом. Появится вопрос "а был ли мальчик?". Или про суслика.

Отличия каждого соседнего элемента исходника в результате скользящего среднего, особенно без потерь, остаются. То есть воссоздание каждого элемента - возможно, с небольшим сглаживанием, чтобы обрезать (в т.ч. паразитную) Fs/2. А по НЧ/постоянке вся информация сохранилась.

Интересный велосипед тут изобретают...
Автор проговорился, что он подозревает, что датчик осуществляет скользящее среднее, а он хочет узнать "правильный" сигнал.
Так вот. Давным-давно озабоченные спектроскописты используют теорему о свертке (Ну, все знают про Фурье-преобразование от свертки). Только нужно знать аппаратную функцию. Ее экспериментально получают с помощью априорных знаний или подают дельта-функцию. Так вот - предположение о скользящем среднем - это только предположение. Еще нужно правильно обрезать лишнее. Тоже вопрос.

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

Разве скользящее среднее не свертка?
As you may have recognized, calculating the simple moving average is similar to the convolution: in both cases a window is slid along the signal and the elements in the window are summarized.
http://matlabtricks.com/post-11/moving-ave...-by-convolution

Свёртка, конечно. С функцией, где все элементы еденичные.
Порождаемые ею алиасы при увеличении периода сигнала в спектральном представлении исчезают. Но не "бесплатно". А частоты, кратные Fs/(2N+1) (в обозначениях автора темы) вообще будут недостоверны в спектре.

Свертка линейная операция? или "нелинейно искажает высокие частоты".

По методу оптимальной гомоморфной фильтрации все просто:
Из кепстра обработанной последовательности вычитаем кепстр прямоугольного окна и восстанавливаем из результата исходную последовательность, посредством потенцирования и обратного преобразования Фурье.

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


В процессе дискуссии возник интересный диалог "экспериментаторов" с "теоретиками". Почти уверен, что рекурентное решение в силу своей простоты блестяще справится с работой. Ради такого решения несложно добавить в программу процедуру калибровки, которая и получит достоверное значение. Имеются и другие способы передать одно оригинальное значение и его позицию в процедуру восстановления.

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

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

А насчет восстановления... Тут, на самом деле, все еще интереснее и, похоже, однозначного метода восстановления нет.

Как пример. Пусть мы имеем нулевой буфер (длиной N) и подаем синус с периодом N. После N тактов этот синус окажется в буфере и, с N+1 такта, сумма перестанет изменяться. Т.е. на выходе мы имеем константу. Как из нее можно восстановить исходный синус? Т.е. наличие сигналов с периодами N, N/2 и т.д. мы никак не сможем установить только по сумме.

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

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

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