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

27

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?

FigBug
źródło
To jest post dla programistów. Powiedziano mi, że jest to lepsze miejsce do zapytania. programmers.stackexchange.com/questions/193795/…
FigBug
Co decyduje o tym, kiedy czujnik zgłosi odczyt? Czy wysyła odczyt tylko w przypadku zmiany odczytu? Prostym podejściem byłoby wybranie „wirtualnego przedziału próbkowania” (T), który jest po prostu mniejszy niż najkrótszy czas między wygenerowanymi odczytami. Na wejściu algorytmu zapisz tylko ostatni raportowany odczyt (CurrentReading). Na wyjściu algorytmu zgłaszaj CurrentReading jako „nową próbkę” co T sekund, aby usługa filtrująca lub graficzna otrzymywała odczyty ze stałą szybkością (co T sekund). Nie mam pojęcia, czy jest to odpowiednie w twoim przypadku.
user2718
Próbuje próbkować co 5 ms lub 10 ms. Ale jest to zadanie o niskim priorytecie, więc może zostać pominięte lub opóźnione. Mam dokładność pomiaru do 1 ms. Przetwarzanie odbywa się na komputerze, a nie w czasie rzeczywistym, więc powolny algorytm jest odpowiedni, jeśli łatwiej go zaimplementować.
FigBug
1
Przyjrzałeś się rekonstrukcji Fouriera? Istnieje transformacja Fouriera oparta na nierównomiernie próbkowanych danych. Typowym rozwiązaniem jest przekształcenie obrazu czterokierunkowego z powrotem do równomiernie próbkowanej dziedziny czasu.
mbaitoff
3
Czy znasz jakieś cechy podstawowego sygnału, z którego próbujesz? Jeśli nieregularnie rozmieszczone dane są nadal z dość wysoką częstotliwością próbkowania w porównaniu do szerokości pasma mierzonego sygnału, wówczas coś prostego, jak interpolacja wielomianowa do równomiernie rozłożonej siatki czasowej, może działać dobrze.
Jason R

Odpowiedzi:

21

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)=k=1Ny(kT)sin(π(tkT)/T)π(tkT)/T=k=1Ny(kT)sinc(tkTT).

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)

[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)].

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 nnTmnn

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:

Sinc vs Cubic Spline Rekonstrukcja próbek Nonuniform

(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()
Datageist
źródło
Niezła metoda i kod. Ale dla z kilkoma (np. [0 1 - 3 4 5 - 7 8] T), które, jak myślę, jest pytaniem OP, czy nie sinc s w macierzy wszystkie 0? Pewnie, że istnieją sposoby, aby to naprawić, ale. tj=jT
denis
Jak rozumiem, pytanie OP dotyczy próbek, które zostaną upuszczone i / lub opóźnione. Ta metoda jest po prostu przesadzonym układem równań, więc upuszczone próbki pojawiają się tylko jako nieznane (nie jako punkty danych o wartości 0). A może nie o to pytasz?
Datageist
Co się stanie, jeśli wszystkie są liczbami całkowitymi (T = 1)? Powiedzmy, że mamy punkty danych [ ] dla , zbiór niezerowych liczb całkowitych, np. {-1 1} lub {-2 -1 1 2}; czy interpolowane , niezależnie od - czy coś przeoczyłem? tjj,yjjJy0=0yj
denis
Jeśli częstotliwości próbkowania są dokładnie identyczne (w / brakujące punkty), wówczas macierz interpolacji będzie rzadka (ponieważ każde wyjście zależy tylko od jednego wejścia). Zasadniczo średni odsetek próbek niejednorodnych musi być większy niż jednolity wskaźnik rekonstrukcji. Innymi słowy, trzeba będzie odbudowywać w niższym tempie, aby „uzupełnić luki” (T> 1 na przykład). Rozumiem jednak twój punkt widzenia.
Datageist
2
Takie odpowiedzi to czyste złoto.
Ahmed Fasih
6

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)=n=x[n]sinc(tnTsTs)

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.Ts

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

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

y[n]=x(nTns )

gdzie jest pożądanym czasem próbkowania.Tns

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

y[m]=n=x[n]sinc(mTnsnTs[n]Ts[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]=n=kkx[n]kernel(mTnsnTs[n]Ts[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.

Jakub
źródło
Użycie funkcji sinc do odtworzenia ciągłej wersji sygnału wymaga, aby próbki były równomiernie rozmieszczone, więc funkcja sinc będzie musiała dostosować się do arbitralnego odstępu między próbkami. Może być raczej trudny do wdrożenia.
user2718
tak, nie byłoby to zbyt skuteczne, aby zrobić dokładnie tak, jak pokazano tutaj. Wymagałoby to obliczenia nowych współczynników jądra dla każdego innego czasu próbkowania. Można jednak obliczyć zbiór kilku jąder, a czas skwantyzować do jednego z nich. Wystąpiłby błąd wydajności związany z błędem kwantyzacji.
Jacob
Można również wstępnie obliczyć pojedynczą tabelę wyszukiwania sinc i interpolować punkty w tej tabeli.
jms
5

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.

Jim Clay
źródło
1
Odpowiedź wydaje się o wiele łatwiejsza do zaimplementowania niż odpowiedź Jacoba, więc najpierw to zrobiłem. Wygląda na to, że działa, ale nie przeprowadziłem jeszcze pełnego zestawu testów.
FigBug
1
@FigBug - Jeśli masz czas, dodaj komentarz ze swoim ostatecznym rozwiązaniem.
user2718,
2

(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
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

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.)f2πftf


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ć).

denis
źródło
1

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
Rhei
źródło
0

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]
   }
}
użytkownik2718
źródło
0

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.
pvv
źródło