Tema 2

Tcnicas de identificacin no paramtricas

Ejercicios y problemas resueltos
Problema 2.1.
Se quiere disear una seal de excitacin Random discreta de media 5, desviacin estndar 2
(o varianza 4), con 200 elementos y periodo de muestreo de 0.2 segundos. Evaluar utilizando
como herramienta el programa Matlab, la variacin temporal y la densidad espectral
estimada. Calcular su energa total o RMS
Solucin
El programa Matlab dispone de una funcin randn que permite generar una secuencia de nmeros
aleatorios de media cero y variancia 1. Por lo tanto aplicando:
>>u=5+2*randn(N,1); con N=200
es posible generar un vector con N filas y una columna de media 5 y desviacin estndar 2.
La funcin spectrum estima la densidad espectral de la seal utilizando el mtodo de Welch.
>> P = spectrum(u)
la visualizacin del espectro se realiza utilizando el comando specplot
>>specplot(P,Fs)
siendo Fs la frecuencia de muestreo del seal.
Observar que la frecuencia mxima es de: (N-1)/N*Fs/2, ya que se evala en la regin entre 0 y .
La figura muestra la variacin temporal y la densidad espectral de la seal diseada. Para la
representacin temporal, generamos el vector tiempo como:
>>t=0:0.2:(200-1)*0.2; %variaciones de 0.2
>>%obtencin de las curvas
>> subplot(211),plot(t,u),title('variacion temporal'), xlabel('tiempo'),ylabel('amplitud')
>> subplot(212),specplot(P,1/0.2)
Aplicando la ecuacin (2.2) se determina la energa total o RMS de la seal diseada:
>> sqrt(sum(u.^2)/N)

Problema 2.2.
Se quiere filtrar la seal diseada en el problema anterior de forma que tenga un ancho de
banda de 1.5Hz. Observar la variacin temporal y espectral de la nueva seal.
Solucin
El programa Matlab dispone de una pantalla grfica para disear filtros de distintas caractersticas,
fdatool. En nuestro caso se pretende disear un filtro discreto IIR pasa bajos de orden 2. El mtodo
utilizado ser butterword.
>> fdatool
Nos aparece la pantalla de diseo (figura
P2.2), en la que debemos seleccionar el
mtodo de diseo del filtro, el orden, la
frecuencia de muestreo de la seal (1/0.2)
y la frecuencia de corte. Una vez definidos
los parmetros podemos disear el filtro,
ver figura.

En file seleccionar export, definiendo la matriz SOS y ganancia G. Obtenindose:

>> SOS
SOS =
1.0000
>> G
G=

2.0000

1.0000

1.0000

0.3695

0.1958

0.3913

La funcin de transferencia del filtro discreto es:

> NUM=G*SOS(1:3)
NUM =
0.3913

0.7827

0.3913

>> DEN=SOS(4:6)
DEN =
1.0000

0.3695

0.1958

Seal filtrada
>> sys=tf(NUM,DEN,0.2)
Transfer function:
0.3913 z^2 + 0.7827 z + 0.3913
-----------------------------z^2 + 0.3695 z + 0.1958

Aplicacin del filtro diseado a

determinada en el problema anterior:

la

seal

>uf=filter(NUM,DEN,u);
>Pf=spectrum(uf);
>specplot(Pf,1/0.2);

Espectro de la seal filtrada

Problema 2.3. Se quiere generar una seal MLBS con las siguientes caractersticas: ancho de
banda 6 Hz, con una frecuencia mnima de 0.1 Hz, 4 muestras por cambio de registro y 2
periodos. Mostrar tanto la variacin temporal como frecuencial.
Solucin
1. Determinar t. Ya que la potencia media del ancho de banda se da aproximadamente en 0.443/t,
para que esta est situada a 6Hz el valor de t debe ser: t=0.443/60.074.
2. Ajuste de M. Al ser fmin0.1, el valor de M debe estar prximo a 135. Segn la tabla 2.1 el M mas
prximo a este valor es cuando en nmero de registros es de 7 y por lo tanto M = 127. El valor de fmin
ser de: fmin = 1/(M*t ) = 0.1064.
3. La frecuencia de muestreo de la seal ser de t=0.074, ya que se toman 4 muestras por cambio
de registro ts = t/4 = 0.0185.
Se ha programado una funcin Matlab de nombre sbina
>> [t,s,f,S]=sbina(n, fmax, sdt,amp,periods,invre)
con entradas
n
nmero de registros
fmax
frecuencia mxima de inters en el diseo de la seal
sdt
numero de muestras a adquirir por registro
amp
amplitud de la seal de entrada (-a,a)
periods
nmero de periodos de la seal de entrada
invre
seal con armnicos impares (odd harmnic)
Como salidas de la funcin tenemos:
t
s
f
S

vector tiempo
seal obtenida
frecuencia de la seal
transformada de Fouries de la seal

En nuestro caso tendremos por lo tanto:

>> [t,s,f,S]=sbina(7, 6, 4,1,2)
Los resultados obtenidos se muestran en la siguiente figura P2.6
MLBS [ fmax: 6 Hz, regs: 7, samples/bit: 4, fund: 0.107143 Hz, fr: 54.4286 Hz, amp: 1.00 ]
1
0.5
0
Amplitude
-0.5
-1
0

8
10
12
Time (s) [ 1016 Samples ]

14

16

18

Spectrum of above
0.2

0.15

0.1
Amplitude
0.05

0
0

5f(-3db)

10

15
Frequency (Hz)

20

25

30

Problema 2.4
Se pretende disear 3 periodos de dos seales multisinusoidal con las siguientes
caractersticas: 0.04 Hz de frecuencia fundamental, de 60 armnicos y una amplitud de 1
parar cada componente, una de ellas se disear utilizando una la frecuencia random y la
segunda utilizando el mtodo de Schroeder. La frecuencia de adquisicin de datos ser de 50
Hz. Se pretende comparar las dos seales obtenidas desde el campo temporal y frecuencial.
Solucin
Definimos las caractersticas de las seales
>>harm=60; fund=0.04; a=1; Fs=50; period=3;
>>Tsamples=round(Fs/fund);
>>samples=Tsamples*period;
>>t=(0:1/Fs:samples/Fs)-1/Fs;
Considerando una fase random
>>FI=2*pi*rand(harm,1);
>>u=zeros(size(t));
>>Fv=fund*(1:harm);
>>for i=1:harm
>>
u=u+cos(2*pi*Fv(i)*t+FI(i)*ones(size(t)));
>>end
>>plot(t,u)
Para hacer el estudio espectral de una seal con mas de un periodo, la frecuencia queda definida
como:
>>n=length(u);
>>f=fund/period*(1:period*harm)';
>>Udat=fft(u,n);
>>Udat=Udat(1:(n/2+1));
>>Udat(2:n/2)=2*Udat(2:n/2)/n;
>>% espectro de la senyal de entrada
>>Guu=conj(Udat).*Udat;
>>plot(f,Guu(2:period*harm+1)), title('Espectro'),xlabel('Frecuencia')
Espectro

Respuesta temporal
1.4

15

1.2

10

5
0.8

0
0.6

-5
0.4

-10

-15
0

0.2

10

20

30

40
tiempo

50

60

70

80

0.5

1.5
Frecuencia

Con un factor de cresta de 2.4774.

2.5

En el caso de utilizar el mtodo de Schroeder, la fase se determina:

>>FI1=pi*((1:harm)').^2/harm;
>>u1=zeros(size(t));
>>for i=1:harm
>>
u1=u1+cos(2*pi*Fv(i)*t+FI1(i)*ones(size(t)));
>>end
>>plot(t,u1)
>>U1dat=fft(u1,n);
>>U1dat=U1dat(1:(n/2+1));
>>U1dat(2:n/2)=2*U1dat(2:n/2)/n;
>>% espectro de la senyal de entrada
>>Guu1=conj(U1dat).*U1dat;
>>plot(f,Guu1(2:period*harm+1))
1.4

10
8

1.2

6
1

4
0.8

2
0

0.6

-2
0.4

-4
-6

0.2

-8
0

-10
0

10

20

30

40

50

60

70

0.5

1.5

2.5

80

Con un factor de cresta de 1.7386.

Problema 2.5.
Se tiene un sistema representado por la funcin de transferencia discreta con un periodo de
muestreo de 1 segundo:
q 1 + 0.5q 2
y (t ) =
u (t ) + v (t )
1 1.5q 1 + 0.7 q 2
a) Comparar la respuesta impulsional del sistema en el caso en que no haya ruido, v(t)=0, y
en el caso en que el ruido sea coloreado con una funcin de transferencia:

v (t ) =

1 q 1 + 0.2q 2
e (t )
1 1.5q 1 + 0.7 q 2

siendo e(t) un ruido blanco de varianza 0.5.

b) Determinar los valores de la respuesta impulsional en cada caso y calcular el error
cuadrado acomunado en presencia de ruido.
Solucin
a) Para determinar la respuesta impulsional del sistema utilizando el programa MATLAB v5 y v6, se
debe proceder de la siguiente forma:

>> % simulamos el sistema con la instruccin dimpulse considerando que v(t)=0

>> % la evaluacin se realizar sobre 40 muestras
>> yimp=dimpulse([0 1 0.5],[1 -1.5 0.7],40);
>> % simulamos la parte estocstica del sistema a una entrad randon
>> e=sqrt(0.5)*randn(40,1);
>> v=dlsim([1 1 0.2],[1 1.5 0.7],e);
>> yv=yimp+v;
Respuesta impulsional sin ruido (-) y con ruido (--)
>> t=1:40;
2.5
>> plot(t,yimp,t,yv,)
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2

10

15

20

25

30

35

40

tiempo

b) Los valores de la respuesta impulso (yimp) en el caso en que no hay ruido son una secuencia de
valores infinita:
g(t)= [ 0 1.0000 2.0000 2.3000 2.0500 1.4650 0.7625 0.1182 -0.3564 -0.6173 ...]
esta secuencia de valores se obtiene dividiendo el numerador por el denominador de la funcin
de transferencia.
En el caso en que si hay ruido (yv), la secuencia que se obtiene es:
g(t)= [0.3955 1.5360 1.7876 2.3558 1.0422 0.3852 -0.2976 -0.2088 -1.6083 -1.9502 ...]
con un error cuadrado acumulado de: > sum((yimp-yv).^2)/40 = 0.6006

Problema 2.6
Se desea estimar la respuesta impulsional del sistema que tienen por funcin de transferencia:

G( s ) =

1
e 2 s
4s + 4s + 1
4

a partir de la respuesta escaln del mismo. Comparar la respuesta impulsonal real con la
respuesta estimada.
Solucin
Para poder estimar la respuesta impulsional del sistema a partir del sistema real se debe disponer de
datos muestreados. Aplicando la funcin step es posible simular la respuesta de un sistema continuo
a un escaln.

>> %sistema sin retardo puro

>> sys1=tf(1,[4 4 1]);
>> %aproximamos el retardo puro por un sistema de orden 10
>> [nd,dd]=pade(2,10);
>> sys2=tf(nd,dd);
>> % realizamos el producto de las dos funciones de transferencia
>> sys=series(sys1,sys2);
>> %simulamos el sistema de forma continua y muestreamos cada segundo
>> t=0:0.1:20;
>> td=0:20;
>> y=step(sys,t);
>> yd=step(sys,td);
>> subplot(211),plot(t,y,td,yd,x),title(Respuesta escaln), xlabel(tiempo)
>> % las etapas indicadas hasta este punto podran ser realizadas directamente en
>> SIMULINK
>>
>>
>>
>>
>>

%estimamos respuesta impulso con la ecuacin (2.38)

g=yd(2:21)-yd(1:20);
yi=impulse(sys,t);
subplot(212),plot(t,yi),hold on,bar(td,g,0.1), title(Respuesta impulsional),xlabel(tiempo)
axis([0 25 0 0.2])
Respuesta escaln
1
0.8
0.6
0.4
0.2
0
0

10
12
14
tiempo
Respuesta impulsional

16

18

20

0.2
0.15
0.1
0.05
0
0

10

15

20

25

tiempo

Problema 2.7
Se desea utilizar el anlisis de correlacin para estimar la respuesta impulso del sistema
descrito en el problema 2.5. Con el fin de realizar un estudio comparativo de lo explicado en
el tema se propone comparar la respuesta impulso real del sistema con la estimada en los
siguientes casos:
a) Cuando v(t)=0 y la seal de excitacin es un ruido blanco de amplitud 1.
b) Cuando v(t) es un ruido blanco, no correlacionado con la entrada, de variancia 0.3. La
seal excitacin es la misma que en el caso anterior.
Solucin
a) y b) En primer lugar se van a generar las seales de excitacin y perturbacin con las que se va a
simular el sistema. Generaremos 500 datos con un muestreo de 1 seg.
>> t=(0:500);
>> randn(seed,5);

>>
>>
>>
>>

u=randn(size(t));
randn(seed,0);
e=sqrt(0.3)*randn(size(t));
subplot(2,2,1), plot(t,u,t,e,), title(Evolucin temporal de u(t) y e(t))

Para evaluar las seales diseadas se har un estudio de correlacin y correlacin cruzada. La
funcin de MATLAB covf permite realizar este estudio.
R=covf(Z,M);
donde: Z es una matriz de N*nz, generalmente Z=[y u] o el conjunto de seales a evaluar;
M es el valor mximo de retardo para el cual se quiere evaluar la funcin de covariancia.
R es la funcin de covariancia E zi(t)*z j(t+k), k tiene los valores de 1 a M.
>> R=covf([e u],20);
>> subplot(2,2,3),plot(R(1,:)), title(Autocorrelacin de la seal e(t))
>> subplot(2,2,4),plot(R(4,:)), title(Autocorrelacin de la seal u(t))
>> subplot(2,2,2),plot(R(2,:)), title(Correlacin cruzada E e(t)*u(t+k))
Evolucin temporal de u(t) y e(t)

Correlacin cruzada

0.15

0.1

0.05

-2

-4

-0.05
0

200

400

600

Autocorrelacin de la seal e(t)

10

15

20

Autocorrelacin de la seal u(t)

0.4

0.3
0.5

0.2
0.1

0
-0.1

10

15

-0.5

20

10

15

20

Una vez generadas la seales y evaluadas, procedemos a la simulacin del sistema.

>> y1=dlsim([0 1 0.5],[1 1.5 0.7],u);
>> y2=y1+e;
>> ti=0:50;
>> yimp=dimpulse([0 1 0.5],[1 1.5 0.7],ti);
Para estimacin de la respuesta impulsional utilizaremos la funcin cra de MATLAB, adems de
calcular la funcin impulso indica con una confianza del 99% el intervalo de valores a partir de los
cuales dos seales pueden ser consideradas no correlacionadas.
[ir,R,cl] = cra(z,M,na,plot);
en donde:
z = [y u], matriz de datos
M = 50, nmero de retardos a evaluar la funcin de covariancia
na = 0, orden del filtre AR a aplicar para blanquear las entradas

plot = 2 o 1, plot = 0 visualizacin de los resultados obtenidos; plot = 1, dibuja la

funcin de transferencia impulso; y plot = 2, dibuja las covariancias y
funcin impulso.
ir = respuesta impulso estimada, ir(1)=g(0)
R=
R(:,1) retardos;
R(:,2) funcin de covariancia de la y(t);
R(:,3) funcin de covariancia de la u(t);
R(:,4) funcin de correlacin cruzada entre u(t) i y(t);
cl = nivel de confianza
>>

[IR,R,cl]=cra([y1 u], 50,0,2);

Covf for y

Covf for u

20

15
0.5

10
5

0
-5
-50

-0.5
-50

50

Correlation from u to y
3

0.4

0.2

50

Impulse response estimate

0.6

-0.2
-50

-1
-50

50

50

Para comparar la respuesta impulsional estimada con la real, podemos aplicar el siguiente comando:
>> plot(ti,IR,ti,yimp,),title(Respuesta impulso comparacin, real (--) y estimada (-))
En la grfica obtenida se observa que la respuesta impulsional estimada se aproxima a la real en los
primeros datos M<7, para valores mayores la diferencia es considerable. Otro aspecto a destacar es
que con la respuesta impulsional se observa una respuesta tpica de un sistema de 2 orden.

10

Respuesta impulso comparacin, real(--) y estimada (-)

2.5

1.5

0.5

-0.5

-1

10

15

20

25

30

35

40

45

50

En el caso b) tendremos:
>> [IRn,Rn,cln]=cra([y2 u], 50,0,2);

Covf for y

Covf for u

20

15
0.5

10
5

0
-5
-50

-0.5
-50

50

Correlation from u to y
3

0.4

0.2

50

Impulse response estimate

0.6

-0.2
-50

-1
-50

50

50

La comparacin entre la respuesta impulso real y la estimada ser:

>> plot(ti,IRn,ti,yimp,),title(Respuesta impulso comparacin, real (--) y estimada (-))

11

Respuesta impulso comparacin, real (--) y estimada (-)

2.5

1.5

0.5

-0.5

-1

10

15

20

25

30

35

40

45

50

La respuesta obtenida con ruido es muy similar a la obtenida sin ruido.

Problema 2.8
Utilizando las tcnicas frecuenciales del anlisis de Fourier y el anlisis espectral determinar
la respuesta frecuencial del sistema descrito en el problema 2.5. Como seal de excitacin
utilizar una seal multisinusoidal de caractersticas: 100 armnicos, 0.002 de frecuencia
fundamental y 3 periodos. Como seal perturbacin utilizar un ruido blanco de variancia 5.
Solucin
Archivo solP2_8. Obsrvese que la respuesta frecuencial se estima correctamente, como resultado
de aplicar la funcin emprica ETFE (figura 5) como el anlisis espectral (figura 6). La figura 7
muestra la respuesta frecuencial terica, resultado de utilizar el comando dbode.

Nota: al utilizar en este caso una seal peridica se ha utilizado directamente la

transformada de Fourier para la estimacin de la funcin de transferencia emprica
(respuesta frecuencial) tanto a partir del anlisis frecuencial como espectral. En el caso en
que las seales utilizadas no sean peridicas, deben utilizarse directamente los comandos
spa y/o etfe.

12

Ejercicios y problemas propuestos

E2.1 Se quiere identificar un modelo matemtico de un horno tubular que calienta el aire a
partir de una resistencia elctrica. La entrada al proceso es la tensin suministrada a la
resistencia calefactora y la salida es la seal enviada por el sensor de temperatura situado a la
mitad de este horno. Como seal de excitacin se ha utilizado una seal binaria. Las seales
de entrada y salida se muestran en las siguientes figuras:
8

entrada

2
240

250

260

250

260

270
tiempo

280

290

300

280

290

300

salida
0.5
0.4
0.3
0.2
0.1
0
240

270
tiempo

Se pide que a la vista de los resultados se comenten los tratamientos previos que deberan
realizarse antes de proceder a la identificacin del modelo.
E2.2. Se desea identificar el modelo discreto de un proceso. Partiendo de un ensayo
preliminar se ha determinado que el intervalo de frecuencias que describen el proceso es:
[fmin , fmax ] = [0.5, 4 ] Hz. Comentar:
- A partir de que ensayos preliminares es posible determinar este intervalo?
- Cul es el periodo de adquisicin de datos que utilizaras para identificar este proceso?
- Si es menester filtrar los datos con un filtro antialias, que frecuencia de corte escogeras
para el diseo de este filtro?
E2.3- Ajustar por el mtodo de Kpfmller y Strejc un modelo que se ajuste a la respuesta
escaln del problema 2.6. Comparar los resultados del modelo determinado con el modelo
real.
E2.4.- Estimar la respuesta impulso del problema 2.7 pero en el caso en que se disponga de
una entrada coloreada. Utilizar para ello el comando cra seleccionando el orden adecuado del
filtro, na, para blanquear la seal de entrada..

13

E2.5- Variar las caractersticas de la seal de entrada (frecuencia fundamental y nmero de

armnicos) y el nivel de ruido del problema 2.8 y observar como se modifica la respuesta
frecuencial del sistema.
E2.6- Se desea estimar la respuesta frecuencial de un sistema, para ello se ha excitado dicho
sistema con una seal de entrada no peridico u(t) obtenindose como resultado del mismo la
seal y(t), ambas se encuentran el fichero adjunto E26.mat. Utilizar las funciones de Matlab
etfa y spa para estimar la funcin de transferencia emprica a partir del anlisis de Fourier y
del anlisis espectral, respectivamente. Comparar la respuesta frecuencial utilizando ambos
mtodos y variando la longitud de la ventana (M o 1/), justificar los resultados obtenidos.
Podrais indicar de que tipo de sistema se trata?, cual es su ancho de banda?.
E2.7- Partiendo de los mismos datos del ejercicio anterior, E26.mat, estimar la respuesta
impulso del sistema. A partir de los resultados de que tipo de sistema se trata?. Justificar la
respuesta.
E2.8- A partir de los datos de entrada y salida de un sistema, se ha estimado la funcin de
transferencia impulso obtenindose los resultados mostrados en la figura.
Respuesta impulso estimada
1

0.8

0.6

0.4

0.2

-0.2
0

10
retards

15

20

Indicar la informacin podemos obtener a partir de ella. Justificar la respuesta.

E2.9- La funcin de transferencia emprica obtenida a partir del estudio espectral de un
sistema excitado mediante un ruido blanco, se observa en la figura siguiente:
AMPLITUDE PLOT

1 01

10

10

-1

1 0- 2
-2
10

-1

1 0 frequency (rad/sec) 1 0

10

PHASE PLOT

-200
phase

-400
10 -2

1 0- 1

1 00
frequency (rad/sec)

14

1 01

Que frecuencia de muestreo escogerais para dicho sistema? Justificar la respuesta. Indicar
tambin de forma aproximada el tipo de sistema de que se trata.

E2.10-. La funcin de auto covarianza de la entrada y la covarianza cruzada entre la entrada y

la salida estn representadas en las siguientes figuras:
Ifunci dauto-covariana

funci de covariana creuada

0.2

0.8

0.15

0.6

0.1

0.4

0.05

0.2

-0.05

-0.2
0

10
lags

15

20

-0.1
0

Sabrais deducir la funcin de transferencia impulso?.

15

10
lags

15

20