Cách thích hợp để thực hiện Sinc Downsampling (DFT Downsampling) cho các tín hiệu rời rạc được lấy mẫu thống nhất với số lượng mẫu hữu hạn
Đưa ra một tín hiệu $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ cách chính xác để giảm mẫu trong miền tần số (nội suy Sinc) là gì?
Trả lời
Nội suy tần số (miền DFT)
Việc thực hiện được nhiều người biết đến. Trong MATLAB, nó sẽ giống như sau:
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
Vì vậy, chúng tôi xử lý 2 trường hợp ở đây:
- Upsample
Chúng tôi thêm 0 mẫu vào phần trung tâm của DFT để khớp với số lượng mẫu của đầu ra (numSamplesO).
Chúng tôi đề phòng trường hợp số lượng mẫu đầu vào (numSamples) là số chẵn. Trong trường hợp đó, chúng tôi tách mẫu Nyquist ($ X \left[ N / 2 \right] $) thành 2 nơi $ N $ là số lượng mẫu đầu vào. - Downsample
Chúng tôi loại bỏ các mẫu của phần trung tâm của DFT để phù hợp với số lượng mẫu của đầu ra (numSamplesO).
Chúng tôi quan tâm đến trường hợp số lượng mẫu đầu ra (numSamplesO) là số chẵn. Trong trường hợp đó, chúng tôi tách mẫu thành mẫu Nyquist ($ X \left[ M / 2 \right] $) thành 2 nơi $ M $ là số lượng mẫu đầu ra.
Câu hỏi đặt ra là tại sao chúng ta lại làm theo cách này? Tại sao lại là hệ số nội suy interpFactor? Hệ số phân tách của$ 0.5 $đến từ?
Để trả lời rằng chúng ta cần nhớ DFT về cơ bản là Chuỗi Fourier rời rạc (DFS).
Điều này có nghĩa là các giả định quan trọng nhất là dữ liệu được định kỳ cả trong miền thời gian và tần số.
Bây giờ, vì DFT về cơ bản là DFS nên cách tự nhiên để nội suy một tín hiệu trong khoảng thời gian của nó sẽ là sử dụng Chuỗi Fourier.
Trước khi đi vào chi tiết, hãy xác định 2 bộ số nguyên sẽ được sử dụng để xác định giá trị của các chỉ số:
$$ \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} $$
Điều này có nghĩa là đối với một tín hiệu có băng thông tối đa là $ \frac{1}{2 T} $ lấy mẫu bằng Định lý Lấy mẫu cho $ t \in \left[ 0, N T \right) $ Ở đâu $ T $ là khoảng thời gian lấy mẫu và $ P = N T $ là chu kỳ hàm:
$$ \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} $$
Công thức trên phù hợp với trường hợp chẵn $ N = 2 l, \; l \in \mathbb{N} $ và đối với trường hợp kỳ quặc $ N = 2 l + 1, \; l \in \mathbb{N} $. Phần trên xác định mối liên hệ giữa hệ số DFT và hệ số chuỗi 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} $$
Nhưng cũng không có gì ngăn cản chúng tôi sử dụng các điểm lấy mẫu khác cho bất kỳ tập hợp nào $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ Ở đâu $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Cái nào cho$ 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} } $ cho $ t \in \left[ 0, N T \right) $. Điều này sẽ hoạt động đối với các tín hiệu thực và phức tạp.
Đối với các tín hiệu thực,$ x \left( t \right) \in \mathbb{R} $chúng ta cũng có thể sử dụng dạng Cosine của 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} $$
Ở đâu $ {\alpha}_{k} = \begin{cases} 1 & \text{ if } k \in \left\{ 0, \frac{N}{2} \right\} \\ 2 & \text{ else } \end{cases} $.
Vì vậy, bây giờ chúng ta cần suy nghĩ về những gì chúng ta đã thấy ở đây và nó liên quan như thế nào đến thuật toán ở trên.
Trước tiên, chúng ta cần chú ý rằng mẹo chính ở đây là dạng gốc của DFT phải là khi chỉ mục đi$ k \in \mathcal{K}_{DFT}^{N} $. Sau đó, dễ dàng thấy mối liên hệ với nguồn gốc Dòng Fourier rời rạc ( DFS ) của DFT .
Ghi chú : Trong thực tế, DFT được định nghĩa (Và được tính toán) với$ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.
Nếu chúng tôi chọn tập hợp của lưới thời gian thống nhất đầu ra $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ ở dạng $ {t}_{m} = m {T}_{s} $ trong đó tỷ lệ lấy mẫu tăng lên (Chúng tôi sẽ xem xét việc lấy mẫu xuống sau này) $ q = \frac{M}{N} \geq 1 $thì rõ ràng những gì cần phải làm bằng cách xem IDFT để khôi phục lưới:
$$ 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}} $$
Bây giờ chúng ta cần làm cho điều này khớp với công thức nội suy từ trên. Vì nó là một phép biến đổi tuyến tính nhân nó với$ q $sẽ chăm sóc hằng số. Chúng tôi cũng có thể nhận thấy rằng$ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ do đó bằng cách thiết lập:
$$ \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} $$
Từ $ N $ tính tuần hoàn của DFT, chúng ta có thể viết phép nội suy cuối cùng cho một lưới thời gian thống nhất với hệ số nội suy là $ 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}} $$
Ở đâu $ \hat{X} \left[ k \right] $ được định nghĩa là:
$$ \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} $$
Mà chính xác những gì chúng tôi đã làm trong upsample mã trên.
Còn về downample thì sao? Chà, chúng ta có thể sử dụng cùng một trực giác trong miền DFT như đoạn mã hiển thị. Điều này về cơ bản là do phép nội suy sử dụng hệ số Fourier Series không là gì khác ngoài phép nhân trong miền tần số với Hạt nhân Dirichlet, tương đương tuần hoàn của$ \operatorname{sinc} \left( \cdot \right) $chức năng. Đây cũng là trực giác cho$ \frac{1}{2} $nhân tố, khi chúng ta nhân với một trực thăng có giá trị 1 ở miền tần số có gián đoạn nhảy . Thật vậy, Fourier Series hội tụ với giá trị trung bình của bước nhảy tại điểm gián đoạn. Vì chúng tôi đi từ$ 1 $ đến $ 0 $, nó có nghĩa là giá trị tại bước nhảy là $ 0.5 $.
Vì vậy, mã downmaplign và upsampling ở trên chỉ áp dụng Dirichlet Kernel cho dữ liệu theo tần số lấy mẫu của đầu vào, trong trường hợp mẫu lên và đầu ra trong trường hợp mẫu xuống.
Một phương pháp khác để lấy mẫu xuống sẽ lấy mẫu lên một hệ số nguyên của số lượng mẫu đầu ra. Sau đó sử dụng phương pháp decimation (Lấy từng ... mẫu) để lấy các mẫu. Giá trị 2 sẽ phù hợp với trường hợp dữ liệu không có năng lượng trong tần số giữa tốc độ thấp và tốc độ lấy mẫu. Nếu có, chúng sẽ không khớp.
Tôi sẽ thêm Mã MATLAB ...
Nhận xét : Câu trả lời này cũng bao gồm Upsampling . Vui lòng xem xét mở một câu hỏi khác về Lấy mẫu cao hoặc mở rộng câu hỏi này.