Вычислите и постройте график оценки передаточной функции между двумя последовательностями, x и y. Последовательность x состоит из белого Гауссова шума. y результаты фильтрации x с lowpass 30-го порядка с нормированной частотой отсечения $$0.2\pi$$ рад/образец. Для создания фильтра используйте прямоугольное окно. Задайте частоту дискретизации 500 Гц и окно Хэмминга длины 1024 для оценки передаточной функции.

h = fir1(30,0.2,rectwin(31));
x = randn(16384,1);
y = filter(h,1,x);

fs = 500;
tfestimate(x,y,1024,[],[],fs)

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

fvtool(h,1,'Fs',fs)

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

[Txy,f] = tfestimate(x,y,1024,[],[],fs);

plot(f,mag2db(abs(Txy)))

Оцените передаточную функцию для простой системы с одним входом/одним выходом и сравните ее с определением.

Одномерное дискретное время система состоит из модуля массы, $\mathit{m}$, прикрепленный к стене пружиной модуль упругой константы. Датчик дискретизирует ускорение, $\mathit{a}$, массы при ${\mathit{F}}_{\mathit{s}}=1$ Гц. Демпфер препятствует движению массы, прикладывая к ней силу, пропорциональную скорости, с постоянной демпфирования $\mathit{b}=0.01$.

Сгенерируйте 2000 временных выборок. Определите интервал дискретизации $\Delta \mathit{t}=1/{\mathit{F}}_{\mathit{s}}$.

Fs = 1;
dt = 1/Fs;
N = 2000;
t = dt*(0:N-1);
b = 0.01;

Система может быть описана моделью пространства состояний

где $\mathit{x}={\left[\begin{array}{cc}\mathit{r}& \mathit{v}\end{array}\right]}^{\mathit{T}}$ - вектор состояния, $\mathit{r}$ и $\mathit{v}$ являются соответственно положением и скоростью массы, $\mathit{u}$ является движущей силой, и $\mathit{y}=\mathit{a}$ - измеренный выход. Матрицы пространства состояний

$\mathit{I}$ является $2\times 2$ единичные, и матрицы пространства состояний в непрерывном времени

Ac = [0 1;-1 -b];
A = expm(Ac*dt);

Bc = [0;1];
B = Ac\(A-eye(size(A)))*Bc;

C = [-1 -b];
D = 1;

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

rng default

u = zeros(1,N);
u(1:N/2) = randn(1,N/2);

y = 0;
x = [0;0];
for k = 1:N
    y(k) = C*x + D*u(k);
    x = A*x + B*u(k);
end

plot(t,y)

Оцените передаточную функцию системы как функцию от частоты. Используйте 2048 точек ДПФ и задайте окно Кайзера с масштабным фактором 15. Используйте значение перекрытия по умолчанию между смежными сегментами.

nfs = 2048;
wind = kaiser(N,15);

[txy,ft] = tfestimate(u,y,wind,[],nfs,Fs);

Функция частотной характеристики системы дискретного времени может быть выражена как Z-преобразование передаточной функции временной области системы, оцененное по модулю кругу. Проверьте, что оценка вычислена tfestimate совпадает с этим определением.

[b,a] = ss2tf(A,B,C,D);

fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);
frf = polyval(b,z)./polyval(a,z);

plot(ft,20*log10(abs(txy)))
hold on
plot(fz,20*log10(abs(frf)))
hold off
grid
ylim([-60 40])

Постройте график оценки с помощью встроенной функциональности tfestimate.

tfestimate(u,y,wind,[],nfs,Fs)

Оцените передаточную функцию для простой системы с несколькими входами/несколькими выходами.

Идеальная одномерная колебательная система состоит из двух масс, ${\mathit{m}}_1$ и ${\mathit{m}}_2$, ограниченный между двумя стенками. Модули такие, что ${\mathit{m}}_1=1$ и ${\mathit{m}}_2=\mu$. Каждая масса присоединена к ближайшей стенке пружиной с упругой константой $\mathit{k}$. Идентичная пружина соединяет две массы. Три демпфера препятствуют движению масс, прикладывая к ним силы, пропорциональные скорости с постоянной демпфирования $\mathit{b}$. Выборка датчиков ${\mathit{a}}_1$ и ${\mathit{a}}_2$, ускорения масс, при ${\mathit{F}}_{\mathit{s}}=50$ Гц.

Сгенерируйте 30000 временных выборок, эквивалентных 600 секундам. Определите интервал дискретизации $\Delta \mathit{t}=1/{\mathit{F}}_{\mathit{s}}$.

Fs = 50;
dt = 1/Fs;
N = 30000;
t = dt*(0:N-1);

Система может быть описана моделью пространства состояний

где $$x={\left[\begin{array}{cccc}{r}_1& v_1& r_2& v_2\end{array}\right]}^T$$ - вектор состояния, $$r_i$$ и $$v_i$$ являются соответственно местоположением и скоростью $$i$$th-я масса, $$u={\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]}^T$$ является вектором входа движущих сил, и $$y={\left[\begin{array}{cc}{a}_1& a_2\end{array}\right]}^T$$ - вектор выхода. Матрицы пространства состояний

$\mathit{I}$ является $4\times 4$ единичные, и матрицы пространства состояний в непрерывном времени

Набор $\mathit{k}=400$, $\mathit{b}=0$, и $\mu =1/10$.

k = 400;
b = 0;
m = 1/10;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/m b/m -2*k/m -2*b/m];
A = expm(Ac*dt);
Bc = [0 0;1 0;0 0;0 1/m];
B = Ac\(A-eye(4))*Bc;
C = [-2*k -2*b k b;k/m b/m -2*k/m -2*b/m];
D = [1 0;0 1/m];

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

rng default
u = randn(2,N);

x = [0;0;0;0];
for kk = 1:N
    y(:,kk) = C*x + D*u(:,kk);
    x = A*x + B*u(:,kk);
end

Используйте входные и выходные данные, чтобы оценить передаточную функцию системы как функцию от частоты. Задайте 'mimo' опция для создания всех четырех передаточных функций. Используйте 5000-выборочное окно Ханна, чтобы разделить сигналы на сегменты. Задайте 2500 выборки перекрытия между смежными сегментами и $2^{14}$ Точки ДПФ. Постройте график оценок.

wind = hann(5000);
nov = 2500;

[q,fq] = tfestimate(u',y',wind,nov,2^14,Fs,'mimo');

Вычислите теоретическую передаточную функцию как Z-преобразование передаточной функции во временной области, оцененное в модуль круге.

nfs = 2^14;

fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

frf(1,:,1) = polyval(b1(1,:),z)./polyval(a1,z);
frf(1,:,2) = polyval(b1(2,:),z)./polyval(a1,z);

frf(2,:,1) = polyval(b2(1,:),z)./polyval(a2,z);
frf(2,:,2) = polyval(b2(2,:),z)./polyval(a2,z);

Постройте график теоретических передаточных функций и их соответствующих оценок.

for jk = 1:2
    for kj = 1:2
        subplot(2,2,2*(jk-1)+kj)
        plot(fq,20*log10(abs(q(:,jk,kj))))
        hold on
        plot(fz*Fs,20*log10(abs(frf(jk,:,kj))))
        hold off
        grid
        title(['Input ' int2str(kj) ', Output ' int2str(jk)])
        axis([0 Fs/2 -50 100])
    end
end

Передаточные функции имеют максимумы в ожидаемых значениях, ${\omega }_{1,2}/2\pi$, где $\omega$ являются собственными значениями модальной матрицы.

sqrt(eig(k*[2 -1;-1/m 2/m]))/(2*pi)

ans = 2×1

    3.8470
   14.4259

Добавьте демпфирование в систему путем установки $\mathit{b}=0.1$. Вычислите временную эволюцию демпфированной системы с теми же движущими силами. Вычислите ${\mathit{H}}_2$ оценка передаточной функции MIMO с использованием того же окна и перекрытия. Постройте график оценок с помощью tfestimate функциональность.

b = 0.1;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/m b/m -2*k/m -2*b/m];
A = expm(Ac*dt);
B = Ac\(A-eye(4))*Bc;
C = [-2*k -2*b k b;k/m b/m -2*k/m -2*b/m];

x = [0;0;0;0];
for kk = 1:N
    y(:,kk) = C*x + D*u(:,kk);
    x = A*x + B*u(:,kk);
end

clf
tfestimate(u',y',wind,nov,[],Fs,'mimo','Estimator','H2')
legend('I1, O1','I1, O2','I2, O1','I2, O2')

yl = ylim;

Сравните оценки с теоретическими предсказаниями.

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

frf(1,:,1) = polyval(b1(1,:),z)./polyval(a1,z);
frf(1,:,2) = polyval(b1(2,:),z)./polyval(a1,z);

frf(2,:,1) = polyval(b2(1,:),z)./polyval(a2,z);
frf(2,:,2) = polyval(b2(2,:),z)./polyval(a2,z);

plot(fz*Fs,20*log10(abs(reshape(permute(frf,[2 1 3]),[nfs/2 4]))))
legend('I1, O1','I1, O2','I2, O1','I2, O2')
ylim(yl)
grid