Документация Engee
Notebook

Теорема Котельникова

Теорема Котельникова (также известная как теорема отсчётов или теорема Найквиста–Шеннона) — это фундаментальное утверждение в теории сигналов, которое определяет условие точного восстановления непрерывного сигнала по его дискретным отсчётам.

Если непрерывный сигнал $s(t)$ имеет ограниченный спектр (т.е. содержит только частоты в диапазоне от $0$ до $F_{\text{max}}$​), то он может быть точно восстановлен по своим отсчётам, взятым с частотой:

$$F_s \geq 2 F_{\text{max}} $$

где:

  • $F_s$— частота дискретизации (частота отсчётов),

  • $F_{\text{max}}$​ — максимальная частота в спектре сигнала.

На что стоит обратить внимание:

  1. Частота Найквиста — это минимальная частота дискретизации, при которой возможно точное восстановление сигнала:

    $$F_{\text{Найквиста}} = 2 F_{\text{max}} $$

    Если $F_s > 2 F_{\text{max}}$​, возникает наложение спектров (алиасинг), и сигнал невозможно восстановить без искажений.

  2. Восстановление сигнала происходит с помощью ряда Котельникова (интерполяционной формулы):

    $$s(t) = \sum_{n=-\infty}^{\infty} s(nT) \cdot \text{sinc}\left(\frac{t - nT}{T}\right) $$

где:

$T = \frac{1}{F_s}$ - интервал между отсчётами,

$\text{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$

Проверка основного условия теоремы

Смоделируем "непрерывный" сигнал (на самом деле такой же дискретный сигнал, но с очень плотно расположенными отсчётами, и затем будем его "оцифровывать":

In [ ]: 
using DSP;
fs0 = 10000;                    # Завышенная частота дискртизации для отображения "аналогового" сигнала
dt0 = 1/fs0;                    # Шаг временной сетки
endtime = 0.04;                 # Конечный момент времени
t0 = 0:dt0:endtime-dt0;         # Вектор отсчётов времени (временная сетка)
Fmax = 100;                    # Основная частота синусоидального сигнала в Гц
analog = cos.(2*pi*Fmax*t0);   # Синусоидальный "непрерывный" сигнал для "оцифровки"

Дискретный сигнал, заведомо удовлетворяющий условию теоремы Котельникова:

In [ ]: 
fs1 = 10 * Fmax;               # Частота дискретизации сигнала, удовлетворяющая теореме Котельникова    
dt1 = 1/fs1;
t1 = 0:dt1:endtime-dt1;         
stem1 = cos.(2*pi*Fmax*t1);    # Дискретный сигнал (отсчёты непрерывного сигнала)
plot(t0, analog, linewidth=4, label="Непрерывный")
plot!(t1, stem1, line=:stem, marker=:circle, linewidth=2, label="ЧД = 1000 Гц")
Out[0]:

А теперь создадим сигнал, удовлетворяющий условию в предельном случае ($F_{\text{Найквиста}} = 2 F_{\text{max}}$):

In [ ]: 
fs2 = 2 * Fmax;
dt2 = 1/fs2;
t2 = 0:dt2:endtime-dt2;
stem2 = cos.(2*pi*Fmax*t2);
plot(t0, analog, linewidth=3, label="Непрерывный")
plot!(t2, stem2, line=:stem, marker=:circle, linewidth=3, label="ЧД = 200 Гц")
Out[0]:

Наконец, проверим, как выглядит "оцифрованный" сигнал, если условие теоремы Котельникова не выполняется:

In [ ]: 
fs3 = 1.5 * Fmax;
dt3 = 1/fs3;
t3 = 0:dt3:endtime-dt3;
stem3 = cos.(2*pi*Fmax*t3);
plot(t0, analog, linewidth=4, label="Непрерывный")
plot!(t3, stem3, line=:stem, marker=:circle, linewidth=3, label="ЧД = 150 Hz")
Out[0]:

Проверим возможность восстановления сигнала с fs = 200 Гц. Для этого нам придётся добавить нули между отсчётами дискретного сигнала, а потом пропустить через формирующий фильтр нижних частот:

In [ ]: 
upsampled_stem2_chunk = [1; zeros(49); -1; zeros(49)];
upsampled_stem2 = repeat(upsampled_stem2_chunk, 4);
plot(t0, analog, linewidth=4, label="Непрерывный")
plot!(t0, upsampled_stem2, line=:stem, marker=:circle, linewidth=3, label="ЧД = 200 Гц, с нулями")
Out[0]:

Смоделируем восстановление сигнала при помощи цифрового фильтра нижних частот (модель ЦАП):

In [ ]: 
myfilt = digitalfilter(Lowpass(2*110/fs0), Elliptic(4, 0.01, 60));
H, w = freqresp(myfilt);
freq_vec = fs0*w/(2*pi);
plot(freq_vec, pow2db.(abs.(H))*2, linewidth=3,
     xlabel = "Частота, Гц",
     ylabel = "Амплитуда, дБ",
     title = "АЧХ фильтра нижних частот")
Out[0]:

Отобразим сигнал на выходе ЦАП:

In [ ]: 
DAC = filt(myfilt, upsampled_stem2);
plot(t0, analog, linewidth=4, label="Непрерывный")
plot!(t0, DAC * 25, linewidth=4, label="Выход ЦАП")
Out[0]:

Мы убедились, что даже в предельном случае соблюдения условия теоремы, восстановление формы возможно.

Алиасинг при несоблюдении условия теоремы Котельникова

Визуализируем эффект "наложения" отсчётов сигнала с большей основной частотой на отсчёты сигнала с меньшей частотой. Дискретные отсчёты "принадлежат" обоим сигналам, но после дискретизации мы увидим лишь сигнал меньшей частоты:

In [ ]: 
tvec = 0:0.01:10;
x1 = cos.(pi/5*tvec);
x2 = cos.(11*pi/5*tvec)
plot(tvec,x2, linewidth=3, label = "Исходный сигнал")

tnew = 0:10;
y1 = cos.(11*pi/5*tnew)
plot!(tnew,y1, line=:stem, marker=:circle, linewidth=3, label = "Дискретизация")
plot!(tvec,x1, linewidth=4, label = "Наложение", title = "Наложение, сигнал с частотой в 11 раз меньше")
Out[0]:

Если алиасинг уже произошёл, бороться с ним какими-либо методами пост-обработки бесполезно. Необходимо не допускать его возникновения, правильно выбирая частоту дискретизации и используя антиалиасинговые фильтры.

Теперь давайте оценим эффект алиасинга на примере звукового сигнала. Добавим вспомогательную функцию для прослушивания аудио:

In [ ]: 
using WAV
using Base64

function audioplayer(s, fs);
    buf = IOBuffer();
    wavwrite(s, buf; Fs=fs);
    data = base64encode(unsafe_string(pointer(buf.data), buf.size));
    markup = """<audio controls="controls" {autoplay}>
                <source src="data:audio/wav;base64,$data" type="audio/wav" />
                Your browser does not support the audio element.
                </audio>"""
    display("text/html", markup);
end

Меняя положение слайдера, можно из "зоны" соблюдения условия теоремы зайти в зону её несоблюдения. В крайнем правом положении при ЧД = 3000 Гц сигнал с основной частотой 2250 Гц на слух неотличим от сигнала с частотой 750 Гц:

In [ ]: 
var_fs = 3000 # @param {type:"slider",min:3000,max:5000,step:500}

freq1 = 750;
freq2 = 1500;
freq3 = 2250;

timevec = 0:1/var_fs:1 - 1/var_fs;

sig1 = cos.(2*pi*freq1*timevec);
sig2 = cos.(2*pi*freq2*timevec);
sig3 = cos.(2*pi*freq3*timevec);

audioplayer(sig1,var_fs)
audioplayer(sig2,var_fs)
audioplayer(sig3,var_fs)

Заключение

Мы познакомились с основами теоремы Котельникова, которая широко применяется в цифровой обработке сигналов (аудио, видео, телекоммуникации) и при проектировании аналого-цифровых преобразователей (АЦП). Для предотвращения алиасинга перед дискретизацией применяют антиалиасинговые фильтры, подавляющие частоты выше $F_s / 2$.

Кстати, теорема была независимо открыта:

  • В 1927 году Владимиром Котельниковым (СССР).

  • В 1928 году Гарри Найквистом (США).

  • Позже обобщена Клодом Шенноном в 1949 году.

В западной литературе её часто называют теоремой Найквиста–Шеннона, но в русскоязычной традиции закрепилось имя Котельникова.