O que é um algoritmo para amostrar novamente de uma taxa variável para uma taxa fixa?

27

Eu tenho um sensor que relata suas leituras com um carimbo de hora e um valor. No entanto, ele não gera leituras a uma taxa fixa.

Acho os dados de taxa variável difíceis de lidar. A maioria dos filtros espera uma taxa de amostragem fixa. Desenhar gráficos é mais fácil com uma taxa de amostragem fixa também.

Existe um algoritmo para reamostrar de uma taxa de amostra variável para uma taxa de amostra fixa?

FigBug
fonte
Esta é uma postagem cruzada dos programadores. Foi-me dito que este é um lugar melhor para perguntar. programmers.stackexchange.com/questions/193795/...
FigBug
O que determina quando o sensor reportará uma leitura? Ele envia uma leitura somente quando a leitura muda? Uma abordagem simples seria escolher um "intervalo de amostra virtual" (T) que seja apenas menor que o menor tempo entre as leituras geradas. Na entrada do algoritmo, armazene apenas a última leitura relatada (CurrentReading). Na saída do algoritmo, relate o CurrentReading como uma “nova amostra” a cada T segundos, para que o serviço de filtro ou gráfico receba leituras a uma taxa constante (a cada T segundos). Não faço ideia se isso é adequado no seu caso.
precisa saber é o seguinte
Ele tenta amostrar a cada 5ms ou 10ms. Mas é uma tarefa de baixa prioridade, por isso pode ser perdida ou atrasada. Eu tenho o tempo exato para 1 ms. O processamento é feito no PC, não em tempo real; portanto, um algoritmo lento é aceitável se for mais fácil de implementar.
precisa
11
Você deu uma olhada em uma reconstrução de Fourier? Há uma transformação de Fourier baseada em dados amostrados de maneira desigual. A abordagem mais comum é transformar uma imagem de fourier de volta ao domínio do tempo de amostragem uniforme.
precisa saber é o seguinte
3
Você conhece alguma característica do sinal subjacente que está sendo amostrado? Se os dados com espaçamento irregular ainda estiverem com uma taxa de amostragem razoavelmente alta em comparação com a largura de banda do sinal que está sendo medido, algo simples como interpolação polinomial para uma grade de tempo com espaçamento uniforme pode funcionar bem.
Jason R

Respostas:

21

A abordagem mais simples é fazer algum tipo de interpolação de splines, como Jim Clay sugere (linear ou não). No entanto, se você tem o luxo do processamento em lote e, especialmente, se você tem um conjunto sobredeterminado de amostras não uniformes, existe um algoritmo de "reconstrução perfeita" extremamente elegante. Por razões numéricas, pode não ser prático em todos os casos, mas vale a pena conhecer conceitualmente. Eu li pela primeira vez sobre isso neste artigo .

O truque é considerar o seu conjunto de amostras não uniformes como já tendo sido reconstruído a partir de amostras uniformes por meio de interpolação sincera . Após a notação no artigo:

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

Observe que isso fornece um conjunto de equações lineares, uma para cada amostra não uniforme , onde as incógnitas são as amostras igualmente espaçadas , da seguinte forma: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)].

Na equação acima, é o número de amostras uniformes desconhecidas, é o inverso da taxa de amostragem uniforme, e é o número de amostras não uniformes (que podem ser maiores que ). Calculando a solução de mínimos quadrados desse sistema, as amostras uniformes podem ser reconstruídas. Tecnicamente, apenas amostras não uniformes são necessárias, mas dependendo de quão "dispersas" elas estão no tempo, a matriz de interpolação pode estar terrivelmente mal condicionada . Nesse caso, o uso de mais amostras não uniformes geralmente ajuda.T m n nnTmnn

Como um exemplo de brinquedo, veja uma comparação (usando numpy ) entre o método acima e a interpolação de spline cúbico em uma grade levemente instável:

Reconstrução Sinc vs Cubine Spline de amostras não uniformes

(O código para reproduzir o gráfico acima está incluído no final desta resposta)

Tudo o que foi dito, para métodos robustos de alta qualidade, começando com algo em um dos seguintes documentos provavelmente seria mais apropriado:

A. Aldroubi e Karlheinz Grochenig, Amostragem não uniforme e reconstrução em espaços invariantes a turnos , SIAM Rev., 2001, no. 4, 585-620. ( link em pdf ).

K. Grochenig e H. Schwab, Métodos rápidos de reconstrução local para amostragem não uniforme em espaços invariantes por turnos , 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()
datagrama
fonte
Bom método e código. Mas para com algumas desistências (por exemplo, [0 1 - 3 4 5 - 7 8] T), que eu acho que é a questão dos OPs, os sinc s na matriz não são todos 0? Claro, existem maneiras de consertar isso, mas. tj=jT
Denis
Pelo que entendi, a pergunta do OP é sobre amostras que são descartadas e / ou atrasadas. Esse método é basicamente apenas um sistema superdeterminado de equações; portanto, as amostras descartadas aparecem apenas como desconhecidas (não como pontos de dados com o valor 0). Ou talvez não seja isso que você está perguntando?
datageist
O que acontece se são todos números inteiros (T = 1)? Digamos que temos pontos de dados [ ] para , um conjunto de números inteiros diferentes de zero, por exemplo, {-1 1} ou {-2 -1 1 2}; não é o interpolado , independentemente do - ou eu perdi alguma coisa? j , y jtjj,yjy 0 = 0 y jjJy0=0yj
Denis
Se as taxas de amostragem forem exatamente idênticas (com pontos ausentes), a matriz de interpolação será esparsa (porque cada saída depende de apenas uma entrada). Em geral, a taxa média de amostra das amostras não uniformes precisa ser maior que a taxa de reconstrução uniforme. Em outras palavras, você precisará reconstruir em uma taxa mais baixa para "preencher as lacunas" (T> 1 no seu exemplo). Entendo o seu ponto.
datageist
2
Respostas como esta são ouro puro.
Ahmed Fasih,
6

Isso parece um problema de conversão de taxa de amostragem assíncrona. Para converter de uma taxa de amostragem para outra, podemos calcular a representação em tempo contínuo do sinal executando uma interpolação sinc e, em seguida, reamostrar em nossa nova taxa de amostragem. O que você está fazendo não é muito diferente. Você precisa reamostrar seu sinal para ter tempos de amostra fixados.

O sinal de tempo contínuo pode ser calculado convocando cada amostra com uma função sinc. Como a função sinc continua indefinidamente, usamos algo mais prático como um sinc com janelas de um comprimento finito prático. A parte complicada é que, como suas amostras se movem com o tempo, um sinc com um deslocamento de fase diferente pode precisar ser usado para cada amostra ao reamostrar.

Sinal de tempo contínuo do sinal amostrado:

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

onde é o seu tempo de amostra. No seu caso, no entanto, seu tempo de amostra não é fixo. Então, acho que você precisa substituí-lo pelo tempo da amostra nessa amostra.Ts

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

A partir disso, você pode reamostrar o sinal:

y[n]=x(nTns )

onde é o tempo de amostra desejado.Tns

Juntando tudo, você obtém:

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

Como isso não é causal ou tratável, a função sinc pode ser substituída por uma função de suporte finito e os limites de soma ajustados de acordo.

Seja o kernel (t) uma função sinc de janela ou outra função similar de comprimento 2k, então:

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

Espero que isso ajude ..., mas posso ter cometido um erro ao longo do caminho e pode ser um pouco intensivo em matemática. Eu recomendaria pesquisar a conversão da taxa de amostra para obter mais informações. Talvez alguém aqui também possa dar uma explicação ou solução melhor.

Jacob
fonte
O uso de uma função sinc para reconstruir uma versão contínua de um sinal requer que as amostras sejam igualmente espaçadas, de modo que a função sinc terá que se adaptar ao espaçamento arbitrário da amostra. Pode ser bastante difícil de implementar.
usar o seguinte comando
Sim, isso não seria muito eficiente para fazer exatamente como visto aqui. Isso exigiria o cálculo de novos coeficientes do kernel para cada tempo de amostra diferente. Uma coleção de vários núcleos pode ser computada, no entanto, e o tempo quantificado para um deles. Haveria um impacto no desempenho em relação ao erro de quantização.
Jacob
Você também pode pré-calcular uma única tabela de pesquisa sinc e interpolar entre os pontos dessa tabela de pesquisa.
JMS
5

Eu acho que a resposta de Jacob é muito viável.

Um método mais fácil que provavelmente não é tão bom em termos de introdução de distorção é fazer a interpolação polinomial. Eu usaria interpolação linear (fácil, não como um bom sinal em termos de desempenho) ou splines cúbicos (ainda não muito difíceis, melhor desempenho de sinais) para produzir amostras a qualquer momento que você quiser de suas amostras de tempo arbitrárias.

Jim Clay
fonte
11
Você responde parece muito mais fácil de implementar do que o de Jacob, então eu o segui primeiro. Parece estar funcionando, mas ainda não executei um conjunto completo de testes.
FigBug
11
@FigBug -Se você tiver tempo, adicione um comentário com sua solução final.
usar o seguinte comando
2

(Um mês depois), existem duas opções principais para qualquer método de interpolação:
1) o número de pontos de dados mais próximo do ponto ausente a ser usado; 2 4 6 ... 2) a classe de funções a serem usadas: linear, polinomial, seno-cosseno (Fourier), cúbico por partes (spline B ou spline de interpolação), tipo sinc ... (A opção 0 é se você deve usar o método e o código de outra pessoa ou fazer você mesmo.)Nnear

Ajustar uma linha reta aos pontos é fácil: 2 pontos [-1, ], [1, ]: estimam pontos com média : média geral : consulte, por exemplo, Receitas Numéricas p. 781: ajuste uma linha e . Pode-se ajustar quadráticos, cúbicos, seno-cossenos ... da mesma maneira.Nnear
y1y1
[ x i , y i ] x i = 0y0(y1+y1)/2
[xi,yi]xi=0
y i [ x i , y i ]y0yi
[xi,yi]
y 0aa+bxy0a

Entendo que você tem dados espaçados uniformemente com alguns pontos faltando, não é mesmo?
Como funciona a interpolação linear neste caso?
Bem, vamos tentar cos com = 0.25: 1 0 -1 0 1 0 -1 0 ... 2 vizinhos de qualquer ponto médio a 0, terrível. 4 vizinhos: média de [1 0 (faltando -1) 0 1] = 1/2, terrível. (Experimente o filtro de 4 vizinhos [-1 3 3 -1] / 4 sobre isso.)f2πftf


A interação linear com 4, 6 ou 8 vizinhos pode funcionar bem o suficiente para seus dados.
Eu sugiro começar com um método que você entenda bem antes de mergulhar em splines, tipo sinc ... embora esses também possam ser divertidos.


Outro método bastante diferente é a ponderação de distância inversa . É fácil de implementar (consulte idw-interpolation-with-python no SO), funciona em 2d 3d e superior, mas é muito difícil de analisar teoricamente.

(Obviamente, NÃO método de interpolação único pode, eventualmente, ajustar os zilhões de combinações de
[sinal, ruído, métrica de erro, a função de teste] que ocorrem na realidade.
Há mais métodos do mundo, com mais botões, de funções de teste.
No entanto uma galeria métodos e funções de teste podem ser úteis.)

denis
fonte
1

Se você trabalha com o matlab, pode fazê-lo trabalhando com séries temporais.

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
fonte
0

Antes de iniciar um processamento exótico, você pode tentar algo simples como este (pseudo-código - sem interpolação, mas isso pode ser adicionado)

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]
   }
}
user2718
fonte
0

A resposta do IMHO Datageist está correta, a resposta de Jacob não está. Uma maneira fácil de verificar isso é que é garantido que o algoritmo sugerido pelo datagrafista interpole pelas amostras originais (assumindo precisão numérica infinita), enquanto a resposta de Jacob não.

  • Para o caso de amostragem uniforme, o conjunto de funções sinc é ortogonal: se cada função sinc deslocada for discretizada sobre as amostras de entrada, elas formarão uma matriz de identidade infinita. Isso ocorre porque sin (n pi) / (n pi) é zero para todos os n, exceto n = 0.
  • No entanto, essa propriedade não pode simplesmente ser extrapolada para o caso não uniforme: um conjunto semelhante de funções sinc, discretizadas sobre as amostras de entrada, produzirá uma matriz não trivial. Portanto, as contribuições das amostras circundantes não serão nulas e a reconstrução não interpolará mais através das amostras de entrada.
pvv
fonte