¿Qué es un algoritmo para volver a muestrear de una tasa variable a una tasa fija?

27

Tengo un sensor que informa sus lecturas con una marca de tiempo y un valor. Sin embargo, no genera lecturas a una velocidad fija.

Encuentro los datos de tasa variable difíciles de manejar. La mayoría de los filtros esperan una frecuencia de muestreo fija. Dibujar gráficos también es más fácil con una frecuencia de muestreo fija.

¿Existe un algoritmo para volver a muestrear de una frecuencia de muestreo variable a una frecuencia de muestreo fija?

FigBug
fuente
Esta es una publicación cruzada de los programadores. Me dijeron que este es un mejor lugar para preguntar. programmers.stackexchange.com/questions/193795/…
FigBug
¿Qué determina cuándo el sensor informará una lectura? ¿Envía una lectura solo cuando la lectura cambia? Un enfoque simple sería elegir un "intervalo de muestra virtual" (T) que sea más pequeño que el tiempo más corto entre lecturas generadas. En la entrada del algoritmo, almacene solo la última lectura informada (CurrentReading). En la salida del algoritmo, informe el CurrentReading como una "nueva muestra" cada T segundos para que el servicio de filtro o gráfico reciba lecturas a una velocidad constante (cada T segundos). Sin embargo, no tengo idea si esto es adecuado en su caso.
user2718
Intenta muestrear cada 5 ms o 10 ms. Pero es una tarea de baja prioridad, por lo que puede perderse o retrasarse. Tengo el tiempo exacto a 1 ms. El procesamiento se realiza en la PC, no en tiempo real, por lo que un algoritmo lento está bien si es más fácil de implementar.
FigBug
1
¿Has echado un vistazo a una reconstrucción de Fourier? Hay una transformación de Fourier basada en datos muestreados de manera desigual. El método habitual es transformar una imagen de Fourier en un dominio de tiempo muestreado uniformemente.
mbaitoff
3
¿Conoces alguna característica de la señal subyacente que estás muestreando? Si los datos con espacios irregulares todavía se encuentran a una frecuencia de muestreo razonablemente alta en comparación con el ancho de banda de la señal que se está midiendo, entonces algo simple como la interpolación polinómica a una cuadrícula de tiempo uniformemente espaciada podría funcionar bien.
Jason R

Respuestas:

21

El enfoque más simple es hacer algún tipo de interpolación de spline como sugiere Jim Clay (lineal o de otro tipo). Sin embargo, si tiene el lujo del procesamiento por lotes, y especialmente si tiene un conjunto sobredeterminado de muestras no uniformes, hay un algoritmo de "reconstrucción perfecta" que es extremadamente elegante. Por razones numéricas, puede no ser práctico en todos los casos, pero al menos vale la pena conocerlo conceptualmente. Lo leí por primera vez en este artículo .

El truco consiste en considerar que su conjunto de muestras no uniformes ya ha sido reconstruido a partir de muestras uniformes mediante interpolación sinc . Siguiendo la notación en el documento:

y(t)=k=1Ny(kT)sin(π(tkT)/T)π(tkT)/T=k=1Ny(kT)sinc(tkTT).

Tenga en cuenta que esto proporciona un conjunto de ecuaciones lineales, una para cada muestra no uniforme , donde las incógnitas son las muestras igualmente espaciadas , así:y ( k T )y(t)y(kT)

[y(t0)y(t1)y(tm)]=[sinc(t0TT)sinc(t02TT)sinc(t0nTT)sinc(t1TT)sinc(t12TT)sinc(t1nTT)sinc(tmTT)sinc(tm2TT)sinc(tmnTT)][y(T)y(2T)y(nT)].

En la ecuación anterior, es el número de muestras uniformes desconocidas, es el inverso de la frecuencia de muestreo uniforme, es el número de muestras no uniformes (que pueden ser mayores que ). Al calcular la solución de mínimos cuadrados de ese sistema, se pueden reconstruir las muestras uniformes. Técnicamente, solo se necesitan muestras no uniformes, pero dependiendo de cuán "dispersas" estén en el tiempo, la matriz de interpolación puede estar terriblemente mal acondicionada . Cuando ese es el caso, usar más muestras no uniformes generalmente ayuda.T mnTmnnn

Como ejemplo de juguete, aquí hay una comparación (usando numpy ) entre el método anterior y la interpolación de spline cúbico en una cuadrícula ligeramente nerviosa:

Reconstrucción Sinc vs Spline cúbico de muestras no uniformes

(El código para reproducir la trama anterior se incluye al final de esta respuesta)

Dicho todo esto, para métodos robustos de alta calidad, comenzar con algo en uno de los siguientes documentos probablemente sería más apropiado:

A. Aldroubi y Karlheinz Grochenig, Muestreo no uniforme y reconstrucción en espacios invariables por turnos , SIAM Rev., 2001, no. 4, 585-620. ( enlace pdf ).

K. Grochenig y H. Schwab, Métodos de reconstrucción local rápida para muestreo no uniforme en espacios invariables por desplazamiento , SIAM J. Matrix Anal. Appl., 24 (2003), 899-913.

-

import numpy as np
import pylab as py

import scipy.interpolate as spi
import numpy.random as npr
import numpy.linalg as npl

npr.seed(0)

class Signal(object):

    def __init__(self, x, y):
        self.x = x
        self.y = y

    def plot(self, title):
        self._plot(title)
        py.plot(self.x, self.y ,'bo-')
        py.ylim([-1.8,1.8])
        py.plot(hires.x,hires.y, 'k-', alpha=.5)

    def _plot(self, title):
        py.grid()
        py.title(title)
        py.xlim([0.0,1.0])

    def sinc_resample(self, xnew):
        m,n = (len(self.x), len(xnew))
        T = 1./n
        A = np.zeros((m,n))

        for i in range(0,m):
            A[i,:] = np.sinc((self.x[i] - xnew)/T)

        return Signal(xnew, npl.lstsq(A,self.y)[0])

    def spline_resample(self, xnew):
        s = spi.splrep(self.x, self.y)
        return Signal(xnew, spi.splev(xnew, s))

class Error(Signal):

    def __init__(self, a, b):
        self.x = a.x
        self.y = np.abs(a.y - b.y)

    def plot(self, title):
        self._plot(title)
        py.plot(self.x, self.y, 'bo-')
        py.ylim([0.0,.5])

def grid(n): return np.linspace(0.0,1.0,n)
def sample(f, x): return Signal(x, f(x))

def random_offsets(n, amt=.5):
    return (amt/n) * (npr.random(n) - .5)

def jittered_grid(n, amt=.5):
    return np.sort(grid(n) + random_offsets(n,amt))

def f(x):
    t = np.pi * 2.0 * x
    return np.sin(t) + .5 * np.sin(14.0*t)

n = 30
m = n + 1

# Signals
even   = sample(f, np.r_[1:n+1] / float(n))
uneven = sample(f, jittered_grid(m))
hires  = sample(f, grid(10*n))

sinc   = uneven.sinc_resample(even.x)
spline = uneven.spline_resample(even.x)

sinc_err   = Error(sinc, even)
spline_err = Error(spline, even)

# Plot Labels
sn = lambda x,n: "%sly Sampled (%s points)" % (x,n)
r  = lambda x: "%s Reconstruction" % x
re = lambda x: "%s Error" % r(x)

plots = [
    [even,       sn("Even", n)],
    [uneven,     sn("Uneven", m)],
    [sinc,       r("Sinc")],
    [sinc_err,   re("Sinc")],
    [spline,     r("Cubic Spline")],
    [spline_err, re("Cubic Spline")]
]

for i in range(0,len(plots)):
    py.subplot(3, 2, i+1)
    p = plots[i]
    p[0].plot(p[1])

py.show()
Datageist
fuente
Buen método y código. Pero para con algunos abandonos (por ejemplo, [0 1 - 3 4 5 - 7 8] T), que creo que es la pregunta de los OP, ¿no son todos los sinc s en la matriz todos 0? Claro que hay formas de arreglar eso, pero. tj=jT
denis
Según tengo entendido, la pregunta del OP es sobre muestras que se descartan y / o se retrasan. Este método es básicamente un sistema de ecuaciones sobredeterminado, por lo que las muestras descartadas solo aparecen como incógnitas (no como puntos de datos con un valor de 0). ¿O tal vez eso no es lo que estás preguntando?
Datageist
¿Qué sucede si son todos enteros (T = 1)? Digamos que tenemos puntos de datos [ ] para , un conjunto de enteros distintos de cero, por ejemplo, {-1 1} o {-2 -1 1 2}; ¿No es el interpolado , independientemente de , o me he perdido algo? tj j J y 0 = 0 y jj,yjjJy0=0yj
Denis
Si las frecuencias de muestreo son exactamente idénticas (con puntos faltantes), la matriz de interpolación será escasa (porque cada salida depende de una sola entrada). En general, la tasa de muestra promedio de las muestras no uniformes debe ser mayor que la tasa de reconstrucción uniforme. En otras palabras, necesitaría reconstruir a una velocidad menor para "llenar los vacíos" (T> 1 para su ejemplo). Veo su punto, sin embargo.
Datageist
2
Respuestas como esta son de oro puro.
Ahmed Fasih
6

Esto suena como un problema de conversión asincrónica de frecuencia de muestreo. Para convertir de una frecuencia de muestreo a otra, podemos calcular la representación de tiempo continuo de la señal realizando una interpolación sinc y luego volver a muestrear a nuestra nueva frecuencia de muestreo. Lo que estás haciendo no es muy diferente. Debe volver a muestrear su señal para tener tiempos de muestra fijos.

La señal de tiempo continuo se puede calcular convolucionando cada muestra con una función sinc. Dado que la función sinc continúa indefinidamente, utilizamos algo más práctico como un sinc con ventana de una longitud finita práctica. La parte difícil es que debido a que sus muestras se mueven a tiempo, puede ser necesario utilizar un sinc con un desfase de fase diferente para cada muestra al volver a muestrear.

Señal de tiempo continuo de la señal muestreada:

x(t)=n=x[n]sinc(tnTsTs)

donde es su tiempo de muestra. En su caso, sin embargo, su tiempo de muestreo no es fijo. Así que creo que debe reemplazarlo con el tiempo de muestra en esa muestra.Ts

x(t)=n=x[n]sinc(tnTs[n]Ts[n])

A partir de esto, puede volver a muestrear la señal:

y[n]=x(nTns )

donde es el tiempo de muestra deseado.Tns

Poniendo todo junto, obtienes:

y[m]=n=x[n]sinc(mTnsnTs[n]Ts[n])

Como esto no es causal o manejable, la función sinc se puede reemplazar con una función de soporte finito y los límites de suma se ajustan en consecuencia.

Deje que kernel (t) sea una ventana sinc u otra función similar de longitud 2k y luego:

y[m]=n=kkx[n]kernel(mTnsnTs[n]Ts[n])

Espero que esto ayude ..., pero podría haber cometido un error en el camino y podría ser un poco intensivo en matemáticas. Recomendaría investigar la conversión de frecuencia de muestreo para obtener más información. Quizás alguien más aquí podría dar una mejor explicación o solución también.

Jacob
fuente
El uso de una función sinc para reconstruir una versión continua de una señal requiere que las muestras estén igualmente espaciadas, por lo que la función sinc tendrá que adaptarse al espaciado de la muestra arbitraria. Podría ser bastante difícil de implementar.
usuario2718
Sí, esto no sería muy eficiente para hacer exactamente como se ve aquí. Requeriría calcular nuevos coeficientes del núcleo para cada tiempo de muestra diferente. Sin embargo, se podría calcular una colección de varios núcleos y cuantificar el tiempo en uno de estos. Habría un impacto en el rendimiento en relación con el error de cuantificación.
Jacob
También puede calcular previamente una sola tabla de búsqueda sinc e interpolar entre los puntos de esa tabla de búsqueda.
jms
5

Creo que la respuesta de Jacob es muy viable.

Un método más fácil que probablemente no sea tan bueno en términos de introducción de distorsión es hacer una interpolación polinómica. Usaría interpolación lineal (fácil, no tan buena señal de rendimiento) o splines cúbicos (todavía no demasiado duro, mejor rendimiento de señal) para producir muestras en cualquier momento que desee de sus muestras de tiempo arbitrario.

Jim Clay
fuente
1
Su respuesta parece mucho más fácil de implementar que la de Jacob, así que lo hice primero. Parece estar funcionando, pero aún no he realizado un conjunto completo de pruebas.
FigBug
1
@FigBug: si tienes tiempo, agrega un comentario con tu solución final.
user2718
2


Nnear

Nnear
y1y1
y0(y1+y1)/2
[xi,yi]xi=0
y0yi
[xi,yi]
a+bxy0a
Uno puede ajustar cuadráticos, cúbicos, seno-cosenos ... de la misma manera.



2πftf


La interpretación lineal con 4 u 6 u 8 vecinos podría funcionar lo suficientemente bien para sus datos.
Sugeriría comenzar con un método que comprenda a fondo antes de sumergirse en splines, como sinc ... aunque también puede ser divertido.


Otro método bastante diferente es la ponderación de distancia inversa . Es fácil de implementar (ver idw-interpolation-with-python en SO), también funciona en 2d 3d y superior, pero es difícil de analizar teóricamente.

(Obviamente, NINGÚN método de interpolación único puede adaptarse a los miles de millones de combinaciones de
[señal, ruido, métrica de error, función de prueba] que ocurren en la realidad.
Hay más métodos en el mundo, con más botones, que funciones de prueba.
Sin embargo, una galería de métodos y funciones de prueba pueden ser útiles).

denis
fuente
1

Si trabaja con matlab, puede hacerlo trabajando con series de tiempo.

time  % is your starting vector of time

data % vector of data you want to resample 

data_TS = timeseries(data,time); % define the data as a timeseries 

new_time = time(0):dt:time(end); % new vector of time with fixed dt=1/fs

data_res = resample(data_TS,new_time); % data resampled at constant fs
Rhei
fuente
0

Antes de comenzar a realizar un procesamiento exótico, puede probar algo simple como esto (pseudocódigo, sin interpolación, pero podría agregarse)

TimeStamp[]  //Array of Sensor TimeStamps -NULL terminated – TimeStamp[i] corresponds to Reading[i]
Reading[]      //Array of Sensor Readings       -NULL terminated

AlgorithmOut   //Delimited file of of readings in fixed sample time (5ms) 
CurrentSavedReading = Reading[0]

SampleTime=TimeStamp[0] //ms virtual sample time, 5ms fixed samples

i = 0 // loop index
While(TimeStamp[i] != NULL)
{
   FileWrite (CurrentSavedReading, AlgorithmOut)//write value to file
   SampleTime = SampleTime + 5//ms
   if(SampleTime > TimeStamp[i])
   {
      i++
      CurrentSavedReading = Reading[i]
   }
}
usuario2718
fuente
0

La respuesta de IMHO Datageist es correcta, la respuesta de Jacob no lo es. Una manera fácil de verificar esto es que el algoritmo sugerido por datageist está garantizado para interpolar a través de las muestras originales (suponiendo una precisión numérica infinita), mientras que la respuesta de Jacob no.

  • Para el caso de muestreo uniforme, el conjunto de funciones sinc es ortogonal: si cada función sinc desplazada se discretiza sobre las muestras de entrada, forman una matriz de identidad infinita. Esto se debe a que sin (n pi) / (n pi) es cero para todos los n excepto n = 0.
  • Sin embargo, esta propiedad no puede simplemente extrapolarse al caso no uniforme: un conjunto similar de funciones sinc, discretizadas sobre las muestras de entrada, producirá una matriz no trivial. Por lo tanto, las contribuciones de las muestras circundantes no serán cero y la reconstrucción ya no se interpolará a través de las muestras de entrada.
pvv
fuente