Detectando trenes de pulso en aumento

Question

Tengo un proceso de punto unidimensional que representa los tiempos de eventos y que también se mezcla con mucha información que considero ruido. El intervalo entre los eventos en el proceso de puntos es geométricamente creciente.

El número total de puntos está en el orden del millón. Me gustaría detectar (aproximadamente) procesos de esta forma en los datos. ¿Cuál es un buen método para hacer esto?

Como ejemplo, si mis puntos son 1, 7, 8, 10, 16, 25, 28, 30, 31, 50, 52 me gustaría detectar 7, 10, 16, 28, 52. Dado que mis valores son números reales y provienen de mediciones, necesito que se toleren pequeños errores en los datos.

Como otro ejemplo, los puntos en tal tren de impulsos estarán en la forma $t_0 + \sum_{i=1}^n a*k^{i-1}$. En el ejemplo anterior, $t_0 = 7$, $a = 3$, $k = 2.

Answer 1

Puedes intentar buscar cuartetos $(a, b, c, d)$ donde $a$ y $d$ no estén demasiado separados, y $(a-b)(c-d)-(b-c)^2$ sea relativamente pequeño. Luego ve si el cuarteto se extiende a un quinto, y así sucesivamente. ¿Te perderás alguna vez el término ocasional del conjunto de datos?

Answer 2

Esta respuesta está incompleta, pero es donde comenzaría si me asignaran este problema. Sabemos la señal que estamos buscando:

$$ x(t) = \delta( t - t_0 ) \ast \sum_{n=0}^{\lfloor \log(N)/\log(p) \rfloor - 1} \delta \left( t - (p^{n}-1)k \right) $$

donde $\ast$ es la convolución, $N$ es el número de muestras, $p,k\in\mathbb{Z}_{+}$, y $t_0$ es el inicio del tren de impulsos geométricos. Para el ejemplo en el OP, estos valores son $t_0 = 7$, $k=3$, y $p=2$, lo que da lugares de impulso en $7+\{0,3,9,21,45,\dots\}$.

El término $\delta(t-t_0)$ en lo anterior representa un desplazamiento desconocido en la secuencia, así que lo ignoraremos por ahora. La transformada de Fourier del resto de lo anterior es fácil de encontrar ya que es solo la suma de exponenciales complejas:

$$ X(\omega,p,k) = \sum_{n=0}^{\lfloor \log(N)/\log(p) \rfloor - 1} e^{-j\omega (p^{n}-1)k } $$

Esto nos permite crear fácilmente un banco de filtros coincidentes. Por supuesto, esto significa que debes establecer límites en el dominio de la señal, pero dado que es de naturaleza geométrica, no hay tantos valores razonables (por ejemplo, si $p=15$, entonces solo habrá $5$ valores distintos de cero de $10^6$ y es probable que la señal esté por debajo del piso de ruido).

Entonces, para un par dado de $p,k$, genera la siguiente secuencia:

$$ z(t) = \operatorname{FFT}^{-1} \left( \operatorname{FFT}(x) \cdot X^{\ast}(\omega,p,k) \right) $$

Si $p,k$ coinciden con los valores en la señal, entonces $z(t)$ tendrá un pico en $t_0$.

Creo que la parte complicada va a ser la detección, pero tal vez no, ya que deberías poder averiguar la amplitud esperada del filtro coincidente ya que cada impulso tiene amplitud unitaria. Sin embargo, varias secuencias de entrada parecerán casi idénticas y en presencia de ruido puede ser difícil distinguirlas. Por ejemplo, el ejemplo que diste de $t_0=7$, $p=2$, y $k=3$ es casi lo mismo que $t_0=10$, $p=2$, y $k=6$.

A continuación se muestra un código de Matlab que podría ayudarte a comenzar. Como dije, esta respuesta está incompleta, así que tenlo en cuenta.

N = 1e6;
x = zeros(N,1);

% Los trenes de impulsos geométricos que estamos tratando de detectar.
% Damos a cada uno un punto de inicio arbitrario.
x( 63 + ( 5.^( 0 : log(N) / log(5) - 1 ) - 1 ) * 2 ) = 1;
x( 7 + ( 2.^( 0 : log(N) / log(2) - 1 ) - 1 ) * 3 ) = 1;

% Agrega algo de ruido a los datos.
randinds = randperm(N);      % Permutación aleatoria para el ruido.
x(randinds(1:10e4)) = 1;     % Agrega algunos puntos de ruido.

Nfft = 2*N;
X = fftshift( fft( x, Nfft ) );    % DFT de los datos de entrada.
f = (-0.5:1/Nfft:0.5-1/Nfft)';     % Eje de frecuencia normalizado.

% Define un banco de filtros coincidentes.
p = ( 2 : 10 );
k = ( 1 : 10 );
Nmf = length( p );
K = length( k );
mf = complex( zeros( N, K, Nmf ) );

% Ahora creamos y aplicamos cada filtro coincidente.
for ii = 1 : Nmf
  for jj = 1 : K

    % Crea el filtro coincidente para este factor geométrico.
    Q = 0;
    for n = 0 : floor( log(N) / log(p(ii)) ) - 1
      Q = Q + exp( -1j*2*pi*f * ( p(ii)^n - 1 ) * k(jj) );
    end

    % Ahora aplica el filtro coincidente y almacena el resultado.
    tmp = ifft( ifftshift( X .* conj(Q/norm(Q)) ) );
    mf(:,jj,ii) = tmp(1:N);

  end

end

% Crea la estadística de prueba tomando el valor máximo de cada
% par (k,p).
T = squeeze( max( abs( mf ), [], 1 ) );