Надлежащий способ выполнения понижающей дискретизации Sinc (DFT понижающей дискретизации) для однородно дискретизированных дискретных сигналов с конечным числом отсчетов
Учитывая сигнал $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ каков был бы правильный способ уменьшить разрешение в частотной области (интерполяция Sinc)?
Ответы
Интерполяция по частоте (область ДПФ)
Реализация хорошо известна. В MATLAB это будет примерно так:
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
Итак, здесь мы позаботимся о двух случаях:
- Upsample
Мы добавляем нулевые отсчеты в центральную часть ДПФ, чтобы соответствовать количеству отсчетов на выходе (numSamplesO).
Мы позаботимся о том, чтобы входное количество выборок (numSamples) было четным. В этом случае мы разбиваем выборку Найквиста ($ X \left[ N / 2 \right] $) в 2, где $ N $ - входное количество выборок. - Даунсэмпл
Мы удаляем отсчеты центральной части ДПФ, чтобы соответствовать количеству отсчетов на выходе (numSamplesO).
Мы позаботимся о том, чтобы выходное количество выборок (numSamplesO) было четным. В этом случае мы разбиваем образец Найквиста ($ X \left[ M / 2 \right] $) в 2, где $ M $ - количество отсчетов на выходе.
Вопрос в том, почему мы так поступаем? Почему коэффициент интерполяции interpFactor? Где коэффициент разделения$ 0.5 $родом из?
Чтобы ответить на этот вопрос, мы должны помнить, что ДПФ - это, по сути, дискретный ряд Фурье (ДПФ).
Это означает, что наиболее важным допущением является периодичность данных как во временной, так и в частотной области.
Теперь, поскольку ДПФ в основном ДФС естественным образом интерполировать сигнал в период его будет с помощью рядов Фурье.
Прежде чем вдаваться в подробности, давайте определим 2 набора целых чисел, которые будут использоваться для определения значений индексов:
$$ \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} $$
Это означает, что для сигнала с максимальной полосой пропускания $ \frac{1}{2 T} $ выборка по теореме выборки для $ t \in \left[ 0, N T \right) $ где $ T $ период выборки и $ P = N T $ период функции:
$$ \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 $ т = п т $} \\ & = \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} $$
Приведенная выше формула работает для четного случая $ N = 2 l, \; l \in \mathbb{N} $ и для нечетного случая $ N = 2 l + 1, \; l \in \mathbb{N} $. Вышеупомянутое определяет связь между коэффициентами ДПФ и коэффициентами ряда Фурье :
$$ {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} $$
Но также ничего не мешает нам использовать другие точки выборки для любого набора. $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ где $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Который дает$ 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} } $ для $ t \in \left[ 0, N T \right) $. Это будет работать для сложных и реальных сигналов.
Для реальных сигналов$ x \left( t \right) \in \mathbb{R} $мы также можем использовать косинусную форму ДПФ :
$$ \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} $$
куда $ {\alpha}_{k} = \begin{cases} 1 & \text{ if } k \in \left\{ 0, \frac{N}{2} \right\} \\ 2 & \text{ else } \end{cases} $.
Итак, теперь нам нужно подумать над тем, что мы здесь видели, и как это соотносится с приведенным выше алгоритмом.
Во-первых, нам нужно обратить внимание на то, что основная хитрость здесь в том, что собственная форма ДПФ должна быть, когда индекс$ k \in \mathcal{K}_{DFT}^{N} $. Тогда легче увидеть связь с истоками дискретных рядов Фурье ( DFS ) ДПФ .
Примечание . На практике ДПФ определяется (и вычисляется) с помощью$ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.
Если выбрать набор выходной единой временной сетки $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ быть в форме $ {t}_{m} = m {T}_{s} $ где частота повышающей дискретизации (о понижении дискретизации мы позаботимся позже) $ q = \frac{M}{N} \geq 1 $тогда становится ясно, что нужно сделать, взглянув на IDFT, чтобы восстановить сетку:
$$ 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}} $$
Теперь нам нужно привести это в соответствие с формулой интерполяции, приведенной выше. Поскольку это линейное преобразование, умножающее его на$ q $позаботится о постоянном. Мы также можем заметить, что$ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ следовательно, установив:
$$ \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} $$
От $ N $ периодичности ДПФ, мы можем записать окончательную интерполяцию для равномерной сетки времени с коэффициентом интерполяции $ 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}} $$
куда $ \hat{X} \left[ k \right] $ определяется как:
$$ \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} $$
Именно это мы и сделали в приведенном выше коде апсэмпла .
А как насчет субдискретизации? Что ж, мы можем использовать ту же интуицию в области DFT, как показывает код. Это в основном потому, что интерполяция с использованием коэффициентов ряда Фурье - это не что иное, как умножение в частотной области ядром Дирихле, которое является периодическим эквивалентом$ \operatorname{sinc} \left( \cdot \right) $функция. Это тоже интуиция для$ \frac{1}{2} $множитель, так как мы умножаем на прямоугольник со значением 1 в частотной области, которая имеет разрыв скачка . Действительно, ряд Фурье сходится к среднему значению скачка на разрывах. Поскольку мы идем от$ 1 $ к $ 0 $, это означает, что значение скачка равно $ 0.5 $.
Таким образом, приведенный выше код downsmaplign и upsampling просто применяет ядро Dirichlet к данным в соответствии с частотой дискретизации входных данных в случае повышающей дискретизации и выходных данных в случае понижающей дискретизации.
Другой метод субдискретизации - это повышающая дискретизация до целочисленного коэффициента выходного количества выборок. Затем используйте прореживание (взять каждый ... образец), чтобы получить образцы. 2 будет соответствовать случаю, когда данные не имеют энергии на частоте между низкой частотой и частотой дискретизации. Если это так, они не будут совпадать.
Я добавлю код MATLAB ...
Замечание : этот ответ также касается повышающей дискретизации . Пожалуйста, подумайте об открытии другого вопроса о повышении дискретизации или расширении этого.