Необходимо вычислить огибающую дискретного сигнала. Для этого использую преобразование Гильберта следующим образом:
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). Слева: преобразование Гильберта (Зеленый). Справа: ДПФ от преобразования Гильберта.

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

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

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

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

а) Не понял заклинание
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 ) Глюки и выбросы по краям неизбежны, ибо такова суть вещей. Посмотрите на формулу оригинального преобразования Гильберта там интеграл от минус бесконечности до плюс бесконечности....

To DRUID3:

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


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


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

To alex_os:

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


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

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

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

Вот такой фильтр у Вас получается
Нажмите для просмотра прикрепленного файла

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

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

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

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

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


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

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

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

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

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



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

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

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

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

to: petrov


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

Можно и в аналоге на синус-косинус сигнал умножать и фильтровать можно и в цифре.


Сергиенко А.Б.
Цифровая обработка сигналов.

http://lord-n.narod.ru/walla.html#SergijenkoCIS

Всем спасибо за ответы, наконец-то более менее разобрался в вопросе

вывод - кругом одни болваны...


Вы ошибаетесь.


Таааааак... Так Вам на C филтр Гильберта или его расчет?


теперь я спрошу, что есть фильтр Гильберта для Вас?


Хм...И нафига там возведение в квадрат??? А как Вам такой способ - продецимировать до двойной максимальной частоты спектра огибающей и взять модуль .



Квадратурный смеситель и "зачем" оно вообще...

Круто... а йа - подготовленный радиолюбитель

А в чем она - разница между аналоговым и цифровым сигналом?


P.S.: квадратурный смеситель - это хорошо, но вот вопрос зачем нам туда-сюда переходить в область комплексных чисел? Мы же не OFDM принимаем.

Фильтр Гильберт дешевле получается, для него нужен один вещественный фильтр, для всяких квадратурных штуковин нужно два вещественных фильтра.