Jump to content

    

Вычисление АЧХ по переходной функции. Возможно ли?

Вроде бы, теоретически, комплексный коэффициент передачи K(w) = jwG(w), где w - круговая частота, G(w) - преобразование Фурье от переходной функции. Так?

Но вот рассмотрим в Matlab-е переходную функцию RC цепочки - y(t). tau=RC.

fs = 10;
N = 1024;
tau = 10;
t = 0:1/fs:(N-1)/fs;
f = 0:fs/N:fs*(N-1)/N;
y= 1-exp(-t/tau);

Коэффициент передачи всем известен: Kr = 1./(1+1i*2*pi*f*tau);

Как я ни пытался прикладывать ДПФ к массиву y, никак нужная зависимость не получается. Лучшая попытка:

Y = fft(y);

K=2i*pi*(0:N-1)/N.*Y+ y(N-1);

Здесь предполагалось, что до прихода ступеньки на выходе 0, а последний элемент массива y(N-1) уже вышел на постоянное значение 1.

В начале abs(K) и abs(Kr) практически одинаковые, но довольно быстро расходятся. Почему?

Как правильно получить из y массив K с отсчетами комплексного коэффициента передачи?

Share this post


Link to post
Share on other sites

По моему, заменить преобразование Лапласа на Фурье нельзя, замена p на jw делается уже после преобразования для нахождения афчх.

Share this post


Link to post
Share on other sites

ЧХ - преобразование Фурье не от переходной, а от импульсной характеристики. Импульсная - производная от переходной.

Share this post


Link to post
Share on other sites

 

10 hours ago, Alex-lab said:

По моему, заменить преобразование Лапласа на Фурье нельзя, замена p на jw делается уже после преобразования для нахождения афчх.

Можно подробнее? Формула из книжки и там именно Фурье. Лаплас вроде здесь не нужен.

 

8 hours ago, V_G said:

ЧХ - преобразование Фурье не от переходной, а от импульсной характеристики. Импульсная - производная от переходной.

Так вроде в формуле K(w) = jwG(w) умножение на jw это учитывает? Нет?

Edited by GMavr

Share this post


Link to post
Share on other sites

@GMavr

Для аналогового случая у вас переходная характеристика g(t) = 1 - exp(-t/(RC)) -> импульсная характеристика h(t) = 1/(RC)* exp(-t/(RC)), t >= 0. Применяя преобразование Фурье (свойство затухающей экспоненты), получаем K(f) = 1/(1+j*2pi*f*RC).

В дискретной системе: https://www.dummies.com/education/science/science-engineering/working-across-domains-example-take-the-rc-low-pass-filter-to-the-z-domain/

В примере АЧХ аналогового и дискретного фильтров совпадают в соответствующей полосе.

Поскольку у вас коэффициенты импульсной характеристики в дискретной форме известны (2 штуки в числителе, знаменателя нет), то немного изменив код, сможете получить вашу АЧХ с помощью fft: https://ccrma.stanford.edu/~jos/filters/Practical_Frequency_Response_Analysis.html

Share this post


Link to post
Share on other sites

@Grizzly

Спасибо за ответ и за ссылки. Но вопрос состоял немного не в этом.

Есть массив значений переходной функции системы. Полученный экспериментально. Нужно по этим значениям оценить АЧХ системы. С экспериментальными данными ничего не выходит. Попробовал заполнить массив значениями переходной функции RC цепочки, для которой все хорошо известно. И тоже ничего не вышло.

Можно, конечно, по известным входным и выходным данным построить фильтр (Виннеровское оценивание). И потом вычислить АЧХ этого фильтра. Но нет ли пути попроще? Да и результат такого подхода меня как-то не впечатлил, хотя этот способ более или менее рабочий.

6 hours ago, Grizzly said:

Поскольку у вас коэффициенты импульсной характеристики в дискретной форме известны (2 штуки в числителе, знаменателя нет)

Вроде 1 коэффициент в числителе и 2 в знаменателе для цифрового аналога RC фильтра.

Share this post


Link to post
Share on other sites
1 час назад, GMavr сказал:

Вроде 1 коэффициент в числителе и 2 в знаменателе для цифрового аналога RC фильтра.

Прошу прощения. Конечно же, у вас БИХ фильтр.

По виду переходной характеристики вы можете оценить tau = RC, а затем рассчитать частотную характеристику. Или у вас произвольная система?

Как вариант - найти производную от переходной характеристики (diff), а затем посчитать FFT. Пример на стр. 14-16: https://ethz.ch/content/dam/ethz/special-interest/mavt/dynamic-systems-n-control/idsc-dam/Lectures/Signals-and-Systems/Lectures/Fall2018/Lecture11_sigsys.pdf

Но хорошие результаты всё равно не получатся в таком случае.

Share this post


Link to post
Share on other sites

Импульсная х-ка аналогового фильтра - реакция на воздействие дельта-функции дирака. Импульсная х-ка дискретной системы - реакция на воздействие дискретной дельта-функции. Разница в воздействиях, в принципе самих систем порождают когнитивный диссонанс. Ну и умножение на jw этот диссонанс только усугубляет. Не проще почитать про дискретизацию аналоговых фильтров? Многое прояснится.

Share this post


Link to post
Share on other sites
10 часов назад, GMavr сказал:

Есть массив значений переходной функции системы. Полученный экспериментально. Нужно по этим значениям оценить АЧХ системы.

Импульсная характеристика получается из переходной дифференцированием, т.е. подчеркивает (усиливает) верхние частоты. Если есть проблемы с адекватностью АЧХ, вычисленной по экспериментальной ПХ , значит, есть проблемы с точностью дифференцирования и/или с частотой Найквиста.

Попробуйте брать точки на переходной характеристике почаще, эти проблемы постепенно уйдут за точность вычислений

Share this post


Link to post
Share on other sites

Спасибо.

On 12/30/2019 at 1:10 AM, Grizzly said:

Как вариант - найти производную от переходной характеристики (diff), а затем посчитать FFT.

Спасибо. Пока нет шума, полученный данным способом график АЧХ почти совпадает с тем, что должно быть в идеале. Пробовал также численное интегрирование. Результат чуть лучше, но не принципиально. А вот добавление даже небольшого шума заметно портит картину.

fs = 10;
N = 1024;
tau = 10;
t = 0:1/fs:(N-1)/fs;
f = 0:fs/N:fs*(N-1)/N;
y = 1-exp(-t/tau); % Реакция RC цепи на ступеньку (переходная функция)
yn = y + (rand(size(t))-0.5)*0.01; % + шум
% Комплексный коэффициент передачи, вычисленный по переходной функции
K = fft(diff(y));
% Тоже + шум
Kn =  fft(diff(yn));
% Такой коэффициент передачи должен быть
Kr = 1./(1+1i*2*pi*f*tau);
plot(f(1:N/2), abs(K(1:N/2)), 'b',f(1:N/2), abs(Kn(1:N/2)), 'g', f(1:N/2), abs(Kr(1:N/2)), 'r');
legend('K', 'Kn', 'Kr');

afc.png.aeadc70463e91a55da5f02993a7f9a71.png

Тут, наверно, ничего сделать нельзя. Придется ограничится максимальной частотой 1/20 от частоты дискретизации.

Share this post


Link to post
Share on other sites
14 часов назад, GMavr сказал:

А вот добавление даже небольшого шума заметно портит картину.

Изначально у вас про шум ничего не говорилось. Это уже совершенно иная задача. Необходимо использовать различные методы идентификации систем. В принципе вы можете попытаться как-то сгладить полученную нынешнюю оценку.

Как вам посоветовал @V_G, остаётся только увеличивать частоту дискретизации, чтобы получать больше точек (чаще) переходной характеристики.

Share this post


Link to post
Share on other sites
18 часов назад, GMavr сказал:

Спасибо.

Спасибо. Пока нет шума, полученный данным способом график АЧХ почти совпадает с тем, что должно быть в идеале. Пробовал также численное интегрирование. Результат чуть лучше, но не принципиально. А вот добавление даже небольшого шума заметно портит картину.

 

Тут, наверно, ничего сделать нельзя. Придется ограничится максимальной частотой 1/20 от частоты дискретизации.

Преобразование одной функции в другую относится только к линейным системам. Добавление шума убивает линейность.

Share this post


Link to post
Share on other sites
On 1/6/2020 at 6:14 PM, Tanya said:

Преобразование одной функции в другую относится только к линейным системам. Добавление шума убивает линейность.

to TC.
Если шум аддитивный, (а у Вас "yn = y + (rand(size(t))...."), то все с ним боряться усреднением/накоплением. 
Запустите Вашу функцию: "yn = y + (rand(size(t))-0.5)*0.01; % + шум"  в цикле и будет Вам точность...
Кстати, а почему у Вас  "rand" , а не "randn" ?                                                                            
to Tanya: "Добавление шума убивает линейность." -- ну если ОН ну оооочень большой...

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
Sign in to follow this