Осреднение сильно зашумлённых данных

Сейчас использую обычное срденее арифметическое

Среднее по серии осцилограмм, или скользящее среднее по 1 осциллограмме?

мощность лазера гуляет

Обычно ставят опорный датчик исходящей от лазера мощности и вычитают её флуктуации из сигнала с заранее определёнными калибровочными коэффициентами.

шумами/колебаниями в измеряемом объёме

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

Интересующая часть сигнала — центральный пик.

А поподробнее: форма, отношение с/ш, положение, ... ?
Для каждого варианта можно найти оптимальный фильтр, например ©:
фильтр Колмогорова-Винера выделяет форму сигнала из его смеси с помехой,
энергетический фильтр максимизирует отношение сигнал/помеха,
согласованный фильтр используют для вычисления корреляции и положения импульса на временнОй оси.

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

Да, есть 100500 хитростей, описанных в куче книг по ЦОС ©.

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

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

Среднее по серии осцилограмм, или скользящее среднее по 1 осциллограмме?

Среднее по серии. Т.е. что бы получить гладкий пик, усредняем по 64 осциллограммам.

Обычно ставят опорный датчик исходящей от лазера мощности и вычитают её флуктуации из сигнала с заранее определёнными калибровочными коэффициентами.

В данном случае, это просто ещё один источник ошибок. Так же шумит ФЭУ, цепи и пр. Но оно всё даёт некоторый постоянный фон. Тут есть проблема хуже ↴

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

У меня нехорошая ситуация — присутствует алиасинг частот. Лазер светит 8нс-импульсами с частотой 10Гц, отклик плазмы собираем с развёрткой 200нс (на таких временах все параметры исследуемой системы ~ заморожены). А собственные колебания в исследуемом объёме происходят в диапазоне частот 10..100кГц. Т.е. на деле мы захватываем сигнал в попадающий в случайные фазы колебаний + шум всей измерительной системы, отсюда получается большой разброс максимума интенсивности.

Соответственно, хочется обработать сигнал так, что бы он соответствовал среднему по колебаниям уровню. Для этого, вообще-то говоря, требуется что бы измерения были более-менее равномерно распределены по каждой фазе колебаний. Или при осреднении привязать сигнал к функции распределения максимумов интенсивности и использовать её как весовой коэффициент. Я попытался выяснить как она выглядит, и получается что распределение там далеко не равномерное и даже не симметричное: pix.academ.info/images/img/2017/08/29/7d442c887a979f61e33a03053602b7e3.png. Там конечно, маловато измерений для статистики, но тенденция на всех сигналах прослеживается одинаковая — на большей части осциллограмм записан слабый пик. Если получится, попытаюсь всё таки промерять форму этой функции.

А поподробнее: форма, отношение с/ш, положение, ... ?

Я снимаю 2х мерную плоскость (~100x100 точек), в каждой точке снимается серия осциллограмм. Так как измерений очень много, а эксперимент не должен проходить больше суток, то большую серию на каждую точку снять не получается. Максимум ~64 измерений. Форма и положение не очень важны. Собственно, самые искомые данные — это интенсивность сигнала в окрестностях вершины, необходимо что бы эта величина адекватно отражала среднюю интенсивность или хотя бы во всех точках съезжала на один и тот же процент от среднего. Требуемое отношение с/ш ~ 20dB.

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

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

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

«Да, уж...» :) Без привязки к фазе тяжко мерить... Кстати, стробировать измерительную систему (хотя бы ФЭУ) не пробовал?

самые искомые данные — это интенсивность сигнала в окрестностях вершины, необходимо что бы эта величина адекватно отражала среднюю интенсивность

Глянь MathWorks/ Signal Smoothing.

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

Идея правильная, можно модифицировать MathWorks/ Peak Analysis.

Но пока не соображу насколько это корректно.

Если процесс эргодичен на большом периоде, то корректно.

Фильтр Калмана.

1. Согласованная фильтрация. Возможно стоит подумать о нелинейной фильтрации.

2. Подучите теорию, а то вопрос совсем как-то несерьезно сформулирован.

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

Сильно смахивает что вам нужно автоматическое слежение реализованное контуром обратной связи с ПИД-регулятором.

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

Сильно смахивает что вам нужно автоматическое слежение реализованное контуром обратной связи с ПИД-регулятором.

Не понял идеи, на что обратную связь налагать предполагается? Слежение, к сожалению, не сделаешь, всё упирается в оборудование — я не могу запускать лазерную систему от внешнего сигнала, есть только возможность поменять частоту импульсов (в меньшую сторону).

Не понял идеи, на что обратную связь налагать предполагается?

Идея - убрать один из случайных факторов.

«автоматическое слежение реализованное контуром обратной связи с ПИД-регулятором.» - это один из классических способов убирания случайного параметра.

В целом непонятно что у тебя за система и какова твоя роль во всем этом. С моей колокольни это выглядит как попытка исправить глюки Андроида периодическими нажатиями на кнопку «Меню».

Выложи датасет.

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

Пусть у нас есть 400 реализаций случайной величины x(1:400).
Упорядочим данные y=sort(x).
Отбросим нижнюю и верхнюю кварты z=y(101:300).
Усредним m=mean(z).

Вот таким способом рекомендует усреднять атомиздат. Ссылку уже не помню.

Вот таким способом рекомендует усреднять атомиздат. Ссылку уже не помню.

Это называется линейные ранговые статистики ©.

Они удобны для «распределений с тяжёлыми хвостами»: Хьюбер П. Робастность в статистике ©.

А что делать, если я метод велфорда использую?

Тебе набросали много разных методов подходящие под конкретные ситуации, а навелосипедить «для красоты» можно много, начиная от фита сплайнами каждого сигнала, а потом сложения. И заканчивая «сглаживанием» по 3 точкам: на участках монотонности, каждые 3 точки 1,2,3, если |f(3) - f(1)| < |f(2) - f(1)|, то выкинуть или линеаризировать значение в 2.

Проснял прау точек по 2048 измерений и получил Максвелла (или что-то на него похожее): http://pix.academ.info/images/img/2017/08/31/b3a7effd63fdc3cb5e5d29f4408d2c4b... . Не очень понятно откуда он там вылезает, но что есть то есть. Буду думать что с этим делать.

Где датасет сырых данных этой картинки?

http://rgho.st/7bxBDnH4y, в npz пойдёт?

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import sys
import time

def main(args):

	for arg in args[1:]:
		if arg == "-plot":
			_plot = True; continue

		print("open \""+arg+"\"... ",end="")
		obj = np.load(arg)
		print("ok")

		fname = arg.replace(".npz",".csv")

		data = np.load(arg)["data"]
		(ncounts, nchs, npts) = data.shape

		out = np.mean(data,axis=0)
		err = np.std(data,axis=0)
		ms = np.amax(data,axis=2)


		pname = arg.replace(".npz",".png")
		print("save \""+pname+"\"... ",end="")

		fig = plt.figure()
		ax1 = fig.add_subplot(2,1,1)
		ax2 = fig.add_subplot(2,1,2)
		ax1.set_xlabel(r"time")
		ax1.set_ylabel(r"signal")
		ax2.set_xlabel(r"signal_max")
		ax2.set_ylabel(r"number")
		cs = [None,"k","r","g","b"]
		for i in range(ncounts):
			for j in range(1,data.shape[1]):
				ax1.plot(data[i,0],data[i,j],linewidth=0.025,color="k")

		for j in range(1,out.shape[0]):
			ax1.plot(out[0], out[j],color="r",linewidth=2)

		n, bins, patches = ax2.hist(ms[:,1], 25, normed=1, facecolor='grey', alpha=0.75)
		fig.savefig(pname,dpi=300)
		#~ plt.close()
		print("ok")

	return 0

if __name__ == '__main__':
	main(sys.argv)
	plt.show()

Интересующая часть сигнала — центральный пик.

что конкретно? время появления или ширина, или что-то другое?

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

попробуй wavelet transform. Только надо изначально выбрать шаблонную форму.

В целом непонятно что у тебя за система.

Писал же. Есть источник плазмы, на выхлоп светят импульсами лазера (8нс, <=10Гц) и смотрят отклик через оптическую систему с ФЭУ (ага, тоже свой шум вносит) так что бы построить 2д срез. Источник обратных связей не имеет и колеблется на своих внутренних частотах, связанных с выгоранием газа и over 9000 неустойчивостями. Соответственно, единственное что я могу менять — количество измерений на одну точку.

и какова твоя роль во всем этом

Собрать корректно-описывающие _средний_ отклик данные, которые потом будет обрабатывать леди-теоретик.

чо ты туда напаковал?

> library(RcppCNPy)
> mat <- npyLoad("mm.raw.npy")
Ошибка в npyLoad("mm.raw.npy") :Unsupported dimension in npyLoad

будь скромнее, запиши как «sets of arrays»

зы ну или как csv

csv: http://rgho.st/7Qcys6C9v

ты мог просто вывести массив, время идентично во всех датасетах

lsdir <- system("ls ./data", intern=T)
## в формате твоем лишнии пробелы ломающие парсер на character
data.raw <- lapply(lsdir, function(f) data.frame(lapply(read.csv2(paste0("./data/", f), as.is = TRUE), as.numeric)))
## ну и вот только один вариант времени замеров
> nrow(unique(t(sapply(data.raw, function(d) d[,1]))))
> [1] 1
## оставляем только данные 
> data.clean <- t(sapply(data.raw, function(d) d[,2]))
> str(data.clean)
> num [1:2048, 1:800] 0.000503 0.008543 0.008543 0.016583 0.088945 ...

сейчас дальше посмотрю

Вот фильтрация по pca разложению, путем разглядывания графика собственных чисел разложения мне почему то 15 первых оставить показалось :)

> reversepca <- function(nv, include, pca.object) {
+     sapply(nv,
+            function(i)(pca.object$x[,include] %*% pca.object$rotation[i,include])/
+                   (sum(pca.object$rotation[i,include]^2))^0.5)
+ }

> data.pca.reverse <- reversepca(1:800, c(1:15), prcomp(data.clean, center=F))
> plot(data.pca.reverse[1,1:800], type="l", ylim=range(unlist(data.pca.reverse)))
> for(i in 2:nrow(data.pca.reverse)) points(data.pca.reverse[i,1:800], type="l")
> points(colMeans(data.pca.reverse), type="l", col="red")

пичёк похоже двойной

http://pix.academ.info/img/2017/09/01/156941a02f1924e89d649cca7ae484f9.png

Хорошо, если я тебе посоветую поменять источник питания лазера, и фотоприемник?

Немного по существу:

на выхлоп светят импульсами лазера (8нс, <=10Гц)

как на тему светить импульсами код Баркера?

Соответственно, единственное что я могу менять — количество измерений на одну точку.

Из того что ты описал выше это не следует.

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