Вычислите и постройте оценку передаточной функции между двумя последовательностями, x и y. Последовательность x состоит из белого Гауссова шума. y результаты фильтрации x с 30-м порядком фильтр lowpass с нормированной частотой среза $$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)))

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

Одномерное дискретное время колеблющаяся система состоит из модульной массы, $\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)

Оцените передаточную функцию системы как функция частоты. Используйте точки ДПФ 20:48 и задайте окно Кайзера с масштабным фактором 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' опция, чтобы произвести все четыре передаточных функции. Используйте окно Hann с 5000 выборками, чтобы разделить сигналы на сегменты. Задайте 2 500 выборок перекрытия между смежными сегментами и $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