Здравствуйте. Данный вопрос возник в процессе выполнения контрольной по ЦОС, кучу литературы перекопала, но следующий вопрос так и остался непонятен. Уже совсем запуталась, а преподаватель толком не объясняет.

я записала передаточную функцию аналогового фильтра прототипа H(s), теперь мне надо из нее получить H(z) методом инвариантного преобразования импульсной характеристики.

функция записана в виде произведения 4 множителей (порядок получился N=8), причем в числителе и знаменателе многочлен второй степени. (конкретный H(s) приложен).

собственно вопрос: согласно заданному методу найденную H(s) нужно представить в виде суммы простых дробей. нужно ли мне из всего произведения делать разложение? или каждый множитель - это отдельный фильтр, тогда я беру первый множитель и раскладываю его на две простые дроби, далее нахожу h(t) - h(nT) - H1(z), аналогично H2(z),H3(z) и H4(z), а итоговую нахожу как их произведение H(z)=H1(z)H2(z)H3(z)H4(z), т.е. структурная схема цифрового фильтра будет как на рис.2
или это неправильно?

дело в том, что программу надо написать в маткаде, я пробовала H(s) представлять в виде суммы простых дробей, правда в знаменателе так и оставался многочлен второй степени. и возникли трудности на этапе h(nT) - H1(z) при использовании функции z-преобразования...

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

заранее огромное спасибо, очень надеюсь на наставление на правильный путь

Признаться, всё забыл.
Если s=сигма+j*омега, где омега - круговая частота, то для прямого перехода на Z - плоскость используют билинейное z - преобразование. Тогда и получаются рекурсивные фильтры, вроде нарисованного вами каскадного.
Но насколько помню, инвариантность импульсной характеристики подразумевает расчёт импульсной характеристики фильтра по известной частотной, затем расчёт выборок импульсной характеристики фильтра с заданной частотой дискретизации от момента времени 0 до момента затухания импульсной характеристики (у вас фильтр с бесконечной импульсной характеристикой, поэтому до какого-то очень малого уровня - определяется требуемой точностью), затем просто выборки импульсной характеристики берутся в качестве коэффициентов нерекурсивного (КИХ) фильтра.

Есть книжка Ричард Лайонс. Цифровая обработка сигналов. Там про это написано. В интернете была. И здесь есть. Если модератора упросите...

Структура результирующего фильтра в любом случае будет как на рис.
Что касается самого расчета:
Вся дробь разлагается на сумму простейших.
Дроби вида a/(s+ конвертятся в a/(1-exp(-b/Fd)*z^-1)
Дроби второго порядка разлагаются в суммы
a/(s++conj(a)/(s+conj()
Здесь a и b уже комплексные.
Чтобы не париться:
a=g+j*h
b=sigma+j*omega
Тогда
a/(s++conj(a)/(s+conj() = (2*g*s+2*(sigma*g+omega*h))/(s^2+2*sigma*s+(sigma^2+omega^2))
Отсюда можно найти sigma omega g h сразу из дроби 2-го порядка.
Затем получить уже цифровую дробь 2-го порядка:

(2*g - exp(-sigma/Fd)*(2*g*cos(omega/Fd)-2*h*sin(omega/Fd))*z^-1)/(1 - 2*exp(-sigma/Fd)*cos(omega/Fd)*z^-1 + exp(-2*sigma/Fd)*z^-2)

Затем полученные дроби приводятся к общему знаменателю, суммируются. Ищутся нули - полюсы.
Комплексно сопряженные пары нулей-полюсов дают звенья 2-го порядка, объединяемые последовательно.

Либо реализовывать фильтр параллельным включением полученных биквадратных блоков.

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


сумма получилась вот такого вида (при помощи маткада) (эскиз внизу сообщения)
а чтобы получить цифровую дробь - это ведь обратное преобразование Лапласа? то есть результат будет тот же, если к этому выражению я применю готовую функцию invlaplace,s ?
в моем выражение вроде не совсем так получилось, как Вы в общем виде записали...
нули-полюсы ищутся из полученной дроби, приведенной к общему знаменателю?

спасибо, Ваш ответ для меня был наиболее исчерпывающий


Согласно Айфичеру, получила H(s) затем, чтобы записать H(z) надо к полученной H(s) применить обратное преобразование Лапласа, но надо предварительно раскладывать на простые дроби, поэтому и возник вышеописанный вопрос

ну и каждое из слагаемое при помощи того же маткада разложилось на множители:

сейчас попробую сделать замену, которую предложили, результаты приложу

Вот вам шпаргалка. маткад 14

Обратное преобразование Лапласа, применённое к передаточной функции линейной системы H(s) (которую для фильтров принято называть частотной характеристикой, поскольку сигма полагают равной нулю и преобразование Лапласа становится эквивалентно преобразованию Фурье), даст функцию времени, называемую импульсной характеристикой аналогового фильтра - прототипа.
Импульсная характеристика фильтра - сигнал на выходе фильтра, на вход которого подали дельта - функцию (импульс с длительностью, стремящейся к нулю, и амплитудой, стремящейся к бесконечности).
Если рассчитать значения импульсной характеристики через интервалы времени, равные периоду частоты дискретизации, начиная с момента подачи дельта-функции на вход, то таким образом будут вычислены отсчёты импульсной характеристики искомого нерекурсивного цифрового фильтра.
Отсчёты импульсной характеристики являются коэффициентами искомого нерекурсивного цифрового фильтра.
Так обеспечивается эквивалентность импульсной характеристики аналогового фильтра - прототипа и цифрового фильтра. Поэтому и метод называется "инвариантность импульсной характеристики".
Если от Вас хотят, чтобы Вы его использовали - так используйте, он несложный.
Есть масса других методов синтеза цифровых фильтров, куда более популярных. Этот метод более естественен в задачах моделирования поведения динамических систем.

Все верно. За тем исключением, что метод инвариантности их - метод синтеза именно рекурсивных фильтров со всеми вытекающими отсюда последствиями.
Все не так просто, короче говоря...

спасибо. завтра покажу, что у меня получилось.


да, нужно использовать именно этот метод. другие - в других вариантах.
а вот сигма получилась отличной от нуля. (если вы говорим об одном и том же) - это ведь множитель в знаменателе при s ?
и еще - при использовании обратного преобразования Лапласа в маткаде появляется функция Dirac(t) , что с ней то делать?...

Должен извиниться. Наверное, thermit прав. Хотя алгоритм, о котором я Вам рассказал, работает, но Ваш преподаватель скорее всего имел ввиду другой, его в современных книжках по ЦОС приводят.
Обратное преобразование Лапласа не стоит делать в Mathcad. Дроби в изображении Лапласа имеют хорошо известный оригинал (посмотрите мат. справочник - или книжку по ЦОС с этим методом). А функции Дирака у Вас не должно было получиться, если изображения по Лапласу - абсолютно интегрируемые функции. Вообще, положительные действительные части полюсов H(s) соответствуют неустойчивым системам - крайне сомнительно, что у Вас этот случай. Проверьте ещё раз.

Наташ, то что вам тут пытаются советовать не имеет ничего общего с тем, что вам нужно сделать.
Сигма, про которую вам рассказывает Дмитрий_Б это другая сигма.


Вам нужно синтезировать цифровой рекурсивный фильтр по заданному аналоговому прототипу методом инвариантности импульсной характеристики. Внимательно посмотрите мою писанину в маткаде.
Там на вашем примере по шагам расписана методика такого синтеза.

наконец установился 14-й маткад (у меня была другая версия, этот файл не открывался) открыла, посмотрела.
ого!!! Вы такую огромную работу проделали, спасибо большущее!
ну и в этом свете еще один вопрос: как построить карту нулей полюсов синтезированного фильтра?
ФЧХ сама построила (ну это совсем просто )


а можете посмотреть этот файл и подсказать где тут ошибка?
ну вот, хотела маткадовский файл прикрепить, не получилоась, пришлось картинками...
это уже метод взвешивания. теоретически определили,что используем вырезающую функцию Кайзера, для нее записали импульсную характеристику, она даже построилась, и вроде как правильно (хотя здесь тоже есть вопрос - по формуле для ФНЧ импульсная характеристика hd(n)=2fc*sin(n*w среза)/(n*w среза), если подставлять именно круговую частоту, все, кроме первого составляющего равны нулю..., потому как появляется 2*пи и синус становится равным нулю, а вот если просто fc , т.е. частоту среза, то все нормально).
а вот дальше, когда нужно выполнить преобразование H(z)= sum(h(n)*z^-n) не получатеся.
попробовала вручную записать эту функцию, т.е. соответствующее значение h(n) умножила на z в соответствующей степени, но АЧХ какая-то непонятная получается...

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


возник вопрос: по Айфичеру перед тем, как применять метод инвариантного преобразования импульсной характеристики, следует масштабировать нормированную передаточную функцию. Для этого s меняется на s/a, где a=2∙π∙fp, в результате чего фильтр получит желаемую характеристику.
так вот, та характеристика, которую я привела изначально - не масштабированная, если в Вашем файле я ее масштабирую (а сделала я это предельно внимательно), ПФ превратилась в прямую...

вот такой график (масштабированное P1 и R1) так может не надо масштабировать? тогда зачем написано...?

Наталия, присмотрелся я повнимательнее к Вашей H(s). Где Вы её такую взяли? Это фильтр с неограниченной полосой пропускания. Отсюда и функция Дирака (дельта-функция) во временной области.
Посмотрите сами: при омега (s=i*омега, i - мнимая единица) стремящейся к бесконечности (частота стремится к бесконечности) H(s)=const.
Крайне сомнительный фильтр - прототип. Лучше перепроверьте условие задачи.

Фильтр - как фильтр. Ничего выдающегося в нем нет. Фильтр устойчивый. Чего еще надо? Все вполне решаемо. Не парьте мозг ни себе ни студенту. Это просто учебная задача.

фильтр чебышева II типа, исходные данные: Fd=12000 Гц fp =2300 Гц fs =3000 Гц as =3дБ ap =55дБ
H(s) вывела сама, по известным формулам, там трудно было ошибиться.

Thermit , а что по поводу масштабирования скажете? не надо его делать выходит?

Это не фильтр чебышева 2. Никаких резонансов этот фильтр иметь не должен. Где-то вы ошиблись.



Будет время - гляну. Сейчас некогда. Вообще надо приводить нормированный фильтр к нужной полосе, конечно.

Наталия, фильтр Чебышева - неудачный выбор. Именно из-за наличия дельта-функций в импульсной характеристике непрерывный ("аналоговый") фильтр Чебышева считается физически нереализуемым (цифровой фильтр Чебышева - реализуем). Так зачем же брать такой фильтр в качестве аналогового прототипа, да ещё при методе инвариантности импульсной характеристики?
К стати, разложить на простейшие дроби дробно-рациональную передаточную функцию, где числитель и знаменатель имеют равные порядки, как я понимаю, нельзя (в действительности это проявление той же самой трудности, что я описал выше).
Возьмите фильтр Баттерворта. Он идеально подходит на роль прототипа. Примеры решений по методу инвариантности импульсной характеристики для него в учебниках найти можно.

Даже не знаю, как прокомментировать это утверждение...
Т е фильтр с передаточной функцией H(s) = const реализовать нельзя? Ведь их такого фильтра дельта-функция дирака... Чудеса, однако.




Как я уже сказал в ваших расчетах ошибка. Чтобы получился чебышев-2 вам надо в знаменателе к-ты при первой степени s удвоить. Смотрите приложенный файл. Там фильтр приведен к частоте 300 гц. Сами задайте какую надо частоту и руками скопируйте к-ты из разложения дроби на простейшие.
Как построить карту нулей/полюсов в маткаде - не знаю. Вообще - по оси абсцисс - вещественная часть, ординат - мнимая... Маткад, честно говоря, кроме ненависти никакой любви не вызывает.

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


по поводу этого могу только показать как все делала, кстати, в знаменателе при первой степени s коэффициент на 2 умножался.
ну вот так:
все формулы брала отсюда: http://www.dsplib.ru/content/filters/ch5/ch5.html


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

вот только теперь такой вопрос: почему масштабировав относительно fp у меня получилось, что на этой величине фильтр уже не пропускает, т.е. это значение превратилось в fs? сейчас еще раз перепроверила, все верно, написано a=2*pi*fc, но в моем случае fc не задано, а задано fp, но это однозначно не fs...


смотрите какая интересная вещь получилась, когда я масштабировала относительно fs, для R(f) получилась требуемая характеристика (сам прототип), а для P(f) G(f) графики получились одинаковые и оба уползли наверх.
а при масштабировании относительно fp P(f) и G(f) практически совпадают с R(f). все таки как правильно?

по поводу масштабирования - посмотрела Л.Рабинера,Б.Гоулда "теория и практика ЦОС", Гадзиковского, Солонину"Основы ЦОС", там ничего про масштабирование нет. только в Айфичере нашла. но там написано относительно частоты среза на уровне 3дБ...

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

Не. Не все. Когда копировали к-ты из разложения забыли про знаки...

не, не забыла вот:

ой, один коэффициент не туда скопировала
а что по поводу относительно чего масштабировать и почему не совпадает с матлабом?

Все совпадает. Вот скрипт.


Другое дело, что соотношение частоты дискретизации и полосы фильтра маленькое, поэтому получается не очень хорошая х-ка.
Матлаб ругается. Попробуйте частоту дискретизации раз в 10 увеличить...

надо то по заданию... я делала так (по образу и подобию Айфичера)

Fs=12000; % Sampling frequency
Ap=3;
As=55;
wp=2300/6000;
ws=3000/6000;
[N,Wc]= cheb2ord(wp,ws,Ap,As,'s'); % Определить порядок фильтра
%
% Create an analogue filter
%
[b, a]=cheby2(N, As, Wc, 's'); % Determine filter coeffs
[z, p, k]=cheby2(N, As,Wc, 's'); % Determine poles and zeros
%
% Convert analogue filter into Discrete IIR filter
%
[bz, az]=impinvar(b, a, Fs); % Determine coeffs of IIR filter

subplot(4,1,1) % Plot magnitude freq. response
[H, f]=freqz(bz, az, 512, Fs);
plot(f, 20*log10(abs(H)))
xlabel('Частота (Гц)')
ylabel('Амплитудная характеристика (дБ)')
subplot(4,1,2) % Вывести на экран ФЧХ
phase=(angle(H));
plot(f,phase)
xlabel('Частота (Гц)')
ylabel('Фаза (радианы)')

subplot(4,1,3) % Вывести на экран диаграмму нулей и полюсов
zplane(bz, az)
zz=roots(bz); % Determine poles and zeros
pz=roots(az);

subplot(4,1,4) % Вывести на экран импульсную характеристику
impz(bz,az)


так выходит я изначально H(s) неверно определила? или я неправильно в матлабе сделала?

Так должно быть



Кстати, фильтр с такими параметрами имеет 10-й порядок а не восьмой.

Лучше для разбора взять импульсную характеристику фильтра Чебышёва (с АЧХ того именно типа, что здесь обсуждается), имеющую аналитическое выражение. В частном случае ФНЧ обратимся к весовой функции Дольфа-Чебышёва, рассматривая ее как импульсную характеристику.
Нетрудно видеть, что АЧХ такого фильтра оптимальна (это спектр весовой функции).
Теперь, как того требует метод инвариантности имп. х-ки, возьмём отсчёты представленной имп. х-ки: хотя бы один - в начале (t=-Т/2) и один в середине (t=0). Отсчёт в середине - какое-то положительное ограниченное число. Первый отсчёт - плюс бесконечность. Совершенно ясно, что такое соотношение между отсчётами не позволяет представить их в конечной разрядной сетке, равно как и реализовать в аналоговом виде - ввиду конечного значения динамического диапазона и полосы пропускания любого физического устройства.
Впрочем, всё это было известно с самого появления функции Дольфа - Чебышёва, исследованной как амплитудное распределения поля в антеннах.
Источник: Трухачёв А.А. Радиолокационные сигналы и их применения. -М.: Воениздат, 2005.

и правда, в матлабе 10-й, даже по тому, что у меня было... почему же у меня получился в маткаде 8-й? тут то вроде все формулы правильно подставила... а графики ведь совпали, когда Вы первый раз подставили в матлабе по полученной H(s) для 8-го порядка... с тем, что потом написали, где N=10

Не-а. С fp/fs 2300/3000 порядок 10. 8-ой 2100/3000.
Собственно матлаб все по тем же формулам считает.

я имела в виду, что когда Вы построили первый раз график, он же строился для H(s) c 8 порядком. а потом то Вы подставляли исходные данные, получился N=10, а график такой же.

да, и все-таки масштабируем относительно fp или fs? в первый раз у Вас f0=fs=3000

я уж стала пробовать баттерворта, там порядок 13 получился. а то мне еще ФВЧ на основе прототипа рассчитывать, а там исходные данные другие, и если здесь при разложении получилось, что в знаменателе s^2, то для моего ФВЧ получилось две дроби с простыми полюсами и одна с комплексно-сопряженными. там уже по Вашему шаблону сделать не получится.


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

Фильтр чебышева однозначно задается любыми 3-мя из 4-х параметров:
1 Неравномерность в пп
2 Ослабление в пз
3 Минимальная частота на которой есть требуемое ослабление в пз
4 Порядок фильтра



Параметр 3 можно задать как конец пп или как начало пз. Т е что считать 1? Wp или Ws?
Ну и нормирование отсюда и следует. Ваш фильтр расчитан с Ws=1. Поэтому денормировка через fs.

А можно и так:

а-а-а, понятно теперь.

с баттервортом тоже не все так просто... dirac тут конечно нет, но порядок почти в 2 раза больше. ну вобщем что-то так и опять же без H(z)

в маткаде сделала как в матлабе, N получилось 10. в разных источниках одним и тем же разные выражения обозначают, запуталась, поэтому неправильно получилось

я сейчас объединяю свое начало с Вашим файлом, с учетом, что N=10, попозже его Вам покажу, проверьте пожалуйста. (маткад часто подвисает, долго получается), постараюсь в течение ближайшего часа

Вы в корне ошибаетесь.
ИХ фильтра с аппроксимацией ЧХ дробями чебышева да и не только чебышева, а любыми другими рациональными дробями выражается в виде суммы произведений экспонент на гармонические функции. ДФ дирака появляется при равенстве порядков полиномов числителя и знаменателя как дополнительное слагаемое (см обратное преобразование лапласа рациональной дроби). То, что вы привели никакого отношения к типа нашему фильтру не имеет.



А почему вас не устраивает чебышев? ДФ в ИХ испугались? Зря. Ничего страшного в ней нет.

вопросик - там где k менялоась от 0 до 3 это ведь 8/2=4? т.е. если N=10, то k должно меняться до 4? а где было до 8 соответственно до 10?

нет, все-таки продожаю с Чебышевым, с баттервортом ради интереса попробовала


а при построении h(n) я просто эту составляющую выбросила (иначе график не могу построить) - это не будет грубой ошибкой?

ну вот так получилось. сейчас попробую для ФВЧ (там ведь придется еще преобразовывать ФНЧ - ФВЧ)...

Угу.



Будет. Замените ее на 0^t и постройте вторым наложенным графиком с каким-нить маркером. Так будет нормально.

зы
Вы учитесь? Где, если не секрет?

вопрос: Thermit, у Вас в программе используется переменная i и я ее определила как мнимую единицу - эти записи друг другу не мешают?


я уже закончила, сестренке помогаю с контрольной - зо самарской ак.связи - у них всего одна установочная лекция была и никакого методического материала кроме нескольких книг по ЦОС в электронном виде - как хочешь, так и делай.. вот уж пришлось так, зато много интересного для себя почерпнула

Вполне могут. Вообще-то в маткаде мнимая единица пишется как 1i.



Сложный предмет. Без пива не осилить.

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

Похоже на правду.

для ФВЧ получился нечетный порядок фильтра, там без Вашей помощи опять никак....

и почему то при масштабировании функции не раскладывается на простые дроби

Позже.

я кажется поняла: для нечетного порядка фильтра нужно будет в формулу для ПФ для цифрового эквивалента добавить одно слагаемое - 1/(1-z^-1 exp(-sigma0*T)), где sigma0 - это то, что было в скобке (s-sigma0) правильно? но реализовать это программно боюсь, не смогу. я попробую, а вдруг... но обязательно пришлю на проверку.

ну вот только так получилось.
не получается разложить функцию на простые дроби в случае масштабирования (ну и значит дальше дело остановилось)
не получается рассчитать последовательную форму (это я попробовала переделать формулу под нечетное N)
ну и не факт, что паралельную форму правильно подкорректировала (опять же под нечетное N)

разложить на простые дроби получилось, нашла причину. а вот остальное в силе

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

Параметры фильтра какие?

Вы про ФВЧ спрашиваете? в первой строке они прописаны

вот только поменяны местами fp и fs, потому что сначала делаю ФНЧ-прототип, потом его преобразую в ФВЧ.

Есть числовой пример: А. Оппенхейм, Р. Шафер. Цифровая обработка сигналов. Москва: Техносфера, 2007.
На стр.449-451: Пример 7.2 "Импульсная инвариантность и фильтр Баттерворта". Там, правда, фильтр 6 порядка.
Вообще, если найдёте книжку, обратите внимание на ограничения метода инвариантности (раздел 7.2.1).
- "Однако если характеристика (частотная) непрерывного фильтра становится практически нулевой при больших частотах, то наложение спектров может оказаться пренебрежимо малым, и в результате дискретизации импульсной характеристики непрерывного фильтра получится полезный дискретный фильтр."
-"...техника импульсной инвариантности работает только для фильтров с ограниченной полосой, например фильтров нижних частот или полосовых фильтров.
Попытка использовать данный метод для разработки режекторных фильтров или фильтров верхних частот требует дополнительных ограничений полосы (пропускания фильтра), предохраняющих от серьёзного искажения характеристик за счёт ложных частот."
Курсив мой.

Вы совершенно правы. Однако, в жизни всегда есть место пофигу...

книжку искать уже времени нет, до 9-го надо все сдать. а если я сейчас переключусь на другой фильтр, я точно не успею. тут уже ФНЧ сделан, осталось с ФВЧ до конца разобраться, да и с 17-м порядком не хочется, если честно признаться, возиться.

да, это знаю. но этот метод был задан, выбирала я только тип фильтра. ну соответственно это и будет написано в выводах.

Thermit, что скажете по поводу нечетного порядка? правильно ли я изменила формулу для параллельной схемы и где ошибка в последовательной?
с полюсами если получится - это уж совсем идеал будет, так хотя бы АЧХ построить. с ИХ вроде как проблем нет

Thermit, я нашла, что нужно поменять чтобы все графики построились. Вы мне пожалуйста только на такой вопрос ответьте - правильно ли я записала G(f) и P(f)?

я попробовала просто ФНЧ сделать с параметрами моего ФВЧ (чтобы проверить правильность записи P(f), G(f)) - неправильные графики получились. Значит, я все-таки их неправильно задала. Thermit, подскажите пожалуйста как правильно их записать?