#Lectures #DSP

Преобразование Гильберта — конспект лекции

16 Apr

by Vladimir Leonidov 3

Аналитический сигнал

В жизни мы обычно имеем дело с действительными сигналами (ток, напряжение, мощность и т.п.), однако, иногда удобно представлять их в виде набора комплексных чисел. Это позволяет точно оценить амплитуду сигнала, его фазу, знак частоты (положительная или отрицательная). Также, для обработки комплексных чисел имеется мощный математический аппарат. Комплексные сигналы широко применяются в цифровой обработке сигналов: получение огибающей, определение мнгновенной частоты, квадратурная модуляция, системы связи, обработка радиолокационных сигналов и другое.

Давайте сразу перейдём к примеру. Создадим сигнал x[n], который представляет из себя косинусоиду частотой 1 кГц, затем создадим второй сигнал, сдвинутый относительно x[n] на $90^\circ$ (или синусоиду) и назовём его y[n]:

clear, clc, close all

Fs = 40000; % частота дискретизации
ts = 0 : 1/Fs : 0.005-1/Fs; % временные отсчёты

x = cos(2*pi*1000*ts);
y = sin(2*pi*1000*ts);

subplot(2,1,1)
plot(x), grid on, title('x[n]')
xlabel('Время')
ylabel('Амплитуда')

subplot(2,1,2)
plot(x), grid on, hold on, title('x[n] и y[n]')
plot(y), grid on
xlabel('Время')
ylabel('Амплитуда')

Полученные сигналы показаны ниже:

Исходный сигнал и сигнал, сдвинутый на 90°

Теперь построим график в трёхмерной системе координат, по оси x будет исходный сигнал x[n] (действительная часть), по оси y сигнал y[n] (мнимая часть), по оси z — время. Также на отдельных графиках покажем проекции полученного сигнала на каждую из плоскостей трёхмерной системы координат:

figure
subplot(2,2,1)
plot3(x,y,ts), grid on, title('Аналитический сигнал')
xlabel('Действительная часть')
ylabel('Мнимая часть')
zlabel('Время')

subplot(2,2,2)
plot(x,y), grid on, title('Проекция x-y')
xlabel('Действительная часть')
ylabel('Мнимая часть')

subplot(2,2,3)
plot(x,ts), grid on, title('Проекция x-z')
xlabel('Действительная часть')
ylabel('Время')

subplot(2,2,4)
plot(y,ts), grid on, title('Проекция y-z')
xlabel('Мнимая часть')
ylabel('Время')

Запускаем скрипт и анализируем полученные результаты:

Рассмотрим первый график: мы видим спираль, вращающуюся вокруг начала координат в плоскости x-y и поднимающуюся “вверх” вдоль оси z. Спираль вращается против часовой стрелки, это значит, что частота сигнала положительная. Правее показана окружность — это проекция данной спирали на плоскость x-y. Радиус данной окружности не изменяется и равен единице, а значит и амплитуда нашего сигнала на протяжении всего времени также не изменяется и равна единице. Ниже показаны проекции спирали на плоскости y-z и x-z: это и есть наши исходные косинусоидальный и синусоидальный сигналы.

Таким образом, любой дискретный сигнал можно представить в виде набора комплексных чисел:

(1) $\begin{equation*} S[n] = I[n]+jQ[n] \end{equation*}$

где:

$I[n]$ — In-Phase — синфазная составляющая (действительная часть сигнала)
$Q[n]$ — Quadrature — квадратурная составляющая (мнимая часть сигнала)

Мы уже работали с комплексными числами, когда анализировали результаты ДПФ. Например, для расчёта амплитуды сигнала необходимо найти модуль комплексного числа $S[n]$ :

(2) $\begin{equation*} |S[n]| = \sqrt{I^2[n]+Q^2[n]} \end{equation*}$

А для расчёта фазы сигнала необходимо рассчитать арктангенс:

(3) $\begin{equation*} \phi[n] = arctg{\frac{Q[n]}{I[n]}} \end{equation*}$

Теперь предположим что,

(4) $\begin{equation*} S[n] = A \cdot cos(2\pi f_0 t + \phi) \end{equation*}$

где $f_0$ — несущая частота сигнала, $\phi$ — его фаза. Вспомним курс тригонометрии:

(5) $\begin{equation*} \cos(\alpha+\beta) = \cos(\alpha)\cos(\beta) - \sin(\alpha)\sin(\beta), \end{equation*}$

И преобразуем выражение (4):

(6) $\begin{equation*} A \cdot cos(2\pi f_0 t + \phi) = A \cdot \cos(2\pi f_0 t) \cdot \cos(\phi) - A \cdot \sin(2\pi f_0 t) \cdot \sin(\phi) \end{equation*}$

Т.к. для комплексных чисел справедливо:

(7) $\begin{equation*} I = A \cdot \cos(\phi) \end{equation*}$

(8) $\begin{equation*} Q = A \cdot \sin(\phi) \end{equation*}$

Делаем замену в выражении (6):

(9) $\begin{equation*} A \cdot cos(2\pi f_0 t + \phi) = I \cdot cos(2\pi f_0 t) - Q \cdot sin(2\pi f_0 t) \end{equation*}$

И получаем ещё одну форму записи аналитического сигнала:

(10) $\begin{equation*} S[n] = I[n] \cdot cos(2\pi f_0 t) - Q[n] \cdot sin(2\pi f_0 t) \end{equation*}$

Что такое аналитический сигнал — немного разобрались, теперь нужно понять, как его получить из привычного нам действительного сигнала. Эту задачу помогает решить дискретное преобразование Гильберта.

Дискретное преобразование Гильберта

Дискретное преобразование Гильберта представляет собой процедуру, которая на основе исходного сигнала $x[n]$ формирует сигнал $x_h[n]$ , сдвинутый по фазе относительно $x[n]$ на $90^\circ$ . Иными словами: данное преобразование представляет собой идеальный фазовращатель на $90^\circ$ , который называется фильтром Гильберта. Двигать фазу будем с помощью КИХ-фильтра. На предыдущих лекциях мы рассмотрели операцию свёртки, в случае с преобразованием Гильберта её можно записать следующим образом:

(11) $\begin{equation*} x_h[n] = x[n]*h[n] \end{equation*}$

где $h[n]$ — импульсная характеристика фильтра Гильберта, а “ $*$ ” обозначает операцию свёртки. Тогда, по теореме о свёртке, ДПФ выходного сигнала преобразования Гильберта есть произведение ДПФ исходного сигнала и ДПФ фильтра Гильберта:

(12) $\begin{equation*} X_h[\omega] = X[\omega] \cdot H[\omega] \end{equation*}$

Исходя из задачи сдвинуть фазу сигнала на $90^\circ$ , можно сформулировать следующие требования к фильтру $H[\omega]$ :

Все компоненты $X[\omega]$ с положительными частотами должны быть сдвинуты на $-90^\circ$
Все компоненты $X[\omega]$ с отрицательными частотами должны быть сдвинуты на $+90^\circ$

Или, если записать это в виде математического выражения:

(13) $\begin{equation*} H[\omega]=\begin{cases} j, & \text{при $\omega<0$}.\\ 0, & \text{при $\omega=0$}.\\ -j, & \text{при $\omega>0$}. \end{cases} \end{equation*}$

Графическое представление АЧХ и ФЧХ $H[\omega]$ из выражения (13):

Следует обратить внимание, что т.к. $H[0]=0$ , преобразование Гильберта убирает из сигнала постоянную составляющую.

Возьмём аналитический сигнал:

(14) $\begin{equation*} s[n] = x[n] + j x_h[n] \end{equation*}$

Спектр данного сигнала с учётом (12) будет иметь вид:

(15) $\begin{equation*} S[\omega] = X[\omega] + j X_h[\omega] = X[\omega] + j X[\omega] \cdot H[\omega] \end{equation*}$

Подставим в это выражение $H(\omega)$ из выражения (13), в результате получим:

(16) $\begin{equation*} S[\omega]=\begin{cases} X[\omega]+j \cdot j \cdot X[\omega], & \text{при $\omega<0$}.\\ X[0], & \text{при $\omega=0$}.\\ X[\omega]-j \cdot j \cdot X[\omega], & \text{при $\omega>0$}. \end{cases} \end{equation*}$

Упростим выражение (16):

(17) $\begin{equation*} S[\omega]=\begin{cases} 0, & \text{при $\omega<0$}.\\ X[0], & \text{при $\omega=0$}.\\ 2 \cdot X[\omega], & \text{при $\omega>0$}. \end{cases} \end{equation*}$

Из (17) можно сделать вывод: спектр аналитического сигнала при отрицательных частотах равен нулю, а при положительных частотах он равен спектру исходного сигнала, умноженному на два. Т.к. постоянная составляющая нас не интересует, $S[0]$ вычислять не будем.

На основе вышесказанного разработаем алгоритм. Чтобы получить аналитический сигнал $s[n]$ , заданный в выражении (14), из действительного сигнала $x[n]$ , необходимо выполнить следующую последовательность действий:

Вычислить ДПФ от сигнала $x[n]$ , в результате чего будет получен сигнал $X[\omega]$ .
Определить индексы отрицательных $f_n$ и положительных $f_p$ частот $X[\omega]$ .
Умножить на ноль (или приравнять к нулю) элементы с отрицательными частотами: $S[f_n] = X[f_n] \cdot 0$ .
Умножить на 2 элементы с положительными частотами: $S[f_p] = X[f_p] \cdot 2$ .
Вычислить ОДПФ от полученного в п.3 и п.4 сигнала $S[\omega]$ , в результате будет получен искомый сигнал $s[n]$ .

Проверим на практике? Разработаем скрипт, реализующий данный алгоритм (привязка кода к пунктам указана в комментариях):

clear, clc, close all

Fs = 40000;
ts = 0 : 1/Fs : 0.005-1/Fs;
N = length(ts);

% исходный сигнал
x = cos(2*pi*1000*ts);

% вычисление ДПФ, п. 1
X = fft(x);

% поиск индексов положительных и отрицательных частот, п. 2:
fp = 2:floor(N/2)+mod(N,2); % первая половина частот ДПФ - положительные
fn = ceil(N/2)+1+~mod(N,2):N; % вторая половина частот ДПФ - отрицательные

% применение выражения (17), п. 3 и п. 4:
S(fn) = X(fn)*0;
S(fp) = X(fp)*2;

% вычисление ОДПФ, п. 5:
s = ifft(S);

% строим графики:
subplot(2,1,1)
plot(x), grid on, title('Исходный сигнал x[n]')
xlabel('Время')
ylabel('Амплитуда')

subplot(2,1,2)
plot(real(s)), grid on, hold on
plot(imag(s)), grid on, title('Аналитический сигнал s[n]')
xlabel('Время')
ylabel('Амплитуда')
legend({'Действительная часть';'Мнимая часть'})

В результате выполнения скрипта получим графики:

Вычисление аналитического сигнала с помощью ДПФ

Сверху — исходный сигнал $x[n]$ , снизу — действительная и мнимая часть полученного аналитического сигнала $s[n]$ . Можно сравнить графики, полученные в данном примере, с графиками, полученными в начале лекции. Результаты получились идентичные, из чего можно сделать вывод о корректности работы алгоритма.

Дискретный фильтр Гильберта

Для того, чтобы можно было использовать фильтр Гильберта как классический цифровой фильтр, нужно рассчитать его импульсную характеристику, а затем подставить её в выражение (11). Расчёт импульсной характеристики цифрового фильтра мы уже проводили на одной из предыдущих лекций: для этого нам нужно взять частотную характеристику фильтра $H[\omega]$ из выражения (13) и взять от неё ОДПФ.

Запишем выражение ОДПФ для непрерывного сигнала:

(18) $\begin{equation*} x(t) = \int\limits_{-\infty}^{+\infty} X(f) e^{j 2 \pi f t} df \end{equation*}$

Вместо частоты в Герцах будем использовать круговую частоту: $\omega = 2 \pi f$ , тогда $df = d\omega\frac{1}{2\pi}$ . Т.к. частотная характеристика дискретной системы периодична и повторяется с периодом частоты дискретизации $f_s$ (или $\omega_s$ ), пределы интегрирования установим от $-\omega_s/2$ до $+\omega_s/2$ :

(19) $\begin{equation*} h(t) = \frac{1}{2\pi}\int\limits_{-\omega_s/2}^{+\omega_s/2} H(\omega) e^{j \omega t} d\omega \end{equation*}$

Теперь разобьём интервал интегрирования на две области: $-\omega_s/2$ до $0$ и от $0$ до $+\omega_s/2$ :

(20) $\begin{equation*} \begin{aligned} h(t) = \frac{1}{2\pi}\int\limits_{-\omega_s/2}^{0} j e^{j \omega t} d\omega + \frac{1}{2\pi}\int\limits_{0}^{+\omega_s/2} -j e^{j \omega t} d\omega = \\ = \frac{j}{2\pi j t}\left[e^{j \omega t}|_{-\omega_s/2}^0 - e^{j \omega t}|_0^{+\omega_s/2}\right] = \\ = \frac{1}{2\pi t}\left[ e^{j0} - e^{-j \omega_s t/2} - e^{j \omega_s t/2} + e^{j0} \right] =\\ = \frac{1}{2\pi t}\left[ 2 - 2\cos(\omega_s t/2) \right] = \frac{1}{\pi t}\left[ 1 - \cos(\omega_s t/2) \right] \end{aligned} \end{equation*}$

Следует обратить внимание, что при $t=0$ получается неопределённость $0/0$ , решить её можно с помощью правила Бернулли—Лопиталя:

(21) $\begin{equation*} \lim_{t \to 0} \frac{\left[ \cos(\omega_s t/2) \right]'}{(\pi t)'} = \lim_{t \to 0} \frac{\omega_s t sin(\omega_s t/2)}{2\pi} = 0 \end{equation*}$

Откуда можно сделать вывод, что $h(0) = 0$ . Для дискретного сигнала результат выражения (20) будет выглядеть следующим образом:

(22) $\begin{equation*} h[n] = \frac{1}{\pi n t_s}\left[ 1 - \cos(\omega_s n t_s/2) \right] \text{ при } n \neq 0; h[n]=0 \text{ при } n = 0 \end{equation*}$

Т.к. $\omega_s = 2 \pi f_s$ , а $t_s = 1/f_s$ , получим:

(23) $\begin{equation*} h[n] = \frac{f_s}{\pi n}\left[ 1 - \cos(\pi n) \right] \text{ при } n \neq 0; h[n]=0 \text{ при } n = 0 \end{equation*}$

Построим график импульсной характеристики фильтра Гильберта $h[n]$ для 32 отсчётов и $f_s=1$ :

clear,clc,close all

N = 32;
n = -N/2:N/2-1;

h = 1./(pi*n).*(1-cos(pi*n));
h(N/2+1) = 0;

stem(n,h), grid on, title('Импульсная характеристика фильтра Гильберта')
xlabel('n'), ylabel('h[n]')
xlim([-N/2 N/2-1])

Получим результат:

Импульсная характеристика фильтра Гильберта из 32 отсчётов

Теперь проанализируем АЧХ и ФЧХ полученного фильтра:

figure;
% расчёт ДПФ и частотной оси
H = fft(h);
f = (-N/2:N/2-1)*1/N;
% так как мы сдвинули h по оси x из отрицательной области
% на половину его длины вправо,
% применяем теорему о сдвиге ДПФ:
H = H.*exp(-1i*2*pi*f*N/2);

subplot(2,1,1)
plot(f,abs(H),'o-'), grid on
title('АЧХ фильтра Гильберта')
xlabel('\omega/\pi')
ylabel('|H[\omega]|')

phases=[-90;-45;0;45;90];
subplot(2,1,2)
plot(f,angle(H)*180/pi,'o-'), grid on
title('ФЧХ фильтра Гильберта')
xlabel('\omega/\pi')
ylabel('\phi(\omega)')
set(gca, 'YTick', phases)

Получим результат:

АЧХ и ФЧХ спроектированного фильтра Гильберта из 32 отсчётов

Следует обратить внимание, что $H[0] = 0$ — это особенность КИХ-фильтра с антисимметричной импульсной характеристикой. В нашем случае количество отсчётов фильтра чётное, в случае нечётного количества отсчётов мы также будем наблюдать $H[\omega_s/2] = 0$ . В связи с этим, на низких частотах АЧХ фильтра заваливается, что будет плохо сказываться на преобразовании низкочастотных составляющих. Также, в области полосы пропускания мы наблюдаем пульсации АЧХ — это эффект Гиббса, который мы изучали на лекции по цифровым фильтрам. Как с этим бороться, мы уже знаем: увеличить порядок фильтра и использовать взвешивание окном. Что мы и сделаем: ниже показаны АЧХ и ФЧХ фильтра Гильберта, взвешенного окном Хэмминга для $N = 32$ , $N = 64$ и $N = 512$ отсчётов.

АЧХ и ФЧХ фильтра Гильберта, взвешенного окном Хэмминга для разного количества отсчётов

Из графиков видно, что увеличение количества отсчётов и взвешивание импульсной характеристики фильтра Гильберта окном приводит к улучшению АЧХ и ФЧХ фильтра. Теперь полученные коэффициенты можно смело подставлять в выражение (11) и вычислять преобразование Гильберта $x_h[n]$ через операцию свёртки.

Однако, есть ещё один момент: сигнал, прошедший через КИХ-фильтр, задерживается на $D = (K-1)/2$ отсчётов, где $K$ — порядок цифрового фильтра. Поэтому, чтобы скомпенсировать линейную фазовую задержку выходного сигнала $x_h[n]$ относительно входного $x[n]$ , входной сигнал необходимо задержать на $D$ отсчётов. Рассмотрим пример для фильтра Гильберта 128 порядка:

clear,clc,close all
Nh= 128;
n = -Nh/2:Nh/2-1;
h = 1./(pi*n).*(1-cos(pi*n));
h(Nh/2+1) = 0;
h = h.*hamming(Nh)';

% формируем сигнал с частотой 1кГц
fs = 40000;
ts = 0 : 1/fs : 0.01-1/fs;
N = length(ts);
x = cos(2*pi*1000*ts);

% операция свёртки
y = conv(x,h);
% расчёт необходимой задержки основного сигнала
D = round((Nh-1)/2);
% временной сдвиг основного сигнала
xd = [zeros(1, D), x];

% отображение результатов
subplot(2,1,1)
plot(x), grid on, title('Исходный сигнал')
xlabel('Время')
ylabel('Амплитуда')
xlim([0 500])

subplot(2,1,2)
plot(xd), grid on, hold on
plot(y), title('Применение разработанного фильтра Гильберта к сигналу 1кГц')
xlabel('Время')
ylabel('Амплитуда')
legend({'Действительная часть';'Мнимая часть'})
xlim([0 500])

Результат выполнения скрипта показан ниже:

Применение разработанного фильтра Гильберта к сигналу 1кГц и компенсация фазовой задержки

И графика видно, что сформированная мнимая часть сигнала (оранжевый график снизу) имеет задержку относительно начала координат на 63 отсчёта, на которую также была сдвинута действительная часть сигнала (синий график снизу). В результате получили два синусоидальных сигнала, сдвинутых друг относительно друга на $90^\circ$ , которые являются компонентами аналитического сигнала.

Получение огибающей сигнала

Что такое преобразование Гильберта, разобрались. Но у многих наверняка до сих пор есть вопросы: “Зачем оно нужно? Ну получили сигнал, сдвинутый на $90^\circ$ , и что?”. Рассмотрим один из примеров применения: получение огибающей модулированного сигнала.

Создадим несущий сигнал частотой 1 кГц, который будет модулирован частотой 50 Гц и построим его график:

clear, clc, close all

fs = 40000;
ts = 0 : 1/fs : 0.05-1/fs;
N = length(ts);

% несущая частота
fc = cos(2*pi*1000*ts);
% модулирующий сигнал
fm = sin(2*pi*50*ts);
% модулированный сигнал
x = fc.*fm;

plot(ts,x), grid on, title('Амплитудно-модулированный сигнал')
xlabel('Время'), ylabel('Амплитуда')

График полученного сигнала:

Несущая 1 кГц, модулированная сигналом 50 Гц

Теперь разработаем ФНЧ со следующими параметрами:

Fs = 40000 Гц
Fpass = 500 Гц
Fstop = 2000 Гц
Astop = 60 дБ

Затем возьмём модуль модулированного сигнала x и применим к нему вышеуказанный фильтр (функцию фильтра я здесь показывать не буду, она генерируется с помощью пакета filterDesigner):

%% Использование ФНЧ
lp_f = lp_filter;
y = filter(lp_f.Numerator,1,abs(x));

figure
plot(ts,x), grid on, hold on
plot(ts,y), ylim([-1 1])
title('Получение огибающей с помощью ФНЧ')
xlabel('Время'), ylabel('Амплитуда')
legend({'Несущая частота';'Огибающая'})

Получим график:

Из графика видно, что огибающая (оранжевый сигнал) напоминает реальную огибающую нашего сигнала только отдалённо. Если мы будем менять параметры фильтра, мы всё равно не сможем получить идеальный сигнал, повторяющий модуль синусоиды 50 Гц с единичной амплитудой. Попробуем применить полученные сегодня знания и получить огибающую сигнала с помощью преобразования Гильберта:

%% Преобразование Гильберта
z = hilbert(x); % стандартная функция SP Toolbox

figure
plot(ts,x), grid on, hold on
plot(ts,abs(z)), ylim([-1 1])
title('Получение огибающей с помощью преобразования Гильберта')
xlabel('Время'), ylabel('Амплитуда')
legend({'Несущая частота';'Огибающая'})

Результат показан ниже:

Получение огибающей с помощью преобразования Гильберта

Другое дело, правда?

Выводы

Итак, мы рассмотрели два способа расчёта аналитического сигнала:

Через прямое и обратное ДПФ.
Через операцию свёртки с использованием импульсной характеристики фильтра Гильберта.

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

Также на практике показали, как можно применить преобразование Гильберта для получения огибающей модулированного сигнала.

Скачать конспект в pdf: Hilbert Transform Lecture – V.V. Leonidov.pdf

Vladimir Leonidov

Visit Card