tfestimate

Оценка передаточной функции

Описание

txy = tfestimate(x,y) находит оценку передаточной функции, txy, учитывая входной сигнал, x, и выходной сигнал, y.

  • Если x и y оба векторы, у них должна быть та же длина.

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

  • Если x и y матрицы с одинаковым числом строк, но различными количествами столбцов, затем txy мультивход/мультивыход (MIMO) передаточная функция, которая комбинирует все сигналы ввода и вывода. txy 3D массив. Если x имеет столбцы m и y имеет столбцы n, затем txy имеет столбцы n и страницы m. Смотрите Передаточную функцию для получения дополнительной информации.

  • Если x и y матрицы равного размера, затем tfestimate действует по столбцам: txy(:,n) = tfestimate(x(:,n),y(:,n)). Чтобы получить оценку MIMO, добавьте 'mimo' к списку аргументов.

пример

txy = tfestimate(x,y,window) использование window разделить x и y в сегменты и выполняют работу с окнами.

txy = tfestimate(x,y,window,noverlap) использование noverlap выборки перекрытия между смежными сегментами.

txy = tfestimate(x,y,window,noverlap,nfft) использование nfft выборка точек, чтобы вычислить дискретное преобразование Фурье.

txy = tfestimate(___,'mimo') вычисляет передаточную функцию MIMO для матричных входных параметров. Этот синтаксис может включать любую комбинацию входных параметров от предыдущих синтаксисов.

[txy,w] = tfestimate(___) возвращает вектор из нормированных частот, w, в котором оценивается передаточная функция.

пример

[txy,f] = tfestimate(___,fs) возвращает вектор из частот, f, описанный в терминах частоты дискретизации, fs, в котором оценивается передаточная функция. fs должен быть шестой числовой вход к tfestimate. Чтобы ввести частоту дискретизации и все еще использовать значения по умолчанию предыдущих дополнительных аргументов, задайте эти аргументы как пустой [].

[txy,w] = tfestimate(x,y,window,noverlap,w) возвращает оценку передаточной функции на нормированных частотах, заданных в w.

[txy,f] = tfestimate(x,y,window,noverlap,f,fs) возвращает оценку передаточной функции на частотах, заданных в f.

[___] = tfestimate(x,y,___,freqrange) возвращает оценку передаточной функции по частотному диапазону, заданному freqrange. Допустимые опции для freqrange 'onesided', 'twosided', и 'centered'.

пример

[___] = tfestimate(___,'Estimator',est) оценочные передаточные функции с помощью средства оценки est. Допустимые опции для est 'H1' и 'H2'.

tfestimate(___) без выходных аргументов строит оценку передаточной функции в окне текущей фигуры.

Примеры

свернуть все

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

Figure contains an axes. The axes with title Transfer Function Estimate via Welch contains an object of type line.

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

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

Figure Filter Visualization Tool - Magnitude Response (dB) contains an axes and other objects of type uitoolbar, uimenu. The axes with title Magnitude Response (dB) contains an object of type line.

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

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

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

Figure contains an axes. The axes contains an object of type line.

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

Одномерное дискретное время колеблющаяся система состоит из модульной массы, m, присоединенный к стенке к пружине модуля эластичная константа. Датчик производит ускорение, a, из массы в Fs=1 Гц. Демпфер препятствует движению массы путем порождения на него силы, пропорциональной, чтобы ускориться с затуханием постоянного b=0.01.

Сгенерируйте в 2000 раз выборки. Задайте интервал выборки Δt=1/Fs.

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

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

x(k+1)=Ax(k)+Bu(k),y(k)=Cx(k)+Du(k),

где x=[rv]T вектор состояния, r и v соответственно положение и скорость массы, u движущая сила, и y=a измеренный выход. Матрицы пространства состояний

A=exp(AcΔt),B=Ac-1(A-I)Bc,C=[-1-b],D=1,

I 2×2 идентичность и матрицы пространства состояний непрерывного времени

Ac=[01-1-b],Bc=[01].

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)

Figure contains an axes. The axes contains an object of type line.

Оцените передаточную функцию системы в зависимости от частоты. Используйте точки ДПФ 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])

Figure contains an axes. The axes contains 2 objects of type line.

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

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

Figure contains an axes. The axes with title Transfer Function Estimate via Welch contains an object of type line.

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

Идеальная одномерная колеблющаяся система состоит из двух масс, m1 и m2, ограниченный между двумя стенками. Модули таковы что m1=1 и m2=μ. Каждая масса присоединена к самой близкой стенке к пружине с эластичной константой k. Идентичная пружина соединяет эти две массы. Три демпфера препятствуют движению масс путем проявления на них, обеспечивает пропорциональный, чтобы ускориться, с затуханием постоянного b. Выборка датчиков a1 и a2, ускорения масс, в Fs=50 Гц.

Сгенерируйте в 30000 раз выборки, эквивалентные 600 секундам. Задайте интервал выборки Δt=1/Fs.

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

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

x(k+1)=Ax(k)+Bu(k),y(k)=Cx(k)+Du(k),

где x=[r1v1r2v2]T вектор состояния, ri и vi соответственно местоположение и скорость iмасса th, u=[u1u2]T вектор из входных движущих сил, и y=[a1a2]T выходной вектор. Матрицы пространства состояний

A=exp(AcΔt),B=Ac-1(A-I)Bc,C=[-2k-2bkbk/μb/μ-2k/μ-2b/μ],D=[1001/μ],

I 4×4 идентичность и матрицы пространства состояний непрерывного времени

Ac=[0100-2k-2bkb0001k/μb/μ-2k/μ-2b/μ],Bc=[00100001/μ].

Набор k=400, b=0, и μ=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 выборок перекрытия между смежными сегментами и 214 Точки ДПФ. Постройте оценки.

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

Figure contains 4 axes. Axes 1 with title Input 1, Output 1 contains 2 objects of type line. Axes 2 with title Input 2, Output 1 contains 2 objects of type line. Axes 3 with title Input 1, Output 2 contains 2 objects of type line. Axes 4 with title Input 2, Output 2 contains 2 objects of type line.

Передаточные функции имеют максимумы в ожидаемых значениях, ω1,2/2π, где ω собственные значения модальной матрицы.

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

    3.8470
   14.4259

Добавьте затухание в систему путем установки b=0.1. Вычислите эволюцию времени ослабленной системы с теми же движущими силами. Вычислите 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')

Figure contains an axes. The axes with title Transfer Function Estimate via Welch contains 4 objects of type line. These objects represent 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

Figure contains an axes. The axes contains 4 objects of type line. These objects represent I1, O1, I1, O2, I2, O1, I2, O2.

Входные параметры

свернуть все

Входной сигнал в виде вектора или матрицы.

Пример: cos(pi/4*(0:159))+randn(1,160) задает синусоиду, встроенную в белый Гауссов шум.

Типы данных: single | double
Поддержка комплексного числа: Да

Выходной сигнал в виде вектора или матрицы.

Типы данных: single | double
Поддержка комплексного числа: Да

Окно в виде целого числа или как строка или вектор-столбец. Используйте window разделить сигнал на сегменты.

  • Если window целое число, затем tfestimate делит x и y в сегменты длины window и окна каждый сегмент с Окном Хэмминга той длины.

  • Если window вектор, затем tfestimate делит x и y в сегменты той же длины как вектор и окна каждый сегмент с помощью window.

Если длина x и y не может быть разделен точно в целое число сегментов с noverlap перекрывающиеся выборки, затем сигналы являются усеченными соответственно.

Если вы задаете window как пустой, затем tfestimate использует Окно Хэмминга, таким образом что x и y разделены на восемь сегментов с noverlap перекрывающиеся выборки.

Для списка доступных окон смотрите Windows.

Пример: hann(N+1) и (1-cos(2*pi*(0:N)'/N))/2 оба задают окно Hann длины N  + 1.

Типы данных: single | double

Количество перекрытых выборок в виде положительного целого числа.

  • Если window скаляр, затем noverlap должно быть меньшим, чем window.

  • Если window вектор, затем noverlap должно быть меньшим, чем длина window.

Если вы задаете noverlap как пустой, затем tfestimate использует номер, который производит 50%-е перекрытие между сегментами. Если длина сегмента не задана, функция устанавливает noverlap к ⌊N/4.5 ⌋, где N является длиной сигналов ввода и вывода.

Типы данных: double | single

Количество ДПФ указывает в виде положительного целого числа. Если вы задаете nfft как пустой, затем tfestimate наборы этот аргумент к макс. (256,2p), где p = ⌈log2  N  для входных сигналов длины N.

Типы данных: single | double

Частота дискретизации в виде положительной скалярной величины. Частота дискретизации является количеством выборок в единицу времени. Если модуль времени является секундами, то частота дискретизации имеет модули Гц.

Нормированные частоты в виде строки или вектор-столбца по крайней мере с двумя элементами. Нормированные частоты находятся в рад/отсчете.

Пример: w = [pi/4 pi/2]

Типы данных: double

Частоты в виде строки или вектор-столбца по крайней мере с двумя элементами. Частоты находятся в циклах в единицу времени. Единица времени задана частотой дискретизации, fs. Если fs имеет модули выборок/секунда, затем f имеет модули Гц.

Пример: fs = 1000; f = [100 200]

Типы данных: double

Частотный диапазон для передаточной функции оценивает в виде той 'onesided', 'twosided', или 'centered'. Значением по умолчанию является 'onesided' для сигналов с действительным знаком и 'twosided' для сигналов с комплексным знаком.

  • 'onesided' — Возвращает одностороннюю оценку передаточной функции между двумя входными сигналами с действительным знаком, x и y. Если nfft является четным, txy имеет nfft/2 + 1 строка и вычисляется на интервале [0, π] рад/отсчет. Если nfft является нечетным, txy имеет (nfft + 1) строки/2 и интервал [0, π), рад/отсчет. Если вы задаете fs, соответствующие интервалы [0, fs/2] циклы/единица времени для даже nfft и [0, fs/2) циклы/единица времени для нечетного nfft.

  • 'twosided' — Возвращает двухстороннюю оценку передаточной функции между двумя входными сигналами с комплексным знаком или с действительным знаком, x и y. В этом случае, txy имеет nfft строки и вычисляются на интервале [0,2π), рад/отсчет. Если вы задаете fs, интервал [0, fs) циклы/единица времени.

  • 'centered' — Возвращает двухстороннюю оценку в центре передаточной функции между двумя входными сигналами с комплексным знаком или с действительным знаком, x и y. В этом случае, txy имеет nfft строки и вычисляются на интервале (–π, π] рад/отсчет для даже nfft и (–π, π) рад/отсчет для нечетного nfft. Если вы задаете fs, соответствующие интервалы (–fs/2, fs/2] циклы/единица времени для даже nfft и (–fs/2, fs/2) циклы/единица времени для нечетного nfft.

Средство оценки передаточной функции в виде 'H1' или 'H2'.

  • Используйте 'H1' когда шум является некоррелированым с входными сигналами.

  • Используйте 'H2' когда шум является некоррелированым с выходными сигналами. В этом случае количество входных сигналов должно равняться количеству выходных сигналов.

Смотрите Передаточную функцию для получения дополнительной информации.

Выходные аргументы

свернуть все

Оценка передаточной функции, возвращенная как вектор, матрица или 3D массив.

Нормированные частоты, возвращенные как вектор-столбец с действительным знаком.

Циклические частоты, возвращенные как вектор-столбец с действительным знаком.

Больше о

свернуть все

Передаточная функция

Отношение между входом x и выход y моделируется линейной, независимой от времени передаточной функцией txy. В частотном диапазоне, Y (f) = H (f) X (f).

  • Для single-input/single-output системы H 1 оценкой передаточной функции дают

    H1(f)=Pyx(f)Pxx(f),

    где Pyx является взаимной спектральной плотностью мощности x и y, и Pxx является степенью спектральная плотность x. Эта оценка принимает, что шум не коррелируется с системным входом.

    Для мультивхода/мультивыхода (MIMO) системы H 1 средство оценки становится

    H1(f)=PYX(f)PXX1(f)=[Py1x1(f)Py1x2(f)Py1xm(f)Py2x1(f)Py2x2(f)Py2xm(f)Pynx1(f)Pynx2(f)Pynxm(f)][Px1x1(f)Px1x2(f)Px1xm(f)Px2x1(f)Px2x2(f)Px2xm(f)Pxmx1(f)Pxmx2(f)Pxmxm(f)]1

    для входных параметров m и n выходные параметры, где:

    • Pyixk является взаимной спектральной плотностью мощности k th вход и i th выход.

    • Pxixk является взаимной спектральной плотностью мощности k th и i th входные параметры.

    Для двух входных параметров и двух выходных параметров, средство оценки является матрицей

    H1(f)=[Py1x1(f)Px2x2(f)Py1x2(f)Px2x1(f)Py1x2(f)Px1x1(f)Py1x1(f)Px1x2(f)Py2x1(f)Px2x2(f)Py2x2(f)Px2x1(f)Py2x2(f)Px1x1(f)Py2x1(f)Px1x2(f)]Px1x1(f)Px2x2(f)Px1x2(f)Px2x1(f).

  • Для single-input/single-output системы H 2 оценками передаточной функции дают

    H2(f)=Pyy(f)Pxy(f),

    где Pyy является степенью спектральная плотность y и Pxy = P *yx является сопряженным комплексным числом взаимной спектральной плотности мощности x и y. Эта оценка принимает, что шум не коррелируется с системой выход.

    Для систем MIMO H 2 средства оценки четко определены только для равных количеств вводов и выводов: n = m. Средство оценки становится

    H2(f)=PYY(f)PXY1(f)=[Py1y1(f)Py1y2(f)Py1yn(f)Py2y1(f)Py2y2(f)Py2yn(f)Pyny1(f)Pyny2(f)Pynyn(f)][Px1y1(f)Px1y2(f)Px1yn(f)Px2y1(f)Px2y2(f)Px2yn(f)Pxny1(f)Pxny2(f)Pxnyn(f)]1,

    где:

    • Pyiyk является взаимной спектральной плотностью мощности k th и i th выходные параметры.

    • Pxiyk является сопряженным комплексным числом взаимной спектральной плотности мощности i th вход и k th выход.

Алгоритмы

tfestimate усредненный метод периодограммы валлийцев использования. Смотрите pwelch для деталей.

Ссылки

[1] Vold, Håvard, Джон Кроули и Г. Томас Роклин. “Новые Способы Оценить Функции Частотной характеристики”. Звук и Вибрация. Издание 18, ноябрь 1984, стр 34–38.

Расширенные возможности

Смотрите также

| | |

Представлено до R2006a
Для просмотра документации необходимо авторизоваться на сайте