Co to jest algorytm ponownego próbkowania ze zmiennej stawki do stałej dawki?

Mam czujnik, który raportuje swoje odczyty ze znacznikiem czasu i wartością. Jednak nie generuje odczytów ze stałą szybkością.

Trudno mi poradzić sobie z danymi o zmiennej stopie procentowej. Większość filtrów oczekuje stałej częstotliwości próbkowania. Rysowanie wykresów jest również łatwiejsze przy stałej częstotliwości próbkowania.

Czy istnieje algorytm do ponownego próbkowania ze zmienną częstotliwością próbkowania na stałą częstotliwość próbkowania?

Najprostszym podejściem jest wykonanie interpolacji spline, jak sugeruje Jim Clay (liniowy lub inny). Jeśli jednak masz luksus przetwarzania wsadowego, a zwłaszcza jeśli masz zbyt określony zestaw niejednorodnych próbek, istnieje niezwykle elegancki algorytm „idealnej rekonstrukcji”. Z powodów numerycznych może nie być praktyczne we wszystkich przypadkach, ale przynajmniej warto wiedzieć o tym koncepcyjnie. Najpierw przeczytałem o tym w tym artykule .

Sztuką jest rozważenie zestawu niejednorodnych próbek jako już zrekonstruowanych z jednorodnych próbek poprzez interpolację sinc . Zgodnie z notacją w artykule:

$$y(t) = \sum_{k=1}^{N}{y(kT)\frac{\sin(\pi(t - kT)/T)}{\pi(t - kT)/T}} = \sum_{k=1}^{N}{y(kT)\mathrm{sinc}(\frac{t - kT}{T})}.$$

Zauważ, że zapewnia to zestaw równań liniowych, po jednym dla każdej niejednorodnej próbki , gdzie nieznane są próbki o równych odstępach , tak:y ( k T $y(t)$$y(kT)$

$$\begin{bmatrix} y(t_0) \\ y(t_1) \\ \cdots \\ y(t_m) \end{bmatrix} = \begin{bmatrix} \mathrm{sinc}(\frac{t_0 - T}{T}) & \mathrm{sinc}(\frac{t_0 - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_0 - nT}{T}) \\ \mathrm{sinc}(\frac{t_1 - T}{T}) & \mathrm{sinc}(\frac{t_1 - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_1 - nT}{T}) \\ \cdots & \cdots & \cdots &\cdots \\ \mathrm{sinc}(\frac{t_m - T}{T}) & \mathrm{sinc}(\frac{t_m - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_m - nT}{T}) \end{bmatrix} \begin{bmatrix} y(T) \\ y(2T) \\ \cdots \\ y(nT) \end{bmatrix}.$$

W powyższym równaniu jest liczbą nieznanych jednorodnych próbek, jest odwrotnością częstości jednolitych próbek, a jest liczbą niejednolitych próbek (która może być większa niż ). Obliczając najmniejszych kwadratów rozwiązania tego systemu, jednorodne próbki mogą być odtworzone. Technicznie konieczne jest tylko próbek niejednorodnych, ale w zależności od tego, jak „rozproszone” są w czasie, macierz interpolacji może być strasznie źle uwarunkowana . W takim przypadku zwykle pomaga użycie większej liczby niejednolitych próbek.T m n $n$$T$$m$$n$$n$

Jako zabawkowy przykład, oto porównanie (przy użyciu numpy ) między powyższą metodą a interpolacją splajnu sześciennego na lekko zniekształconej siatce:

(Kod do odtworzenia powyższej fabuły znajduje się na końcu tej odpowiedzi)

Biorąc to wszystko pod uwagę, w przypadku wysokiej jakości, solidnych metod, rozpoczęcie od czegoś w jednym z następujących dokumentów prawdopodobnie byłoby bardziej odpowiednie:

A. Aldroubi i Karlheinz Grochenig, Pobieranie próbek i rekonstrukcja Nonuniform w przestrzeniach niezmiennych , SIAM Rev., 2001, nr. 4, 585–620. ( link pdf ).

K. Grochenig i H. Schwab, Szybkie lokalne metody rekonstrukcji dla niejednorodnego próbkowania w przestrzeniach niezmiennych z przesunięciem , 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()

Brzmi to jak problem z asynchroniczną konwersją częstotliwości próbkowania. Aby przekonwertować jedną częstotliwość próbkowania na drugą, możemy obliczyć ciągłą reprezentację czasową sygnału, wykonując interpolację sinc, a następnie ponownie próbkować z naszą nową częstotliwością próbkowania. To, co robisz, nie różni się zbytnio. Musisz ponownie próbkować sygnał, aby mieć ustalone czasy próbkowania.

Sygnał ciągłego czasu można obliczyć przez zwoje każdej próbki funkcją sinc. Ponieważ funkcja sinc trwa w nieskończoność, używamy czegoś bardziej praktycznego, jak cynk okienny o praktycznej skończonej długości. Problem polega na tym, że ponieważ próbki poruszają się w czasie, może być konieczne użycie cynku z innym przesunięciem fazowym dla każdej próbki podczas ponownego próbkowania.

Ciągły sygnał czasu z próbkowanego sygnału:

$x(t) = \sum\limits_{n=-\infty}^\infty x[n]sinc({t-nT_s\over T_s})$

gdzie to czas próby. W twoim przypadku twój czas próbkowania nie jest jednak ustalony. Myślę więc, że musisz zastąpić go czasem próby w tej próbce.$T_s$

$x(t) = \sum\limits_{n=-\infty}^\infty x[n]sinc({t-nT_s[n]\over T_s[n]})$

Z tego można ponownie próbkować sygnał:

$y[n] = x(nT_{ns}$ )

gdzie jest pożądanym czasem próbkowania.$T_{ns}$

Zestawiając to w całość otrzymujesz:

$y[m] = \sum\limits_{n=-\infty}^\infty x[n]sinc({mT_{ns}-nT_s[n]\over T_s[n]})$

Ponieważ nie jest to przyczynowe ani podatne na działanie, funkcję sinc można zastąpić funkcją skończonego wsparcia i odpowiednio dostosować granice sumowania.

Niech jądro (t) będzie oknem okienkowym lub inną podobną funkcją o długości 2k, a następnie:

$y[m] = \sum\limits_{n=-k}^k x[n]kernel({mT_{ns}-nT_s[n]\over T_s[n]})$

Mam nadzieję, że to pomaga ... ale po drodze mogłem popełnić błąd i może to być trochę intensywne z matematyki. Aby uzyskać więcej informacji, zaleciłbym zbadanie konwersji częstotliwości próbkowania. Być może ktoś tutaj mógłby udzielić lepszego wyjaśnienia lub rozwiązania.

Myślę, że odpowiedź Jacoba jest bardzo wykonalna.

Prostszą metodą, która prawdopodobnie nie jest tak dobra pod względem wprowadzania zniekształceń, jest interpolacja wielomianowa. Użyłbym interpolacji liniowej (łatwej, niezbyt dobrej pod względem wydajności sygnału) lub splajnu sześciennego (wciąż niezbyt twardego, lepszej wydajności sygnału), aby wytwarzać próbki w dowolnym momencie z dowolnych próbek czasu.

(Miesiąc później) istnieją dwie główne opcje dla dowolnej metody interpolacji:
1) liczba punktów danych najbliższych brakującego punktu do użycia, 2 4 6 ... 2) klasa funkcji bazowych do użycia: liniowa, wielomianowa, sinus-cosinus (Fouriera), kawałek sześcienny (splajn B lub splajn interpolujący), sinc-like ... (Wybór 0 dotyczy tego, czy użyć metody i kodu innej osoby, czy też zrób to sam.)$Nnear$

Dopasowanie linii prostej do punktów jest łatwe: 2 punkty [-1, ], [1, ]: oszacowanie punkty ze średnią : średnia ogólna : patrz np. Przepisy numeryczne str. 781: dopasuj linię i oszacuj . Można dopasować kwadratykę, sześcienne, sinus-cosinus ... w ten sam sposób.$Nnear$
$y_{-1}$$y_1$
$\qquad$[ x i , y i ] x i = $y_0 \sim (y_{-1} + y_1) / 2$
$[x_i, y_i]$$x_i = 0$
y i [ x i , y i$\qquad y_0 \sim$$y_i$
$[x_i, y_i]$
$\qquad$y 0 ∼ $a + b x$$y_0 \sim a$
$\qquad$

Rozumiem, że macie jednolicie rozmieszczone dane z kilkoma brakującymi punktami, prawda?
Jak dobrze interpolacja liniowa działa w tym przypadku?
Cóż, spróbujmy cos z = 0,25: 1 0 -1 0 1 0 -1 0 ... 2 sąsiadów dowolnego punktu średnio do 0, strasznie. 4 sąsiadów: średnia z [1 0 (brak -1) 0 1] = 1/2, strasznie. (Wypróbuj filtr 4-sąsiedni [-1 3 3 -1] / 4.)$2 \pi f t$$f$

Interpolacja liniowa z 4, 6 lub 8 sąsiadami może działać wystarczająco dobrze dla twoich danych.
Proponuję zacząć od metody, którą dokładnie rozumiesz, zanim zanurzysz się w splajny, podobne do cynku ... chociaż one też mogą być zabawne.

---

Inną, zupełnie inną metodą jest odwrotne ważenie odległości . Jest łatwy do wdrożenia (patrz: idw-interpolation-with-python na SO), działa w 2d 3d i nowszych , ale jest trudny do analizy teoretycznej.

(Oczywiście, NO metoda pojedynczy interpolacja może ewentualnie dopasować zillions kombinacji
[sygnału, hałas błędu metryki, funkcja TEST], które występują w rzeczywistości.
Są to kolejne sposoby na świecie, z większą ilością pokręteł, niż funkcji testowych.
Niemniej galerię metod i funkcji testowych może się przydać).

Jeśli pracujesz z matlabem, możesz to zrobić, pracując z timeseries.

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

Zanim zaczniesz jakieś egzotyczne przetwarzanie, możesz wypróbować coś tak prostego (pseudo kod - bez interpolacji, ale można to dodać)

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]
   }
}

Odpowiedź IMHO Datageist jest poprawna, odpowiedź Jacoba nie jest. Łatwym sposobem na sprawdzenie tego jest to, że sugerowany algorytm bazy danych gwarantuje interpolację oryginalnych próbek (przy założeniu nieskończonej precyzji liczbowej), podczas gdy odpowiedź Jacoba nie.

• W przypadku równomiernego próbkowania zestaw funkcji sinc jest ortogonalny: jeśli każda przesunięta funkcja sinc jest dyskretyzowana na próbkach wejściowych, tworzą one nieskończoną macierz tożsamości. Jest tak, ponieważ sin (n pi) / (n pi) wynosi zero dla wszystkich n, z wyjątkiem n = 0.
• Jednak tej właściwości nie można po prostu ekstrapolować na przypadek niejednorodny: podobny zestaw funkcji sinc, dyskretny w stosunku do próbek wejściowych, da macierz niebanalną. Stąd wkład z otaczających próbek nie będzie wynosił zero, a rekonstrukcja nie będzie już interpolować przez próbki wejściowe.