A maneira adequada de fazer redução da amostragem de Sinc (redução da amostragem DFT) para sinais discretos amostrados uniformemente com número finito de amostras

Jan 04 2021

Dado um sinal $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ qual seria a maneira correta de reduzir a resolução no domínio da frequência (interpolação Sinc)?

Respostas

3 Royi Jan 04 2021 at 19:58

Interpolação em frequência (domínio DFT)

A implementação é bem conhecida. No 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

Portanto, cuidamos de 2 casos aqui:

Upsample
Adicionamos zero amostras à parte central do DFT para corresponder ao número de amostras da saída ( numSamplesO).
Cuidamos caso o número de entrada de amostras ( numSamples) seja par. Nesse caso, dividimos a amostra Nyquist ($ X \left[ N / 2 \right] $) em 2 onde $ N $ é o número de entrada de amostras.
Downsample Removemos
amostras da parte central do DFT para corresponder ao número de amostras da saída ( numSamplesO).
Cuidamos caso o número de saída de samples ( numSamplesO) seja par. Nesse caso, dividimos a amostra de Nyquist ($ X \left[ M / 2 \right] $) em 2 onde $ M $ é o número de saída de amostras.

A questão é: por que fazemos assim? Por que o fator de interpolação interpFactor? Onde está o fator de divisão de$ 0.5 $vem de onde?
Para responder a isso, precisamos lembrar que o DFT é basicamente a Série Discreta de Fourier (DFS).
Isso significa que as suposições mais importantes são os dados sendo periódicos no domínio do tempo e da frequência.

Agora, como o DFT é basicamente o DFS, a maneira natural de interpolar um sinal dentro de seu período seria usando a série de Fourier.

Antes de entrar em detalhes, vamos definir 2 conjuntos de números inteiros que serão usados para definir os valores dos í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} $$

Isso significa que, para um sinal com larguras de banda máximas de $ \frac{1}{2 T} $ amostrado pelo Teorema de Amostragem para $ t \in \left[ 0, N T \right) $ Onde $ T $ é o período de amostragem e $ P = N T $ é o período de função:

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

A fórmula acima funciona para o caso par $ N = 2 l, \; l \in \mathbb{N} $ e para o caso estranho $ N = 2 l + 1, \; l \in \mathbb{N} $. O acima define a conexão entre os coeficientes DFT e os coeficientes da série 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} $$

Mas também não há nada que nos impeça de usar outros pontos de amostragem para qualquer conjunto $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ Onde $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Que 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} } $ para $ t \in \left[ 0, N T \right) $. Isso funcionará para sinais complexos e reais.
Para sinais reais,$ x \left( t \right) \in \mathbb{R} $também podemos usar a forma cosseno do 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} $$

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

Portanto, agora precisamos pensar sobre o que vimos aqui e como isso se relaciona com o algoritmo acima.
Primeiro, precisamos prestar atenção que o truque principal aqui é que a forma nativa do DFT deve ser quando o índice vai$ k \in \mathcal{K}_{DFT}^{N} $. Então é mais fácil ver a conexão com as origens da Série Discreta de Fourier ( DFS ) do DFT .

Observação : Na prática, o DFT é definido (e calculado) com$ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.

Se escolhermos o conjunto da grade de tempo uniforme de saída $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ estar na forma $ {t}_{m} = m {T}_{s} $ onde a taxa de aumento da resolução (cuidaremos da redução da resolução mais tarde) $ q = \frac{M}{N} \geq 1 $então fica claro o que precisa ser feito olhando para o IDFT para recuperar uma grade:

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

Agora precisamos fazer isso corresponder à fórmula de interpolação acima. Uma vez que é uma transformação linear multiplicando-o por$ q $vai cuidar da constante. Também podemos notar que$ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ portanto, definindo:

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

De $ N $ periodicidade da DFT, podemos escrever a interpolação final para uma grade uniforme de tempo com fator de interpolação 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}} $$

Onde $ \hat{X} \left[ k \right] $ é definido 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} $$

Exatamente o que fizemos no código upsample acima.

Que tal reduzir a resolução? Bem, podemos usar a mesma intuição no domínio DFT que o código mostra. Isso ocorre basicamente porque a interpolação usando os coeficientes da série de Fourier nada mais é do que multiplicação no domínio da frequência pelo kernel de Dirichlet, que é o equivalente periódico do$ \operatorname{sinc} \left( \cdot \right) $função. Esta também é a intuição para o$ \frac{1}{2} $fator, conforme multiplicamos por um retângulo com valor 1 no domínio da frequência que tem descontinuidade de salto . Na verdade, a Série de Fourier converge para o valor médio do salto nas descontinuidades. Desde que partimos de$ 1 $ para $ 0 $, isso significa que o valor no salto é $ 0.5 $.
Portanto, o código de downsmaplign e upsampling acima apenas aplica Dirichlet Kernel aos dados de acordo com a frequência de amostragem da entrada, no caso de upsample e a saída no caso de downsample.

Outro método para reduzir a resolução seria aumentar a resolução para um fator inteiro do número de saída de amostras. Em seguida, use a decimação (Take every ... sample) para obter as amostras. O 2 corresponderá para o caso de os dados não terem energia na frequência entre a taxa baixa e a taxa amostrada. Se isso acontecer, eles não corresponderão.