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

[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

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

 

[DRUID3]

А просто модуль отсчетов нельзя? Я так понял сигнал то real...

У Вас довольно низкое отношение несущей(частоты заполнения) и максимальной частоты спектра огибающей, потому то метод так и "сбоит"...

Уточните, у Вас научное или техническое задание? Т.е. строго найти огибающую или сделать приемлемый АМ детектор?

[alex_os]

Необходимо вычислить огибающую дискретного сигнала. Для этого использую преобразование Гильберта следующим образом...

а) Не понял заклинание
x'(k) = (2 / F_d) * Im (sum_{i=1}^{N/2}[X(n) * exp^{2*pi*k*n * i/ N}] ), k = 0..N-1, это matcad чтоли?
Но огибающую можно посчитать гораздо проще без ДПФ.
b ) Глюки и выбросы по краям неизбежны, ибо такова суть вещей. Посмотрите на формулу оригинального преобразования Гильберта там интеграл от минус бесконечности до плюс бесконечности....

[iggylike]

To DRUID3:

А просто модуль отсчетов нельзя? Я так понял сигнал то real...

Про модуль отсчетов ничего не знаю Фильтр что ли? Но переделывать, если честно не хотелось бы, разве что, если действительно на много быстрее. Да и FFT использую заточенное под процеccор. Описание метода нашел тут: http://www.mathworks.com/products/demos/sh...ml?product=DS#2. Там описан еще один метод: возведение в квадрат и пропускание через НЧ фильтр, но он дает худший результат.

Уточните, у Вас научное или техническое задание? Т.е. строго найти огибающую или сделать приемлемый АМ детектор?

Это ТЗ, причем, кроме названия Envelope Analysis, в нем, по сути, ничего нет. Я это проинтерпретировал именно как выделение огибающей. Сигнал real.

У Вас довольно низкое отношение несущей(частоты заполнения) и максимальной частоты спектра огибающей, потому то метод так и "сбоит"...

Повысил частоту заполняющего сигнала, и действительно, метод начал давать хорошие результаты. Но все же, интересно, как бороться с неровностями по краям? Если решения не найдется, то прийдется юзать другой метод

To alex_os:

Не понял заклинание

Заклинание это смесь латеха с с++, вырвалось само собой) Там мнимая часть от неполного обратного преобразования Фурье с 1 по N/2 элемент (если положить, что все действительное БПФ лежит в массиве с индексами 0..N), удвоенная и деленная на частоту дискретизации.

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

Я тоже так подумал, но проблема в том, что увеличение размера преобразования не особо то повышает качество работы алгоритма.

Сообщение отредактировал iggylike - Oct 30 2008, 08:57

[sergunas]

Но все же, интересно, как бороться с неровностями по краям? Если решения не найдется, то прийдется юзать другой метод

не буду утверждать на все 100% но наверное можно попробовать побороться.
Смотрите Help Matlab на fdesign.hilbert

[alex_os]

To alex_os:
Заклинание это смесь латеха с с++, вырвалось само собой) Там мнимая часть от неполного обратного преобразования Фурье с 1 по N/2 элемент (если положить, что все действительное БПФ лежит в массиве с индексами 0..N), удвоенная и деленная на частоту дискретизации.
Я тоже так подумал, но проблема в том, что увеличение размера преобразования не особо то повышает качество работы алгоритма.

Ага теперь понял "заклинание" , т.е. фактически вы делаете прямое ДПФ, обнуляете отрицательные частоты отсчеты (отсчеты N/2..N-1) и делаете обратное ДПФ. Так делать не очень хорошо, потому что если Вы построите AЧХ такого фильтра то увидете что имеются не слабые пульсации как в полосе пропускания (на положительных частотах ) так, так и на отрицательных частотах.

Код

s = zeros(1, 128);
s(64) = 1;
S = fft(s);
S( 65:128 ) = zeros(1, 64);
s_gilb = ifft(S);
% дополняем сигнал нулями и делаем FFT большего размера чтобы посмотреть
% частотную хар-ку фильтра.
S_gilb = fft( [zeros(1,128), s_gilb, zeros(1,128)] );
plot(abs(S_gilb))

Вот такой фильтр у Вас получается

[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));

[DRUID3]

Блин, люди, у вас задача сделать тот или иной АМ детектор. Зачем эти преобразователи Гильберта?

Неужели если есть FFT заточенное под процессор(собственно а в чем проблема его "заточить"?) то нужно "тулить" БПФ даже в те места где оно совсем не нужно? Быстрое преобразование Фурье, быстро делает лишь преобразование Фурье

[iggylike]

Кругом, в том числе на www.dsprelated.com, народ говорит о преобразовании Гильберта, как о самом подходящем и универсальном способе выделения амплитудной огибающей. Но теперь есть одно но.

Блин, люди, у вас задача сделать тот или иной АМ детектор. Зачем эти преобразователи Гильберта?

Неужели если есть FFT заточенное под процессор(собственно а в чем проблема его "заточить"?) то нужно "тулить" БПФ даже в те места где оно совсем не нужно? Быстрое преобразование Фурье, быстро делает лишь преобразование Фурье

Я уже сошелся на мысли, что лучше использовать фильтр Гильберта, но есть одна проблема - реализация на си. В природе нашлась только одна библиотека, которая может синтезировать фильтр методом Паркса-МакЛиллана (автор Jake Janovetz), но именно фильтр Гильберта она строит некорректно. Есть еще какая то его модификазия в проекте Octave, но я так и не смог ее портировать в Visual Studio. Времени разбираться и писать самому, к сожалению, нет. Может подскажете, есть ли какие простые способы синтезировать фильтр с асимметричной ЧХ? Остается два варианта: использовать преобразование Гильберта "в лоб" с припиской, что частота заполнения должна быть где то в 50 раз больше самого сигнала, либо использовать метод возведения в квадрат и последующей низкочастотной фильтрации .

З.ы. При следующем вызове функции Джейка Яновеца:

Код

numBands = 3;
numTaps = 128;

desired[0] = 0;
desired[1] = 1;
desired[2] = 0;

weights[0] = 1;
weights[1] = 1;
weights[2] = 1;

bands[0] = 0;
bands[1] = 0.02;
bands[2] = 0.05;
bands[3] = 0.45;
bands[4] = 0.48;
bands[5] = 0.5;
  
remez(h, numTaps, numBands, bands, desired, weights, HILBERT);

получается ЧХ как на нижнем рисунке.

Сообщение отредактировал iggylike - Nov 6 2008, 09:45

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

 

[petrov]

Непонятно в чём проблема то взять и в матлабе в fdatool фильтр посчитать. Можно без гильбертов всяких квадратурным смесителем сигнал на нулевую частоту перенести и ФНЧ профильтровать, получим комплексную огибающую, хоть амплитуду считай хоть фазу.

[iggylike]

to: petrov

Можно без гильбертов всяких квадратурным смесителем сигнал на нулевую частоту перенести и ФНЧ профильтровать, получим комплексную огибающую, хоть амплитуду считай хоть фазу.

что за квадратурный смеситель? ткните носом! по образованию я математик, многих терминов не знаю. везде пишут про квадратурную и синфазную составляющую, а что оно такое нигде найти не могу. Это разве не к аналоговому сигналу относится?

Сообщение отредактировал iggylike - Nov 6 2008, 11:16