Il modo corretto per eseguire il downsampling sincrono (downsampling DFT) per segnali discreti campionati uniformemente con un numero finito di campioni

Jan 04 2021

Dato un segnale $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ quale sarebbe il modo corretto per sottocampionarlo nel dominio della frequenza (interpolazione Sinc)?

Risposte

3 Royi Jan 04 2021 at 19:58

Interpolazione in frequenza (dominio DFT)

L'implementazione è ben nota. In MATLAB sarà qualcosa del tipo:

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

Quindi ci occupiamo di 2 casi qui:

Upsample
Aggiungiamo zero campioni alla parte centrale del DFT per abbinare il numero di campioni dell'output ( numSamplesO).
Ci occupiamo del caso in cui il numero di campioni immessi ( numSamples) sia pari. In tal caso abbiamo diviso il campione di Nyquist ($ X \left[ N / 2 \right] $) in 2 dove $ N $ è il numero di campioni di input.
Downsample Rimuoviamo i
campioni della parte centrale del DFT per abbinare il numero di campioni dell'output ( numSamplesO).
Ci occupiamo del caso in cui il numero di output di samples ( numSamplesO) sia pari. In tal caso abbiamo diviso il campione di essere Nyquist ($ X \left[ M / 2 \right] $) in 2 dove $ M $ è il numero di campioni di output.

La domanda è: perché lo facciamo in questo modo? Perché il fattore di interpolazione interpFactor? Da dove viene il fattore di scissione di$ 0.5 $vieni da?
Per rispondere a ciò dobbiamo ricordare che il DFT è fondamentalmente il Discrete Fourier Series (DFS).
Ciò significa che l'ipotesi più importante è che i dati siano periodici sia nel dominio del tempo che in quello della frequenza.

Ora, poiché il DFT è fondamentalmente il DFS, il modo naturale per interpolare un segnale nel suo periodo sarebbe usare la serie di Fourier.

Prima di entrare nei dettagli definiamo 2 serie di numeri interi che verranno utilizzati per definire i valori degli indici:

$$ \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} $$

Ciò significa che per un segnale con larghezze di banda massime di $ \frac{1}{2 T} $ campionato dal teorema di campionamento per $ t \in \left[ 0, N T \right) $ dove $ T $ è il periodo di campionamento e $ P = N T $ è il periodo della funzione:

$$ \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} $$

La formula sopra funziona per il caso pari $ N = 2 l, \; l \in \mathbb{N} $ e per il caso dispari $ N = 2 l + 1, \; l \in \mathbb{N} $. Quanto sopra definisce la connessione tra i coefficienti DFT e i coefficienti della serie di Fourier :

$$ {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} $$

Ma nulla ci impedisce di utilizzare altri punti di campionamento per qualsiasi set $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ dove $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Che dà$ 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} } $ per $ t \in \left[ 0, N T \right) $. Questo funzionerà per segnali complessi e reali.
Per segnali reali,$ x \left( t \right) \in \mathbb{R} $possiamo anche usare la forma coseno della 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} $$

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

Quindi ora dobbiamo riflettere su ciò che abbiamo visto qui e su come si collega all'algoritmo sopra.
Per prima cosa dobbiamo prestare attenzione al fatto che il trucco principale qui è che la forma nativa del DFT dovrebbe essere quando l'indice va$ k \in \mathcal{K}_{DFT}^{N} $. Quindi è più facile vedere la connessione alle origini Discrete Fourier Series ( DFS ) del DFT .

Nota : in pratica, il DFT è definito (e calcolato) con$ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.

Se abbiamo scelto il set della griglia temporale uniforme di output $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ essere nella forma $ {t}_{m} = m {T}_{s} $ dove il tasso di sovracampionamento (ci occuperemo del downsampling in seguito) $ q = \frac{M}{N} \geq 1 $quindi è chiaro cosa è necessario fare guardando l' IDFT per recuperare una griglia:

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

Ora dobbiamo fare in modo che corrisponda alla formula di interpolazione dall'alto. Poiché è una trasformazione lineare moltiplicandola per$ q $si prenderà cura della costante. Possiamo anche notare quello$ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ quindi impostando:

$$ \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} $$

Dal $ N $ periodicità del DFT possiamo scrivere l'interpolazione finale per una griglia di tempo uniforme con fattore di interpolazione di $ 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}} $$

Dove $ \hat{X} \left[ k \right] $ è definito come:

$$ \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} $$

Che esattamente cosa abbiamo fatto nel codice di sovracampionamento sopra.

E il downsampling? Bene, possiamo usare la stessa intuizione nel dominio DFT come mostra il codice. Questo è fondamentalmente perché l'interpolazione utilizzando i coefficienti della serie di Fourier non è altro che una moltiplicazione nel dominio della frequenza per il Dirichlet Kernel che è l'equivalente periodico del$ \operatorname{sinc} \left( \cdot \right) $funzione. Questa è anche l'intuizione per il$ \frac{1}{2} $fattore, poiché moltiplichiamo per un rettangolo con valore 1 nel dominio della frequenza che ha discontinuità di salto . Infatti la serie di Fourier converge al valore medio del salto alle discontinuità. Da quando andiamo da$ 1 $ per $ 0 $, significa che il valore al salto è $ 0.5 $.
Quindi il codice downsmaplign e upsampling sopra applica solo il Dirichlet Kernel ai dati in base alla frequenza di campionamento dell'ingresso, nel caso dell'upsample e dell'output nel caso del downsample.

Un altro metodo per eseguire il downsampling potrebbe essere l'upsampling a un fattore intero del numero di campioni in uscita. Quindi utilizzare la decimazione (Prendi ogni ... campione) per ottenere i campioni. Il 2 corrisponderà nel caso in cui i dati non abbiano energia nella frequenza tra la velocità bassa e la frequenza campionata. Se lo fa, non corrisponderanno.