Помогите определить реакцию фильтра...

[Altemirx]

Всем привет.

Уважаемые, помогите, пожалуйста, определить реакцию фильтра с передаточной функцией вида H(jw) (ПФ сложная и может быть впоследствии изменена, пример: http://i.exponenta.ru/exponenta/2006/8/24/1707.gif ) на единичное воздействие Ф(t-1) (ступенчатая функция Хэвисайда, сдвинутая по оси времени на 1). Желательно получить аналитическое выражение F(t) этой реакции. Если нет, то хотя бы численное, так, чтобы можно было построить график этой реакции (тем же MathCad-ом, например). Если невозможно будет получить это решение, то необходимо найти хотя бы значение амплитуды на выходе фильтра в установившемся режиме (считаем, что фильтр устойчивый).

Совет по моделированию в Simulink-е не принимается, поскольку уже сделал это.
Пробовал произвести обратное преобразование Лапласа от отношения H(iw)/(exp(-s)/s), но в символьном виде даже Maple сие не может решить до конца Конкретно - нет решений для функций вида 1/(exp(-s)*s^n), где n - натуральные числа

[G^a]

Здравствуйте!

Я вижу примерно следующее быстрое численное (оценочное) решение Вашей задачи:

1) Переводим частотную характеристику $H(j\omega)$ в импульсную переходную характеристику $h(\tau)$.

2) Проводим свёртку между функцией Хевисайда и импульсной переходной характеристикой $h(\tau)$ - получаем решение.

3) Обрабатываем полученную реакцию - считаем амплитуду на выходе.

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

Конечно - это не будет точным решением, а всего лишь оценкой, но может Вам этого хватит?

И ещё два вопроса, если у Вас известна частотная характеристика, то почему Вы не проверите устойчивость фильтра явно, и зачем Вам Лаплас?

[Altemirx]

Цитата:

1) Переводим частотную характеристику в импульсную переходную характеристику .

К сожалению, я недавно этим занимаюсь, потому не нашёл в явной форме формулы для этого перехода. Ранее пользовался для простых передаточных функций: L-1(H(s)/(exp(-s)/s)), получал во временной области реакцию на единичный импульс (ответ на вопрос про Лапласа). Не могли бы Вы привести формулу для указанного Вами перехода.

Цитата:

2) Проводим свёртку между функцией Хевисайда и импульсной переходной характеристикой - получаем решение.

В продолжение предыдущего вопроса - если можно, пример формулы. Какой вид при этом будет иметь решение?

Цитата:

3) Обрабатываем полученную реакцию - считаем амплитуду на выходе.

Что понимается под обработкой?

Цитата:

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

Вы имеете ввиду - расположение его нулей и полюсов? Мне этого недостаточно. Необходимо иметь численное представление реакции фильтра. Далее - не только на функцию Хэвисайда.

Извиняюсь за наивность, но хочется подробного описания решения задачи.

[Altemirx]

Попробовал решить поставленную задачу операторным методом (учитывая, что Ф(t) <-> 1/s), но что-то не то получается (см. рис: http://i.exponenta.ru/exponenta/2006/8/24/f_1709.gif ), т.к. результат MC отличается от моделирования в Simulink. Help!

[G^a]

Добрый день!

По пунктам.

1) Частотная характеристика переводится в импульсную - обратным преобразованием Фурье. Обратное преобразование Фурье производите по формуле:
$h_{k}=2\Delta f\sum\Limits^{N-1}_{i=0}H_{i}e^{j2\pi\Delta f k\Delta t i}$

Выбирая отсчёты по частоте - помните про теорему Котельникова.

2) Свёртку входного сигнала с импульсной производите по формуле:

$y_{k}=\sum\Limits^{N-1}_{i=0}h_{i}x_{k-i}$

Как видно из формулы, решение получаете в виде временного ряда.

3) Под обработкой понимается нахождение всевозможных параметров: среднеквадратичного, эффективного, амплитудного значения вых. сигнала.

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

[Altemirx]

Спасибо за ответы, но несколько непонятно, что есть \Delta f, \Delta t и как получены значения H_{i}?
Как я понимаю: H_{i} - отсчёты передаточной функции (но она же комплексная!), взятые через интервал \Delta f (в Гц), а \Delta t - временнЫе интервалы для проекции на временнУю плоскость. Так? Разъясните пожалуйста. И ещё, как выбирается \Delta t?

[G^a]

Продолжаем!


$H_{i}$ - это отсчёты комплексной частотной характеристики, получаете их посредством следующей операции:

$H_{i}=H(i*\Delta f)$. Здесь $H(\circ)$ - Ваша, аналитически заданная, частотная характеристика.

Задаёте также $T$ - исследуемый временной период (на котором однозначно все переходные процессы закончатся), его значение обычно выбирается исходя из "физики" системы.

Задаёте также $Nt$ - количество отсчётов Вашего временного ряда (во временной области).

Откуда следует

$\Delta t=\frac{T}{Nt}$,

из чего по теореме Котельникова получаем:

$\max\{f\}=\frac{1}{2\Delta t}$.

Далее рассчитываем:

$Nf=2Nt$, где $Nf$ - количество отсчётов Вашего временного ряда (в частотной области).

И окончательно:

$\Delta f=\frac{\max\{f\}}{Nf}$.

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

P.S.
А вообще держите перед глазами теорему Котельникова об отсчётах, тогда Вы сможете самостоятельно перевести преобразования из непрерывной интегральной формы в дискретную, заданную конечными суммами. И конечно не ошибитесь с индексами.

Успехов!

[Altemirx]

G^a, Вам большое спасибо за подробные разъяснения! К сожалению, моя рабочая неделя закончилась, но в самое ближайшее время обязательно проработаю Ваше решение и сообщу о результатах.

Кстати, а насчёт результатов, полученных мной в MC11 и в Simulink никаких замечаний не будет?

С уважением,
Altemirx

[G^a]

Здравствуйте!

Altemirx писал(а):

Кстати, а насчёт результатов, полученных мной в MC11 и в Simulink никаких замечаний не будет?

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

[Altemirx]

Здравствуйте, уважаемый G^a.

Сегодня попробовал по приведённой методике произвести вычисления. Всё более-менее нормально до момента вычисления свёртки. Итак, по-порядку.

1. В приведённых Вами формулах сумма берётся от $N-1$ элементов. Как я понимаю, для временнОй области (и для последующего проведения свёртки) $N=Nt$? Или я неправ?
2. $N$ в формулах для $h_{k}$ и $y_{k}$ - одни и те же?
3. Насколько я знаю, импульсная характеристика сигнала есть $h(t)$ и представлена во временнОй области. Мы получаем $h_{k}$ также во временнОй области? Тогда как учитывать её комплексность? Или для представления на оси времени я должен брать её только по модулю? Правильно ли будет, если для оси времени я в этом случае буду значения времени брать как $\frac{k\cdot T}{Nt}$?
4. Если предположения в п.3 правильны, то $k_{max}$ получаем, исходя из последнего выражения?
5. Если предположения в п.1 и п.3 правильны, то свёртка берётся для временнОй области. Возникают неясности с представлением функции Хэвисайда: если для MC11 её брать как $\Phi(\frac{k\cdot T}{Nt}-i\cdot\Delta t)$, то результат будет некорректен, поскольку он критично зависит от $T, N, Nt$. Если функцию Хэвисайда определять средствами программирования, то возникает неопределённость для $x_{k-i}$, поскольку если $k<i$, то получаются отрицательные индексы. Как в этом случае правильно записать эту функцию?

Файл для MC11: http://forum.exponenta.ru/download.php?id=1720

С уважением,
Altemirx

[Altemirx]

Господа, ау?!

[G^a]

Добрый день!

Altemirx писал(а):

1. В приведённых Вами формулах сумма берётся от $N-1$ элементов. Как я понимаю, для временнОй области (и для последующего проведения свёртки) $N=Nt$? Или я неправ?

Да, для временной области $N=Nt$, но суммируются $N$ элементов, так как есть ещё нулевое значение индекса.

Altemirx писал(а):

2. $N$ в формулах для $h_{k}$ и $y_{k}$ - одни и те же?

Да.

Altemirx писал(а):

3. Насколько я знаю, импульсная характеристика сигнала есть $h(t)$ и представлена во временнОй области. Мы получаем $h_{k}$ также во временнОй области? Тогда как учитывать её комплексность? Или для представления на оси времени я должен брать её только по модулю?

Да, импульсная переходная функция - это скалярная вещественная функция скалярного вещественного аргумента; её физический смысл - это "взвешивание" отсчётов из прошлого, то есть она определяет, как прошлое системы, влияет на её настоящее.

Обратное преобразование Фурье от Вашей комплексной частотной характеристики, должно давать вещественную функцию.

Altemirx писал(а):

Правильно ли будет, если для оси времени я в этом случае буду значения времени брать как $\frac{k\cdot T}{Nt}$?

Считайте их так: $t_{k}=k\Delta t$.

Altemirx писал(а):

5.свёртка берётся для временнОй области.

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

Altemirx писал(а):

Возникают неясности с представлением функции Хэвисайда: если для MC11 её брать как $\Phi(\frac{k\cdot T}{Nt}-i\cdot\Delta t)$, то результат будет некорректен, поскольку он критично зависит от $T, N, Nt$. Если функцию Хэвисайда определять средствами программирования, то возникает неопределённость для $x_{k-i}$, поскольку если $k<i$, то получаются отрицательные индексы. Как в этом случае правильно записать эту функцию?

Зная $\Delta t$, Вы можете определить функию Хевисайда посредством выражения: $x_{k}=\Phi(k-k_{*})$, где: $k_{*}$ - момент "включения" функции; \Phi(\circ)$ - можете определить через обычный "if(...)".

Вы правы, в том, что результат свёртки зависит от $T, Nt$, но я Вам в самом начале говорил, что Вы получите оценку, а не точное значение выходного сигнала, и подход этот годится не для всех случаев, а только для большинства. Если Вы выберите достаточно большими значения параметров $T, Nt$, это будет означать, что величина $\Delta t$ будет маленькой, и тогда фронт переключения функции Хевисайда будет очень крутым и Вы получите достаточно точную инженерную оценку. На практике могу сказать, если $\frac{1}{2\Delta t}$ больше $20 f_g$, где f_g$ - граничная частота Вашей системы, то фактически функция Хевисайда определённая вышеуказанным образом будет "прямоугольной".

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

Успехов!

[Altemirx]

Цитата:

Да, импульсная переходная функция - это скалярная вещественная функция скалярного вещественного аргумента; её физический смысл - это "взвешивание" отсчётов из прошлого, то есть она определяет, как прошлое системы, влияет на её настоящее.

Обратное преобразование Фурье от Вашей комплексной частотной характеристики, должно давать вещественную функцию.

В таком случае, я не совсем понимаю, на каком этапе выполнять обратное преобразование Фурье (если возвращаться к предложенному Вами алгоритму). Между п.1 и п.2? Необходимо тогда обратное дискретное Фурье-преобразование. А не будут ли перемешаны при этом частотная/временнАя области?

Спасибо ещё раз за конструктивные и терпеливые ответы.

[G^a]

Добрый день!

Altemirx писал(а):

В таком случае, я не совсем понимаю, на каком этапе выполнять обратное преобразование Фурье (если возвращаться к предложенному Вами алгоритму). Между п.1 и п.2?

Обратное преобразование Фурье можно проводить, как только Вы зададите $T, \Delta T. Оно необходимо, чтобы из частотной характеристики (полной: АЧХ + ФЧХ; записанной в комплексной форме) получить импульсную переходную функцию.

Altemirx писал(а):

Необходимо тогда обратное дискретное Фурье-преобразование. А не будут ли перемешаны при этом частотная/временнАя области?

Да, дискретное обратное преобразование Фурье. Именно его формулу я Вам и писал в самом начале. Я не совсем понял, что Вы имеет в виду под перемещиванием, но рискну предположить, что Вы говорите об Найквист-эффекте наложения частот, вызванном конечным разрешением по времени Вашего временного ряда. Но именно для подавления этого негативного эффекта, мы с Вами и считали величину: $max\{f\}=\frac{1}{2\Delta t}.

Altemirx писал(а):

Спасибо ещё раз за конструктивные и терпеливые ответы.

Пожалуйста! Но Вам бы стоило книги почитать.

[uni]

У меня такой вопрос к автору поста: "Вы математик или радиотехник?". Все эти тут математические выкладки только запутывают суть дела. Ничего там сложного нету во всех этих преобразованиях. Смущает нестандартный вид выражения для фильтра. Сложная там не сложная, чёрт с ней, это не проблема, была бы запись понятная. Я Вам вроде бы дал совет обратиться на telesys.ru, на формум по DSP, там есть спецы по цифровым фильтрам. В общем и целом, радиотехник мыслит следующим образом:

1) Делаем декомпозицию выражения, описывающего фильтр. Это означает, что нужно нарисовать структурную схему фильтра. Ваше выражение не является просто рациональной дробью, тогда было бы проще. Мы должны привести к виду, который поможет быстрее прийти к решению. Как делать декомпозицию. Если выражения перемножаются, то отдельные фильтры, отвечающие за каждое выражение, стоят последовательно, как например, выражение в скобках в пятой степени - это пять фильтров расположенных последовательно.
2) Если мы хотим вообще найти какую-нить реакцию, то проще всего это сделать численно. Для этого H(f) дискретизируют с помощью z-преобразования. Для этого и делается пункт 1), чтобы потом по таблицам найти H(z) для каждого элементарного фильтра. Получив H(z), можно тоже по таблицам найти и h[n] - это импульсная характеристика работающая по дискретному набору данных, она включает в себя время dt - период дискретизации.

Ваша проблема в том, что на практике никто никогда не видел такой формы записи частотной характеристики фильтра. Это больше похоже на задание "от балды" какого-то выражения. Запись с экспонентой - это ещё куда ни шло. Z-преобразование само по себе не сложное (есть ещё варианты дискретизации), когда рука набита, но я уже плохо помню, может гляну на досуге. Обычно на практике, той практике которой занимаются спецы по цифровой обработке сигналов, задание на разработку фильтра стандартизовано в некотором смысле и никто не придумывает функцию H(f). Ищут сразу коэффициенты фильтра для работы в непосредственно в железе. На что есть спец программы в том же Matlab'е - Filter Desing Tools. Опять же Ваш случай заставляет возвращаться к математическим истокам, что практикам не всем под силу (взять z-преобразование или методом Ципкина-Гольденберга, или Раггазини-Бергена привести систему к дискретному виду), а теретики мало чего смыслят в практике. Я хотел было Вам помочь, но книжку где-то в библиотеке своей затерял, где таблицы есть с уже сделанными преобразованиями. Насколько я помню они также есть в книжке Быкова "Цифровые методы в радиотехнике" - так вроде называется, если мне память не изменяет.
Да, преобразование Лапласа - это весчь. Но, мне кажется, сложновата Ваша функция для нахождения нужного интеграла. Опять же его используют, когда можно разложить H(f) на элементарные дроби, для которых уже известны преобразования - это та же декомпозиция 1). Всё упирается в нестандартный вид. Хотя с точки зрения математика, я думаю, без разницы какой там вид. Ах, да, я забыл, если Ваша H(f) является комплексной, а по-видимому так оно и есть, то Вы вообще от реальности далеко уходите, либо её существенно усложняете. Комплексный фильтр работает с комплексным входным сигналом. Такие сигналы применяются при квадратурной обработке. Например, в современных цифровых системах связи. Там есть или могут быть комплексные фильтры. Так вот, такой фильтр имеет 2 входа и 2 выхода. Когда Вы подаёте на один вход ступеньку, то на другой вход подаётся нолик, это неявно будет присутствовать при математическом описании фильтрации. Но выход может быть комплексным, тогда реакция будет иметь 2 составляющие реальную и мнимую. Определить комплексность филльтра в вашем случае просто: найдите H(f) при какой-то фиксированной частоте. Подставив, мы обнаруживаем, что фильтр комплексный.
И откуда Вы только взяли такое задание? Мне лично кажется оно сложноватым, только посмотрев на вид H(f) - отпадает всякое желание возиться. Дайте хоть пример похожего решения такой же задачи. Если оставить только скобки с пятой степенью, то задача упроститься. Лучше даже начинать с выражения в скобках. Если я не ошибаюсь, то выражение в скобках описывает дискретный фильтр, выход которого получается путём вычета из текущего отсчёта одного запаздывающего. (1 - exp(-z)) - это в дискретном варианте (если опять же не ошибаюсь). 1-текущий отсчёт, "минус" параллельное соединение фильтров, exp(-z) - запаздыващий отсчёт.

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

P.S. Про комплексный фильтр я тут загнул. |H(f)| у реального фильтра симметричная относительно f=0. Если нет, то фильтр комплексный. У Вас как раз симметричный вариант. Так что приношу заранее извинения, чёрт попутал, комплексность H(f) отвечает как раз за измение амплитуды и фазовую задержку. Симметричность берётся, насколько я помню, из физической реализуемости фильтра, т.е. h(t) при t < 0 должна быть нулевой, иначе сигналы из будущего смогут воздействовать на текущие, что считается невозможным. Когда я писал про 2 входа, я имел в виду программную\железную - цифровую реализацию, т.е. программной нет "земли", а в железной аналоговой части будет 2 входа + "земля". Обычный фильтр на операх - 4-х полюсник - 2 входа и 2 выхода, но это другое, здесь с учётом "земли". Т.е. тут не нужно путать цифровые и аналоговые фильтры.