Вычисление огибающей, Вопрос по преобразованию Гильберта Опции

[iggylike]

Необходимо вычислить огибающую дискретного сигнала. Для этого использую преобразование Гильберта следующим образом:
1). Вычисляю БПФ X(f) по сигналу x(k).
2). Считаю преобразование Гильберта:
x'(k) = (2 / F_d) * Im (sum_{i=1}^{N/2}[X(n) * exp^{2*pi*k*n * i/ N}] ), k = 0..N-1,
где F_d - частота дискретизации.
3). Нахожу аналитический сигнал z = x + ix'.
4). Модуль |z| и будет огибающей.

На практике получается, что в начальные и конечные моменты времени |z| имеет большие скачки, определенно не являющиеся огибающей. Подскажите, является ли такой эффект нормальным результатом или я все же напортачил где-то в вычислениях? И если это нормально, то каким способом лучше склеить огибающие подряд идущих сигналов? Дело в том, что для разных сигналов длина этих "выбросов" разная.

На рисунке слева сигналы во временной области, справа ДПФ. Сигнал x[i]= sin(2*pi*100*i) * exp(-100 *i). Частота дискретизации 1000, длина x - 128.
Построчно:
1). Слева: Входной x (Красный) и расчитанная огибающая |z|(Голубой) На конце значительное расхождение с сигналом. Справа: ДПФ входного.
2). Слева: Обратное преобразование Фурье от предыдущего ДПФ (синий)(посчитал для проверки) и преобразование Гильберта(Зеленый). Справа ДПФ от ОДПФ .
3). Слева: преобразование Гильберта (Зеленый). Справа: ДПФ от преобразования Гильберта.

Спасибо за ответ

Сообщение отредактировал iggylike - Oct 29 2008, 16:47

Эскизы прикрепленных изображений

 

[iggylike]

alex_os:
в последнем вашем ответе, который ушел в небытие, вы говорили про такой способ:
1). Посчитать хороший фильтр гильберта (например, Парксом-МакКлеланом).
2). Взять ДПФ от фильтра. Как я понял, что бы получить частотную характеристику.
3). Преремножить с S, вместо того, что бы обнулять S.
4). Дальше все так же (...как я делал).

Вобщем, не очень понял суть этого метода. Если посчитать Парксом-МакКлеланом фильтр, то разве на выходе этого фильтра мы не получим мнимый сигнал?

Сообщение отредактировал iggylike - Nov 3 2008, 13:02

[alex_os]

alex_os:
в последнем вашем ответе, который ушел в небытие, вы говорили про такой способ:
1). Посчитать хороший фильтр гильберта (например, Парксом-МакКлеланом).
2). Взять ДПФ от фильтра. Как я понял, что бы получить частотную характеристику.
3). Преремножить с S, вместо того, что бы обнулять S.
4). Дальше все так же (...как я делал).

Вобщем, не очень понял суть этого метода. Если посчитать Парксом-МакКлеланом фильтр, то разве на выходе этого фильтра мы не получим мнимый сигнал?

Конечно. Вы же хотели в частотной области считать? Еще, если в частотной области считать может быть меньше MIPS потребуется (особенно если фильтр длинный).

[shf_05]

в своей практике я сделал так- расчитал в Матлаб ф-р Гильберта 33-го пор-ка, округлил коэф-ты, проверил АЧХ "квантованного" ф-ра, загнал туда сигнал, задержал исходный для согласования задержек, получил уровень "зеркальной" ч-ты порядка -25дБ отн-но исходного с-ла. поскольку у ф-ра коэф-ты через 1 нули- считать в 2 р. меньше, не требуется особой мощи от проца
ИМХО если это Вас устроит- то самый простой способ

Сообщение отредактировал shf_05 - Nov 5 2008, 06:07

[iggylike]

На сколько я понимаю, фильтр Гильберта имеет асимметричную частотную характеристику, тоесть обрезает отрицательные частоты, но я ни как не пойму, как его строить в матлабе.

Код

h = remez(30,[0.1 0.9],[1 1],'Hilbert');
f = fft ([zeros(1, 30), h, zeros(1, 30)]);
plot (abs(f));

Характеристика вс равно симметрична. Может потому что надо использовать firpm? Но у меня старый матлаб, и таковой нет

Эскизы прикрепленных изображений

 

[alex_os]

На сколько я понимаю, фильтр Гильберта имеет асимметричную частотную характеристику, тоесть обрезает отрицательные частоты, но я ни как не пойму, как его строить в матлабе.

Код

h = remez(30,[0.1 0.9],[1 1],'Hilbert');
f = fft ([zeros(1, 30), h, zeros(1, 30)]);
plot (abs(f));

Характеристика вс равно симметрична. Может потому что надо использовать firpm? Но у меня старый матлаб, и таковой нет

Чтобы отфильтровать "отрицательные" частоты коэффициенты должны комплексные:

Код

h = remez(30,[0.1 0.9],[1 1],'Hilbert');
h = 1j*h;
h(16) = 1; % единственный вещественный коэффициент =1  
f = fft ([zeros(1, 30), h, zeros(1, 30)]);
plot (abs(f));