Вычислите и постройте график оценки передаточной функции между двумя последовательностями, x и y. Последовательность x состоит из белого гауссова шума. y результаты фильтрации x с фильтром нижних частот 30-го порядка с нормированной частотой отсечения $$\mathrm{}\mathrm{}\mathrm{0,2\delta }$$ рад/выборка. Используйте прямоугольное окно для проектирования фильтра. Укажите частоту дискретизации 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)))

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

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

Создание 2000 временных выборок. Определите интервал выборки $\mathit{\Delta t}\mathrm{=}{\mathit{}}_{\mathit{}}$1/Fs.

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

Система может быть описана с помощью модели state-space

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

$\mathit{I}$ - $\mathrm{}\mathrm{}$тождество 2 × 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;

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

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 точек DFT и укажите окно Кайзера с коэффициентом формы 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{}}_{\mathrm{m1}}$ и ${\mathit{}}_{\mathrm{m2}}$, замкнутых между двумя стенками. Единицы измерения таковы, что ${\mathit{}}_{\mathrm{m1}}\mathrm{=}$ 1 ${\mathit{и}}_{\mathrm{}}\mathrm{m2}$= λ. Каждая масса крепится к ближайшей стенке пружиной с постоянной упругости $\mathit{}$k. Идентичная пружина соединяет две массы. Три демпфера препятствуют движению масс, воздействуя на них силами, пропорциональными скорости, с $\mathit{}$постоянной демпфирования b. Образцы датчиков ${\mathit{}}_{\mathrm{}}$a1 ${\mathit{}}_{\mathrm{}}$и a2, ускорения масс, ${\mathit{}}_{\mathit{}}при\mathrm{}$Fs = 50 Гц.

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

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

Система может быть описана с помощью модели state-space

где $$x{=\begin{array}{cccc}{}_{\mathrm{[}}& {}_{\mathrm{}}& {}_{\mathrm{}}& {}_{\mathrm{}}\end{array}}^{\mathrm{r1v1r2v2}}$$] T - $${\mathrm{вектор\; состояния,}}_{}$$ri $${}_{и}$$ vi - соответственно местоположение и скорость $$$$i-ой $$\mathrm{массы},{\begin{array}{cc}{u}_{\mathrm{}}& {=}_{\mathrm{}}\end{array}[}^{}$$u1u2] T - вектор входных движущих сил$$,{\begin{array}{cc}{}_{\mathrm{и}}& {}_{\mathrm{y}}\end{array}}^{=}$$ [a1a2] T - выходной вектор. Матрицы состояния-пространства:

$\mathit{I}$ - $\mathrm{}\mathrm{}$тождество 4 × 4, а матрицы состояния-пространства непрерывного времени -

Установите $\mathit{k}\mathrm{=}$400$\mathit{,}\mathrm{b}$= 0 $\mathrm{и}\mathrm{\lambda }$= 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];

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

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 образцов перекрытия между соседними сегментами и ${\mathrm{}}^{214}$ точками DFT. Постройте график оценок.

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

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

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

ans = 2×1

    3.8470
   14.4259

Добавьте демпфирование в систему, установив $\mathit{b}\mathrm{=}\mathrm{}$0,1. Вычислите эволюцию времени демпфированной системы с теми же движущими силами. Вычислите ${\mathit{}}_{\mathrm{}}$оценку H2 функции передачи 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