La forma correcta de realizar la reducción de resolución de Sinc (reducción de resolución de DFT) para señales discretas muestreadas uniformemente con un número finito de muestras

Jan 04 2021

Dada una señal $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ ¿Cuál sería la forma correcta de reducirlo en el dominio de la frecuencia (interpolación Sinc)?

Respuestas

3 Royi Jan 04 2021 at 19:58

Interpolación en frecuencia (dominio DFT)

La implementación es bien conocida. En MATLAB será algo como:

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

Por eso nos ocupamos de 2 casos aquí:

Upsample Agregamos
cero muestras a la parte central de la DFT para que coincida con el número de muestras de la salida ( numSamplesO).
Nos ocupamos del caso de que el número de muestras de entrada ( numSamples) sea par. En ese caso, dividimos la muestra de Nyquist ($ X \left[ N / 2 \right] $) en 2 donde $ N $ es el número de entrada de muestras.
Downsample Eliminamos
muestras de la parte central de la DFT para que coincidan con el número de muestras de la salida ( numSamplesO).
Nos ocupamos del caso de que el número de muestras de salida ( numSamplesO) sea par. En ese caso, dividimos la muestra de Nyquist ($ X \left[ M / 2 \right] $) en 2 donde $ M $ es el número de salida de muestras.

La pregunta es, ¿por qué lo hacemos de esta manera? ¿Por qué el factor de interpolación interpFactor? ¿Dónde el factor de división de$ 0.5 $¿viene de?
Para responder a eso, debemos recordar que la DFT es básicamente la Serie Discreta de Fourier (DFS).
Esto significa que la suposición más importante es que los datos son periódicos tanto en el dominio del tiempo como en el de la frecuencia.

Ahora, dado que el DFT es básicamente el DFS, la forma natural de interpolar una señal dentro de su período sería utilizando la Serie de Fourier.

Antes de entrar en detalles, definamos 2 conjuntos de números enteros que se utilizarán para definir los valores de los índices:

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

Esto significa, para una señal con anchos de banda máximos de $ \frac{1}{2 T} $ muestreado por el teorema de muestreo para $ t \in \left[ 0, N T \right) $ dónde $ T $ es el período de muestreo y $ P = N T $ es el período de función:

$$ \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 fórmula anterior funciona para el caso par $ N = 2 l, \; l \in \mathbb{N} $ y por el extraño caso $ N = 2 l + 1, \; l \in \mathbb{N} $. Lo anterior define la conexión entre los coeficientes DFT y los coeficientes de la serie de 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} $$

Pero tampoco hay nada que nos impida utilizar otros puntos de muestreo para cualquier conjunto $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ dónde $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Lo que da$ 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} } $ para $ t \in \left[ 0, N T \right) $. Esto funcionará para señales complejas y reales.
Para señales reales,$ x \left( t \right) \in \mathbb{R} $también podemos usar la forma coseno del 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} $$

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

Entonces, ahora debemos pensar en lo que vimos aquí y cómo se relaciona con el algoritmo anterior.
Primero, debemos prestar atención a que el truco principal aquí es que la forma nativa de la DFT debe ser cuando el índice va$ k \in \mathcal{K}_{DFT}^{N} $. Entonces es más fácil ver la conexión con los orígenes de la serie discreta de Fourier ( DFS ) del DFT .

Observación : En la práctica, la DFT se define (y calcula) con$ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.

Si elegimos el conjunto de la cuadrícula de tiempo uniforme de salida $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ estar en la forma $ {t}_{m} = m {T}_{s} $ donde la tasa de muestreo superior (nos ocuparemos de la reducción de muestreo más adelante) $ q = \frac{M}{N} \geq 1 $entonces está claro lo que se debe hacer mirando el IDFT para recuperar una cuadrícula:

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

Ahora tenemos que hacer que esto coincida con la fórmula de interpolación de arriba. Dado que es una transformación lineal multiplicándola por$ q $se encargará de la constante. También podemos notar que$ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ por lo tanto, estableciendo:

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

Desde el $ N $ periodicidad de la DFT podemos escribir la interpolación final para una cuadrícula uniforme de tiempo con factor de interpolación de $ 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}} $$

Dónde $ \hat{X} \left[ k \right] $ Se define como:

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

Que es exactamente lo que hicimos en el código de muestra anterior.

¿Qué pasa con la reducción de la resolución? Bueno, podemos usar la misma intuición en el dominio DFT que muestra el código. Esto se debe básicamente a que la interpolación que utiliza los coeficientes de la serie de Fourier no es más que una multiplicación en el dominio de la frecuencia por el Kernel de Dirichlet, que es el equivalente periódico del$ \operatorname{sinc} \left( \cdot \right) $función. Ésta es también la intuición del$ \frac{1}{2} $factor, ya que multiplicamos con una rectagle con valor 1 en el dominio de frecuencia que tiene discontinuidad de salto . De hecho, la Serie de Fourier converge al valor medio del salto en las discontinuidades. Desde que pasamos de$ 1 $ a $ 0 $, significa que el valor en el salto es $ 0.5 $.
Por lo tanto, el código de muestreo ascendente y descendente anterior solo aplica Dirichlet Kernel a los datos de acuerdo con la frecuencia de muestreo de la entrada, en el caso de muestra ascendente y la salida en el caso de muestra descendente.

Otro método para reducir la muestra sería aumentar la muestra a un factor entero del número de salida de muestras. Luego use diezmado (Tomar cada ... muestra) para obtener las muestras. El 2 coincidirá en el caso de que los datos no tengan energía en la frecuencia entre la frecuencia baja y la frecuencia de muestreo. Si es así, no coincidirán.