Właściwy sposób na zmniejszenie próbkowania Sinc (DFT Downsampling) dla jednolicie próbkowanych sygnałów dyskretnych o skończonej liczbie próbek

Jan 04 2021

Dano sygnał $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ jaki byłby prawidłowy sposób na zmniejszenie próbkowania w dziedzinie częstotliwości (interpolacja Sinc)?

Odpowiedzi

3 Royi Jan 04 2021 at 19:58

Interpolacja częstotliwości (domena DFT)

Wdrożenie jest dobrze znane. W MATLABIE będzie to coś takiego:

if(numSamplesO > numSamples)
    % Upsample
    halfNSamples = numSamples / 2;
    if(mod(numSamples, 2) ~= 0) % Odd number of samples
        vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); zeros(numSamplesO - numSamples, 1, 'like', vXDft); vXDft((ceil(halfNSamples) + 1):numSamples)];
    else % Even number of samples -> Special Case
        vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; zeros(numSamplesO - numSamples - 1, 1, 'like', vXDft); vXDft(halfNSamples + 1) / 2; vXDft((halfNSamples + 2):numSamples)];
    end
else
    % Downsample
    halfNSamples = numSamplesO / 2;
    if(mod(numSamples, 2) ~= 0) % Odd number of samples
        vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); vXDft((numSamples - floor(halfNSamples) + 1):numSamples)];
    else % Even number of samples -> Special Case
        vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; vXDft((numSamples - halfNSamples + 2):numSamples)];
    end
end

Zajmujemy się więc tutaj 2 przypadkami:

Upsample
Dodajemy zero próbek do środkowej części DFT, aby dopasować liczbę próbek wyjścia ( numSamplesO).
Dbamy o przypadek, gdy wejściowa liczba próbek ( numSamples) jest parzysta. W takim przypadku podzielimy próbkę Nyquist ($ X \left[ N / 2 \right] $) na 2, gdzie $ N $ jest wejściową liczbą próbek.
Downsample Usuwamy
próbki środkowej części DFT, aby dopasować liczbę próbek wyniku ( numSamplesO).
Dbamy o przypadek, gdy wyjściowa liczba próbek ( numSamplesO) jest parzysta. W takim przypadku podzielimy próbkę na Nyquist ($ X \left[ M / 2 \right] $) na 2, gdzie $ M $ jest wyjściową liczbą próbek.

Pytanie brzmi, dlaczego robimy to w ten sposób? Dlaczego współczynnik interpolacji interpFactor? Gdzie współczynnik podziału$ 0.5 $pochodzić z?
Aby odpowiedzieć na to pytanie, musimy pamiętać, że DFT to w zasadzie dyskretna seria Fouriera (DFS).
Oznacza to, że najważniejszym założeniem jest okresowość danych zarówno w dziedzinie czasu, jak i częstotliwości.

Teraz, ponieważ DFT jest zasadniczo DFS, naturalnym sposobem interpolacji sygnału w jego okresie byłoby użycie szeregu Fouriera.

Zanim przejdziemy do szczegółów, zdefiniujmy 2 zestawy liczb całkowitych, które posłużą do zdefiniowania wartości indeksów:

$$ \begin{aligned} \mathcal{K}_{DFS}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lceil \frac{N - 1}{2} \right\rceil \right\} \\ \mathcal{K}_{DFT}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lfloor \frac{N - 1}{2} \right\rfloor \right\} \\ \end{aligned} $$

Oznacza to, że dla sygnału o maksymalnej szerokości pasma wynoszącej $ \frac{1}{2 T} $ próbkowane przez twierdzenie o próbkowaniu dla $ t \in \left[ 0, N T \right) $ gdzie $ T $ jest okresem pobierania próbek i $ P = N T $ to okres funkcji:

$$ \begin{aligned} x \left( t \right) {\Big|}_{t = n T} & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{P} } && \text{By Fourier Series} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{By the period of the function / series} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k n}{N} } && \text{Setting $ t = n T $} \\ & = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k n}{N} } && \text{The DFT} \end{aligned} $$

Powyższy wzór działa dla przypadku parzystego $ N = 2 l, \; l \in \mathbb{N} $ i dla dziwnego przypadku $ N = 2 l + 1, \; l \in \mathbb{N} $. Powyższe definiuje związek między współczynnikami DFT a współczynnikami szeregu Fouriera :

$$ {c}_{k} = \begin{cases} \frac{ X \left[ k \right ] }{2 N} & \text{ if } k = \frac{N}{2} \\ \frac{ X \left[ k \right ] }{2 N} & \text{ if } k = -\frac{N}{2} \\ \frac{ X \left[ k \right ] }{N} & \text{ if } k \notin \left\{\frac{N}{2}, -\frac{N}{2} \right\} \end{cases}, \; k \in \mathcal{K}_{DFS}^{N} $$

Ale nic też nie powstrzymuje nas przed użyciem innych punktów próbkowania dla dowolnego zestawu $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ gdzie $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Co daje$ x \left( t \right) = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k t}{N T} } $ dla $ t \in \left[ 0, N T \right) $. To zadziała w przypadku złożonych i rzeczywistych sygnałów.
W przypadku rzeczywistych sygnałów$ x \left( t \right) \in \mathbb{R} $możemy również użyć formy kosinusowej DFT :

$$ \begin{aligned} x \left( t \right) & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{From the above} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} \left| {c}_{k} \right| \cos \left( 2 \pi \frac{k t}{N T} + \angle {c}_{k} \right) && \text{Fourier series in its Cosine form} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Fourier series in its Cosine form} \\ & = \sum_{k = 0}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} {\alpha}_{k} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Using the DFT conjugate symmetry of a real signal} \end{aligned} $$

Gdzie $ {\alpha}_{k} = \begin{cases} 1 & \text{ if } k \in \left\{ 0, \frac{N}{2} \right\} \\ 2 & \text{ else } \end{cases} $.

Więc teraz musimy przemyśleć to, co tu widzieliśmy i jak to się ma do powyższego algorytmu.
Najpierw musimy zwrócić uwagę, że główna sztuczka polega na tym, że natywna forma DFT powinna występować w momencie, gdy indeks przechodzi$ k \in \mathcal{K}_{DFT}^{N} $. Wtedy łatwiej jest dostrzec związek z początkami dyskretnych szeregów Fouriera ( DFS ) DFT .

Uwaga : W praktyce DFT jest definiowany (i obliczany) za pomocą$ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.

Jeśli wybraliśmy zbiór wyjściowej siatki czasu jednolitego $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ być w formie $ {t}_{m} = m {T}_{s} $ gdzie częstotliwość upsamplingu (zajmiemy się downsamplingiem później) $ q = \frac{M}{N} \geq 1 $wtedy jest jasne, co należy zrobić, patrząc na IDFT, aby odzyskać siatkę:

$$ x \left[ m \right] = \frac{1}{M} \sum_{k = 0}^{M - 1} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} = \frac{1}{M} \sum_{k = - \left\lceil \frac{M - 1}{2} \right\rceil}^{\left\lfloor \frac{M - 1}{2} \right\rfloor} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$

Teraz musimy dopasować to do wzoru interpolacji z góry. Ponieważ jest to transformacja liniowa pomnożona przez$ q $zadba o stałą. Możemy to również zauważyć$ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ stąd poprzez ustawienie:

$$ \tilde{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \mathcal{K}_{DFT}^{N} \setminus \left\{ k \mid k = \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \notin \mathcal{K}_{DFT}^{N} \end{cases} $$

Z $ N $ okresowość DFT możemy zapisać ostateczną interpolację dla jednolitej siatki czasu o współczynniku interpolacji $ q $:

$$ x \left[ m \right] = \frac{q}{M} \sum_{k = 0}^{M - 1} \hat{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$

Gdzie $ \hat{X} \left[ k \right] $ definiuje się jako:

$$ \hat{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \left\{ 0, 1, \ldots, N - 1 \right\} \setminus \left\{ \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \in \left\{ N, N + 1, \ldots, M - 1 \right\} \end{cases} $$

Co dokładnie zrobiliśmy w powyższym kodzie upsample .

A co z downsample? Cóż, możemy użyć tej samej intuicji w domenie DFT, jak pokazuje kod. Dzieje się tak po prostu dlatego, że interpolacja przy użyciu współczynników szeregu Fouriera to nic innego jak mnożenie w dziedzinie częstotliwości przez jądro Dirichleta, które jest okresowym odpowiednikiem$ \operatorname{sinc} \left( \cdot \right) $funkcjonować. Taka jest również intuicja dotycząca$ \frac{1}{2} $współczynnik, ponieważ mnożymy przez prostokąt z wartością 1 w dziedzinie częstotliwości, która ma nieciągłość skoku . Rzeczywiście, szereg Fouriera zbiega się do średniej wartości skoku przy przerwach. Ponieważ wychodzimy z$ 1 $ do $ 0 $oznacza to, że wartość przy skoku wynosi $ 0.5 $.
Tak więc powyższy kod downsmaplign i upsampling stosuje jądro Dirichleta do danych zgodnie z częstotliwością próbkowania wejścia, w przypadku upsamplingu i wyjścia w przypadku downsample.

Inną metodą zmniejszenia próbkowania byłoby upsamplowanie do współczynnika całkowitego wyjściowej liczby próbek. Następnie użyj decymacji (Take every ... sample), aby pobrać próbki. 2 będzie pasować do przypadku, gdy dane nie mają energii w częstotliwości między niską częstotliwością a częstotliwością próbkowaną. Jeśli tak, nie będą pasować.