Main Content

La traducción de esta página aún no se ha actualizado a la versión más reciente. Haga clic aquí para ver la última versión en inglés.

tfestimate

Estimación de la función de transferencia

Descripción

txy = tfestimate(x,y) encuentra una estimación de la función de transferencia entre la señal de entrada x y la señal de salida y evaluada en un conjunto de frecuencias.

  • Si x e y son ambos vectores, deben tener la misma longitud.

  • Si una de las señales es una matriz y la otra es un vector, la longitud del vector debe ser igual al número de filas de la matriz. La función expande el vector y devuelve una matriz de estimaciones de la función de transferencia columna por columna.

  • Si x e y son matrices con el mismo número de filas pero con diferente número de columnas, txy es una función de transferencia multi-entrada y multi-salida (MIMO) que combina todas las señales de entrada y salida. txy es un arreglo tridimensional. Si x tiene m columnas e y tiene n columnas, txy tiene n columnas y m páginas. Para obtener más información, consulte Función de transferencia.

  • Si x e y son matrices del mismo tamaño, tfestimate opera por columnas: txy(:,n) = tfestimate(x(:,n),y(:,n)). Para obtener una estimación MIMO, añada 'mimo' a la lista de argumentos.

ejemplo

txy = tfestimate(x,y,window) utiliza window para dividir x e y en segmentos y aplicar ventanas.

txy = tfestimate(x,y,window,noverlap) utiliza muestras de solapamiento noverlap entre segmentos contiguos.

txy = tfestimate(x,y,window,noverlap,nfft) utiliza puntos de muestreo nfft para calcular la transformada discreta de Fourier.

txy = tfestimate(___,'mimo') calcula una función de transferencia MIMO para entradas matriciales. Esta sintaxis puede incluir cualquier combinación de argumentos de entrada de las sintaxis anteriores.

[txy,w] = tfestimate(___) devuelve un vector de frecuencias normalizadas, w, en las que se estima la función de transferencia.

ejemplo

[txy,f] = tfestimate(___,fs) devuelve un vector de frecuencias, f, expresado en términos de la tasa de muestreo, fs, a la que se estima la función de transferencia. fs debe ser la sexta entrada numérica a tfestimate. Para generar una tasa de muestreo y seguir utilizando los valores predeterminados de los argumentos opcionales anteriores, especifique estos argumentos como vacíos, [].

[txy,w] = tfestimate(x,y,window,noverlap,w) devuelve la estimación de la función de transferencia en las frecuencias normalizadas especificadas en w.

[txy,f] = tfestimate(x,y,window,noverlap,f,fs) devuelve la estimación de la función de transferencia en las frecuencias especificadas en f.

[___] = tfestimate(x,y,___,freqrange) devuelve la estimación de la función de transferencia a lo largo del rango de frecuencia especificado por freqrange. Las opciones válidas para freqrange son 'onesided', 'twosided' y 'centered'.

ejemplo

[___] = tfestimate(___,'Estimator',est) estima las funciones de transferencia con el estimador est. Las opciones válidas para est son 'H1' y 'H2'.

tfestimate(___) sin argumentos de salida representa la estimación de la función de transferencia en la ventana de la figura actual.

Ejemplos

contraer todo

Calcule y represente la estimación de la función de transferencia entre dos secuencias, x e y. La secuencia x consiste en ruido blanco gaussiano. y resulta de filtrar x con un filtro de paso bajo de 30.º orden con frecuencia de corte normalizada de 0.2π rad/muestra. Utilice una ventana rectangular para diseñar el filtro. Especifique una tasa de muestreo de 500 Hz y una ventana de Hamming de longitud 1024 para estimar la función de transferencia.

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 object. The axes object with title Transfer Function Estimate via Welch, xlabel Frequency (Hz), ylabel Magnitude (dB) contains an object of type line.

Utilice fvtool para verificar que la función de transferencia se aproxima a la respuesta en frecuencia del filtro.

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

Figure Figure 1: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Frequency (Hz), ylabel Magnitude (dB) contains an object of type line.

Obtenga el mismo resultado devolviendo la estimación de la función de transferencia en una variable y representando su valor absoluto en decibelios.

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

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

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

Estime la función de transferencia de un sistema simple de una entrada y una salida y compárela con la definición.

Un sistema oscilante unidimensional de tiempo discreto consiste en una masa unitaria, m, unida a una pared por un muelle de constante elástica unitaria. Un sensor muestrea la aceleración, a, de la masa a Fs=1 Hz. Un amortiguador impide el movimiento de la masa ejerciendo sobre ella una fuerza proporcional a la velocidad, con una constante de amortiguación b=0.01.

Genere 2000 muestras de tiempo. Defina el intervalo de muestreo Δt=1/Fs.

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

El sistema puede describirse mediante el modelo de espacio de estados

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

donde x=[rv]T es el vector de estado, r y v son respectivamente la posición y la velocidad de la masa, u es la fuerza motriz y y=a es la salida medida. Las matrices de espacio de estados son

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

I es la identidad 2×2 y las matrices de espacio de estados en tiempo continuo son

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;

Durante la mitad del intervalo de medición, una entrada aleatoria impulsa la masa. Utilice el modelo de espacio de estados para calcular la evolución temporal del sistema partiendo de un estado inicial totalmente nulo. Represente la aceleración de la masa como una función de tiempo.

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 object. The axes object contains an object of type line.

Estime la función de transferencia del sistema como una función de frecuencia. Utilice 2048 puntos de la DFT y especifique una ventana de Kaiser con un factor de forma de 15. Utilice el valor predeterminado de solapamiento entre segmentos contiguos.

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

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

La función frecuencia-respuesta de un sistema de tiempo discreto puede expresarse como la transformada Z de la función de transferencia del sistema en el dominio del tiempo, evaluada en el círculo unitario. Compruebe que la estimación calculada por tfestimate coincide con esta definición.

[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 object. The axes object contains 2 objects of type line.

Represente la estimación utilizando la funcionalidad incorporada de tfestimate.

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

Figure contains an axes object. The axes object with title Transfer Function Estimate via Welch, xlabel Frequency (mHz), ylabel Magnitude (dB) contains an object of type line.

Estime la función de transferencia de un sistema multi-entrada y multi-salida simple.

Un sistema ideal de oscilación unidimensional consiste en dos masas, m1 y m2, confinadas entre dos paredes. Las unidades son tales que m1=1 y m2=μ. Cada masa está unida a la pared más cercana por un muelle con una constante elástica k. Un muelle idéntico conecta las dos masas. Tres amortiguadores impiden el movimiento de las masas ejerciendo sobre ellas unas fuerzas proporcionales a la velocidad, con una constante de amortiguación b. Los sensores muestrean a1 y a2, las aceleraciones de las masas, a Fs=50 Hz.

Genere 30000 muestras de tiempo, equivalentes a 600 segundos. Defina el intervalo de muestreo Δt=1/Fs.

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

El sistema puede describirse mediante el modelo de espacio de estados

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

donde x=[r1v1r2v2]T es el vector de estado, ri y vi son respectivamente la ubicación y la velocidad de la i-ésima masa, u=[u1u2]T es el vector de las fuerzas motrices de entrada y y=[a1a2]T es el vector de salida. Las matrices de espacio de estados son

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

I es la identidad 4×4 y las matrices de espacio de estados en tiempo continuo son

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

Establezca k=400, b=0 y μ=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];

Una entrada aleatoria impulsa las masas a lo largo de la medición. Utilice el modelo de espacio de estados para calcular la evolución temporal del sistema partiendo de un estado inicial totalmente nulo.

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

Utilice los datos de entrada y salida para estimar la función de transferencia del sistema como una función de frecuencia. Especifique la opción 'mimo' para producir las cuatro funciones de transferencia. Utilice una ventana de Hann de 5000 muestras para dividir las señales en segmentos. Especifique 2500 muestras de solapamiento entre segmentos contiguos y 214 puntos de la DFT. Represente las estimaciones.

wind = hann(5000);
nov = 2500;

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

Calcule la función de transferencia teórica como la transformada Z de la función de transferencia en el dominio del tiempo, evaluada en el círculo unitario.

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);

Represente las funciones de transferencia teóricas y sus correspondientes estimaciones.

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 objects. Axes object 1 with title Input 1, Output 1 contains 2 objects of type line. Axes object 2 with title Input 2, Output 1 contains 2 objects of type line. Axes object 3 with title Input 1, Output 2 contains 2 objects of type line. Axes object 4 with title Input 2, Output 2 contains 2 objects of type line.

Las funciones de transferencia tienen máximos en los valores esperados, ω1,2/2π, donde ω son los valores propios de la matriz modal.

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

    3.8470
   14.4259

Añada amortiguación al sistema ajustando b=0.1. Calcule la evolución temporal del sistema amortiguado con las mismas fuerzas motrices. Calcule la estimación de H2 de la función de transferencia MIMO utilizando la misma ventana y solapamiento. Represente las estimaciones con la función 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 object. The axes object with title Transfer Function Estimate via Welch, xlabel Frequency (Hz), ylabel Magnitude (dB) contains 4 objects of type line. These objects represent I1, O1, I1, O2, I2, O1, I2, O2.

yl = ylim;

Compare las estimaciones con las predicciones teóricas.

[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 object. The axes object contains 4 objects of type line. These objects represent I1, O1, I1, O2, I2, O1, I2, O2.

Argumentos de entrada

contraer todo

Señal de entrada, especificada como vector o matriz.

Ejemplo: cos(pi/4*(0:159))+randn(1,160) especifica una sinusoide integrada en ruido blanco gaussiano.

Tipos de datos: single | double
Soporte de números complejos:

Señal de salida, especificada como vector o matriz.

Tipos de datos: single | double
Soporte de números complejos:

Ventana, especificada como un entero o como un vector fila o vector columna. Utilice window para dividir la señal en segmentos.

  • Si window es un entero, tfestimate divide x e y en segmentos de longitud window y aplica una ventana de Hamming de esa longitud a cada segmento.

  • Si window es un vector, tfestimate divide x e y en segmentos de la misma longitud que el vector y aplica ventanas a cada segmento con window.

Si la longitud de x e y no puede dividirse exactamente en un número entero de segmentos con muestras solapadas noverlap, las señales se truncan en consecuencia.

Si se especifica window como vacío, tfestimate utiliza una ventana de Hamming de forma que x e y se dividen en ocho segmentos con muestras solapadas noverlap.

Para obtener una lista de las ventanas disponibles, consulte Ventanas.

Ejemplo: hann(N+1) y (1-cos(2*pi*(0:N)'/N))/2 especifican ambas una ventana de Hann de longitud N + 1.

Tipos de datos: single | double

Número de muestras solapadas, especificado como entero positivo.

  • Si window es escalar, noverlap debe ser menor que window.

  • Si window es un vector, noverlap debe ser menos que la longitud de window.

Si especifica noverlap como vacío, tfestimate utiliza un número que genera un 50% de solapamiento entre segmentos.

Tipos de datos: double | single

Número de puntos DFT, especificado como entero positivo. Si se especifica nfft como vacío, tfestimate establece este argumento en max(256,2p), donde en las señales de entrada de longitud N p = ⌈log2 N y los símbolos ⌈ ⌉ determinan la función de techo.

Tipos de datos: single | double

Tasa de muestreo, especificada como un escalar positivo. La tasa de muestreo es el número de muestras por unidad de tiempo. Si la unidad de tiempo es el segundo, la tasa de muestreo tendrá unidades de Hz.

Frecuencias normalizadas, especificadas como vector fila o vector columna con al menos dos elementos. Las frecuencias normalizadas están en rad/muestra.

Ejemplo: w = [pi/4 pi/2]

Tipos de datos: double | single

Frecuencias, especificadas como vector fila o vector columna con al menos dos elementos. Las frecuencias están en ciclos por unidad de tiempo. La unidad de tiempo viene especificada por la tasa de muestreo, fs. Si las unidades de fs son muestras/segundo, las unidades de f son Hz.

Ejemplo: fs = 1000; f = [100 200]

Tipos de datos: double | single

Rango de frecuencia para la estimación de la función de transferencia, especificado como 'onesided', 'twosided' o 'centered'. El valor predeterminado es 'onesided' para señales de valor real y 'twosided' para señales de valor complejo.

  • 'onesided': Devuelve la estimación unilateral de la función de transferencia entre dos señales de entrada de valor real, x e y. Si nfft es par, txy tiene nfft/2 + 1 filas y se calcula en el intervalo [0,π] rad/muestra. Si nfft es impar, txy tiene (nfft + 1)/2 filas y el intervalo es [0,π) rad/muestra. Si se especifica fs, los intervalos correspondientes son [0,fs/2] ciclos/unidad de tiempo cuando nfft es par y [0,fs/2) ciclos/unidad de tiempo cuando nfft es impar.

  • 'twosided': Devuelve la estimación bilateral de la función de transferencia entre dos señales de entrada de valor real o complejo, x e y. En este caso, txy tiene nfft filas y se calcula en el intervalo [0,2π) rad/muestra. Si se especifica fs, el intervalo es [0, fs] ciclos/unidad de tiempo.

  • 'centered': Devuelve la estimación bilateral centrada de la función de transferencia entre dos señales de entrada de valor real o complejo, x e y. En este caso, txy tiene nfft filas y se calcula en el intervalo (–π,π] rad/muestra cuando nfft es par y (–π,π) rad/muestra cuando nfft es impar. Si se especifica fs, los intervalos correspondientes son (–fs/2, fs/2] ciclos/unidad de tiempo cuando nfft es par y (–fs/2, fs/2) ciclos/unidad de tiempo cuando nfft es impar.

Estimador de la función de transferencia, especificado como 'H1' o 'H2'.

  • Utilice 'H1' cuando el ruido no esté correlacionado con las señales de entrada.

  • Utilice 'H2' cuando el ruido no esté correlacionado con las señales de salida. En este caso, el número de señales de entrada debe ser igual al número de señales de salida.

Para obtener más información, consulte Función de transferencia.

Argumentos de salida

contraer todo

Estimación de la función de transferencia, devuelta como vector, matriz o arreglo tridimensional.

Frecuencias normalizadas, devueltas como vector columna de valor real.

Frecuencias cíclicas, devueltas como vector columna de valor real.

Más acerca de

contraer todo

Función de transferencia

La relación entre la entrada x y la salida y se modela mediante la función de transferencia lineal e invariable en el tiempo txy. En el dominio de la frecuencia, Y(f) = H(f)X(f).

  • En un sistema de una sola entrada y una sola salida, la estimación de H1 de la función de transferencia viene dada por

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

    , donde Pyx es la densidad espectral de potencia cruzada de x e y, y Pxx es la densidad espectral de potencia de x. Esta estimación supone que el ruido no está correlacionado con la entrada del sistema.

    En los sistemas multi-entrada y multi-salida (MIMO), el estimador H1 se convierte en

    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

    para las entradas m y las salidas n, donde:

    • Pyixk es la densidad espectral de potencia cruzada de la entrada k y la salida i.

    • Pxixk es la densidad espectral de potencia cruzada de las entradas k e i.

    En el caso de que haya dos entradas y dos salidas, el estimador es la matriz

    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).

  • En un sistema de una sola entrada y una sola salida, la estimación de H2 de la función de transferencia viene dada por

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

    , donde Pyy es la densidad espectral de potencia de y y Pxy = P*yx es el complejo conjugado de la densidad espectral de potencia cruzada de x e y. Esta estimación supone que el ruido no está correlacionado con la salida del sistema.

    En los sistemas MIMO, el estimador de H2 solo está bien definido para un número de entradas y salidas igual: n = m. El estimador se convierte en

    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,

    , donde:

    • Pyiyk es la densidad espectral de potencia cruzada de las salidas k e i.

    • Pxiyk es el complejo conjugado de la densidad espectral de potencia cruzada de la entrada i y la salida k.

Algoritmos

tfestimate utiliza el método del periodograma promediado de Welch. Para obtener más detalles, consulte pwelch.

Referencias

[1] Vold, Håvard, John Crowley, and G. Thomas Rocklin. “New Ways of Estimating Frequency Response Functions.” Sound and Vibration. Vol. 18, November 1984, pp. 34–38.

Capacidades ampliadas

Historial de versiones

Introducido antes de R2006a

expandir todo

Consulte también

| | |