La bonne façon de procéder à un sous-échantillonnage de Sinc (sous-échantillonnage DFT) pour des signaux discrets uniformément échantillonnés avec un nombre fini d'échantillons
Étant donné un signal $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ quelle serait la bonne façon de le sous-échantillonner dans le domaine fréquentiel (interpolation Sinc)?
Réponses
Interpolation en fréquence (domaine DFT)
La mise en œuvre est bien connue. Dans MATLAB, ce sera quelque chose comme:
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
On s'occupe donc de 2 cas ici:
- Suréchantillonnage
Nous ajoutons zéro échantillon à la partie centrale de la DFT pour correspondre au nombre d'échantillons de la sortie (numSamplesO).
Nous prenons en charge le cas où le nombre d'échantillons d'entrée (numSamples) est pair. Dans ce cas, nous divisons l'échantillon de Nyquist ($ X \left[ N / 2 \right] $) en 2 où $ N $ est le nombre d'échantillons d'entrée. - Sous-
échantillonnage Nous supprimons les échantillons de la partie centrale de la DFT pour correspondre au nombre d'échantillons de la sortie (numSamplesO).
Nous prenons soin du cas où le nombre d'échantillons en sortie (numSamplesO) est pair. Dans ce cas, nous divisons l'échantillon de Nyquist ($ X \left[ M / 2 \right] $) en 2 où $ M $ est le nombre d'échantillons de sortie.
La question est, pourquoi le faisons-nous de cette façon? Pourquoi le facteur d'interpolation interpFactor? D'où vient le facteur de fractionnement de$ 0.5 $viens de?
Pour répondre à cela, nous devons nous rappeler que la DFT est essentiellement la série discrète de Fourier (DFS).
Cela signifie que les hypothèses les plus importantes sont que les données sont périodiques à la fois dans le domaine temporel et fréquentiel.
Maintenant, puisque le DFT est fondamentalement le DFS, le moyen naturel d'interpoler un signal dans sa période serait d'utiliser la série de Fourier.
Avant d'entrer dans les détails définissons 2 ensembles de nombres entiers qui serviront à définir les valeurs des indices:
$$ \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} $$
Cela signifie, pour un signal avec des bandes passantes maximales de $ \frac{1}{2 T} $ échantillonné par le théorème d'échantillonnage pour $ t \in \left[ 0, N T \right) $ où $ T $ est la période d'échantillonnage et $ P = N T $ est la période de fonction:
$$ \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 formule ci-dessus fonctionne pour le cas pair $ N = 2 l, \; l \in \mathbb{N} $ et pour le cas étrange $ N = 2 l + 1, \; l \in \mathbb{N} $. Ce qui précède définit la connexion entre les coefficients DFT et les coefficients de la 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} $$
Mais rien ne nous empêche non plus d'utiliser d'autres points d'échantillonnage pour n'importe quel ensemble $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ où $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Qui donne$ 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} } $ pour $ t \in \left[ 0, N T \right) $. Cela fonctionnera pour les signaux complexes et réels.
Pour des signaux réels,$ x \left( t \right) \in \mathbb{R} $nous pouvons également utiliser la forme cosinus du 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} $$
Où $ {\alpha}_{k} = \begin{cases} 1 & \text{ if } k \in \left\{ 0, \frac{N}{2} \right\} \\ 2 & \text{ else } \end{cases} $.
Alors maintenant, nous devons réfléchir à ce que nous avons vu ici et comment cela se rapporte à l'algorithme ci-dessus.
Tout d'abord, nous devons faire attention au fait que l'astuce principale ici est que la forme native de la DFT doit être lorsque l'index va$ k \in \mathcal{K}_{DFT}^{N} $. Ensuite, il est plus facile de voir la connexion aux origines de la série de Fourier discret ( DFS ) de la DFT .
Remarque : En pratique, la DFT est définie (et calculée) avec$ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.
Si nous avons choisi l'ensemble de la grille de temps uniforme de sortie $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ être dans la forme $ {t}_{m} = m {T}_{s} $ où le taux de suréchantillonnage (nous nous occuperons du sous-échantillonnage plus tard) $ q = \frac{M}{N} \geq 1 $alors il est clair ce qu'il faut faire en regardant l' IDFT pour récupérer une grille:
$$ 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}} $$
Nous devons maintenant faire correspondre la formule d'interpolation ci-dessus. Puisqu'il s'agit d'une transformation linéaire en la multipliant par$ q $prendra soin de la constante. On peut également remarquer que$ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ donc en définissant:
$$ \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} $$
Du $ N $ périodicité de la DFT on peut écrire l'interpolation finale pour une grille de temps uniforme avec un facteur d'interpolation 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}} $$
Où $ \hat{X} \left[ k \right] $ est défini comme:
$$ \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} $$
C'est exactement ce que nous avons fait dans le code de suréchantillonnage ci-dessus.
Qu'en est-il du sous-échantillonnage? Eh bien, nous pouvons utiliser la même intuition dans le domaine DFT que le code montre. C'est essentiellement parce que l'interpolation utilisant les coefficients de la série de Fourier n'est rien d'autre qu'une multiplication dans le domaine fréquentiel par le noyau de Dirichlet qui est l'équivalent périodique du$ \operatorname{sinc} \left( \cdot \right) $fonction. C'est aussi l'intuition du$ \frac{1}{2} $facteur, comme nous multiplions par un rectagle de valeur 1 dans le domaine fréquentiel qui a une discontinuité de saut . En effet, la série de Fourier converge vers la valeur moyenne du saut aux discontinuités. Puisque nous allons de$ 1 $ à $ 0 $, cela signifie que la valeur au saut est $ 0.5 $.
Ainsi, le code de downsmaplign et de suréchantillonnage ci-dessus applique simplement le noyau de Dirichlet aux données en fonction de la fréquence d'échantillonnage de l'entrée, dans le cas du suréchantillonnage et de la sortie dans le cas du sous-échantillon.
Une autre méthode pour sous-échantillonner serait le suréchantillonnage à un facteur entier du nombre d'échantillons de sortie. Ensuite, utilisez la décimation (Take every ... sample) pour obtenir les échantillons. Le 2 correspondra au cas où les données n'ont pas d'énergie dans la fréquence entre le débit bas et le débit échantillonné. Si c'est le cas, ils ne correspondent pas.
J'ajouterai le code MATLAB ...
Remarque : Cette réponse couvre également le suréchantillonnage . Veuillez envisager d'ouvrir une autre question sur le suréchantillonnage ou d'élargir celle-ci.