Prácticas de DSP. Procesadores digitales de señales

Procesadores digitales de señales

1. Práctica 1. Sistemas LTI y análisis de Fourier.
2. Práctica 2. Muestreo.
3. Práctica 3 La Transformada Z.
4. Práctica 4. La DFT y la FFT. y Práctca 5.Diseño de filtros digitales.

Nota: Para la realización de esta memoria se ha utilizado el procesador de textos "Microsoft Word 2002 XP"

PRACTICA 1. Sistemas LTI y análisis de Fourier.

>> help conv

La función conv(x,y) sirve para calcular la convolucion de x e y. El vector que resulta es de longitud Length(a)+length(b)-1. Si A y B son vectores de coeficientes polinómicos, la convolución es equivalente a multiplicar los dos polinómios.

CONV Convolution and polynomial multiplication.

C = CONV(A, B) convolves vectors A and B. The resulting

vector is length LENGTH(A)+LENGTH(B)-1.

If A and B are vectors of polynomial coefficients, convolving

them is equivalent to multiplying the two polynomials.

See also DECONV, CONV2, CONVN, FILTER and, in the Signal

Processing Toolbox, XCORR, CONVMTX.

>> help filter

La función filter() permite calcular la salida de un sistema LTI causal cuando éste es definido mediante una ecuación en diferencias de coeficientes constantes.

Si x es un vector con la secuencia de entrada y a y b son vectores con los coeficientes de la ecuación , la expression y=filter(b,a,x) devuelve un vector y con la respuesta del sistema LTI causal.

FILTER One-dimensional digital filter.

Y = FILTER(B,A,X) filters the data in vector X with the

filter described by vectors A and B to create the filtered

data Y. The filter is a "Direct Form II Transposed"

implementation of the standard difference equation:

a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)

- a(2)*y(n-1) - ... - a(na+1)*y(n-na)

If a(1) is not equal to 1, FILTER normalizes the filter

coefficients by a(1).

FILTER always operates along the first non-singleton dimension,

namely dimension 1 for column vectors and non-trivial matrices,

and dimension 2 for row vectors.

[Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final

conditions, Zi and Zf, of the delays. Zi is a vector of length

MAX(LENGTH(A),LENGTH(B))-1 or an array of such vectors, one for

each column of X.

FILTER(B,A,X,[],DIM) or FILTER(B,A,X,Zi,DIM) operates along the

dimension DIM.

See also FILTER2 and, in the Signal Processing Toolbox, FILTFILT.

Ejercicio 1.2 La función "conv"

c=1; % Control

if c==1

x=ones(1,6); % Secuencia

y=conv(x,x); % Convolución x*x

Nx=0; % Indices

Mx=5;

Ny=Nx+Nx;

My=Mx+Mx;

ny=Ny:My;

figure % Nueva figura

stem(ny,y); % Resultados

title('Figura 1.1')

end

Ejercicio 1.2.1 Promediador Móvil

c=1; % Control

if c==1

n=-30:30; % Tiempo discreto

h=1/7*ones(1,7); % Respuesta impulsional

x=sin(2*pi*n/13); % Secuencia

y=conv(h,x); % Convolución

V=zeros(1,81); % Vector de ceros

Nh=-3; % Indices h[n]

Mh=3;

Nx=-33; % Indices x[n]

Mx=33;

Ny=Nx+Nx; % Indices y[n]

My=Mx+Mx;

ny=Ny:My;

V(8:8+66)=y; % Salida ajustada

nV=-40:40; % Tiempo discreto

figure % Nueva figura

stem(nV,V); % Resultados

title('Figura 1.2')

end

Ejercicio 1.3 La Función Filter

c=1; % Control

if c==1

a=[1 2]; % Coeficientes y[n]

b=[1 -3]; % Coeficientes x[n]

n=0:5; % Tiempo discreto

x=n; % Secuencia x[n]

y=filter(b,a,x); % Filtrado de x[n]

figure % Nueva figura

stem(x,y); % Resultados

title('Figura 1.3')

end

Para ver el gráfico seleccione la opción "Descargar" del menú superior

Ejercicio 1.3.1

c=1; % Control

if c==1

%--- Punto 1

% Coeficientes y[n]

a=[1 (-1.8*cos(pi/16)) 0.81];

b=[1 1/2]; % Coeficientes x[n]

n=0:50; % Tiempo discreto

x=zeros(1,51); % Vector de ceros

x(1)=1;

y=filter(b,a,x); % Filtrado de x[n]

figure % Nueva figura

stem(n,y); % Resultados

title('Figura 1.4')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

%---

n=-10:50; % Tiempo discreto

x=zeros(1,61); % Vector de ceros

x(11)=1;

y=filter(b,a,x); % Filtrado de x[n]

figure % Nueva figura

stem(n,y); % Resultados

title('Figura 1.4 Modificada la x[n]')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

%--- Punto 2

n=-10:50; % Tiempo discreto

x=zeros(1,61); % Vector de ceros

x(11:61)=1;

y=filter(b,a,x); % Filtrado de x[n]

figure % Nueva figura

stem(n,y); % Resultados

title('Figura 1.5')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

%---

n=-10:50; % Tiempo discreto

x=zeros(1,61); % Vector de ceros

x(11:61)=1;

y=filter(b,a,x); % Filtrado de x[n]

yc=cumsum(y,1); % "Filtrado" de y[n]

figure % Nueva figura

stem(n,yc); % Resultados

title('Figura 1.5 con Cumsum')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

%---

n=-10:50; % Tiempo discreto

x=zeros(1,61); % Vector de ceros

x(11:61)=1;

y=filter(b,a,x); % Filtrado de x[n]

y1=diff(y); % Recuperación

n1=-9:50; % Nuevo intervalo

figure % Nueva figura

stem(n1,y1); % Resultados

title('Respuesta al impulso con diff')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

%--- Punto 3

a=[1]; % Coeficientes y[n]

b=1/7*ones(1,7); % Coeficientes x[n]

n=-30:30; % Tiempo discreto

x=sin(2*pi*n/13); % Secuencia

y=filter(b,a,x); % Filtrado de x[n]

figure % Nueva figura

stem(n,y); % Resultados

title('Figura 1.2 con función filter')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

end

Ejercicio 1.4 Cancelación de ecos

c=1; % Control

if c==1

% Lectura WAV

[voz,fs]=wavread('ho1.wav');

sound(voz,fs) % Reproducción

% [voz fs]=loadw16('ho1.wav');

pause % Control de reproducción

%-- Reproducción con ecos 1.4.1

N=1000; % Retardo

alfa=0.5; % Factor de atenuación

% Sistema

b=[1 zeros(1,999) alfa];

% Filtrado

ve=filter(b,1,voz);

sound(ve,fs) % Reproducción con eco

%---

figure % Nueva figura

stem(b); % Resultados

title('Figura 1.6')

pause % Control de reproducción

%-- Sistema inverso 1.4.2

% Coeficientes y[n]

a=[1 zeros(1,999) alfa];

b=[1]; % Coeficientes x[n]

y=filter(b,a,ve); % Señal sin ecos

sound(y,fs) % Reproducción sin eco

%--- Cálculo del sistema inverso

n=0:4000; % Intervalo discreto

a=[1 zeros(1,999) alfa]; % Coeficientes de y[n]

b=1; % Coeficientes de x[n]

x=zeros(1,4000); % Secuencia impulso

x(1)=1;

hi=filter(b,a,x); % Filtrado

figure % Nueva figura

subplot(211)

plot(n,hi) % Respuesta al impulso

title('Respuesta al impulso nuevo sistema')

b=[1 zeros(1,999) alfa]; % Sistema con eco

he=filter(b,1,x); % Filtrado

ht=conv(hi,he); % Convolución

subplot(212)

plot(n,ht) % Respuesta al impulso

title('Convolución sistema inverso y con eco')

end

Ejercicio 1.5.1 Señales de duración finita

c=1; % Control

if c==1

p=ones(1,15); % Pulso de duración 15

ep=fft(p,256); % Espectro del pulso

figure % Nueva figura

subplot(221)

plot(abs(ep)) % Módulo del espectro

title('Espectro del pulso ( 0 y 1 )')

ep=fftshift(ep); % Espectro del pulso

% Frecuencia

f=linspace(-0.5,0.5-(1/256),256);

subplot(222)

plot(f,abs(ep)) % Módulo del espectro

title('Espectro del pulso ( -0.5 y 0.5 )')

subplot(223)

plot(f,ep) % Módulo y fase del espectro

title('Módulo y Fase del espectro')

subplot(224)

plot(f,imag(ep)) % Módulo y fase del espectro

title('Fase del espectro')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

ep=ep.*exp(j*2*pi*f*7);

figure % Nueva figura

subplot(211)

plot(f,real(ep)) % Módulo del espectro

title('Módulo del espectro')

subplot(212)

plot(f,imag(ep)) % Módulo del espectro

title('Fase del espectro')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

Ejercicio 1.5.2 Ecuaciones en diferencia

c=1; % Control

if c==1

% Espectro señal exponencial

[H,W]=freqz(1,[1 -0.9],256);

figure % Nueva figura

subplot(221)

plot(W/(2*pi),abs(H)) % Módulo del espectro

title('Módulo del espectro ( 0 y pi )')

% Espectro señal exponencial

[H,W]=freqz(1,[1 -0.9],256,'whole');

subplot(222)

plot(W/(2*pi),abs(H)) % Módulo del espectro

title('Módulo del espectro ( 0 y 2*pi )')

W=W-pi; % Espectro entre -0.5 y 0.5

H=fftshift(H);

subplot(223)

plot(W/(2*pi),abs(H)) % Módulo del espectro

title('Módulo del espectro ( -0.5 y 0.5 )')

subplot(224)

plot(W/(2*pi),angle(H)) % Fase del espectro

title('Fase del espectro ( -0.5 y 0.5 )')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

end

PRACTICA 2. Muestreo.

>> help fftshift

FFTSHIFT es útil para visualizar la Transformada de Fourier con la componente de la frecuencia cero en el centro del espectro. Por vectores, FFTSHIFT(X) intercambia las mitades izquierdas y derechas de X. Por matrices, FFTSHIFT(X) intercambia los primeros y terceros cuadrantes y los segundos y cuartos cuadrantes.

FFTSHIFT Shift zero-frequency component to center of spectrum.

For vectors, FFTSHIFT(X) swaps the left and right halves of

X. For matrices, FFTSHIFT(X) swaps the first and third

quadrants and the second and fourth quadrants. For N-D

arrays, FFTSHIFT(X) swaps "half-spaces" of X along each

dimension.

FFTSHIFT(X,DIM) applies the FFTSHIFT operation along the

dimension DIM.

FFTSHIFT is useful for visualizing the Fourier transform with

the zero-frequency component in the middle of the spectrum.

See also IFFTSHIFT, FFT, FFT2, FFTN.

Ejercicio 2.2 SEÑALES SINUSOIDALES

c=1; % Control

if c==1

fs=8000; % Frecuencia de muestreo

T=1/fs; % Período de muestreo

f0=1000; % Frecuencia señal sinusoidal

n=0:7999; % Tiempo discreto

x=sin(2*pi*f0*n*T); % Secuencia x[n]

sound(x,fs) % Reproducción de x[n]

figure % Nueva figura

t=n*T; % Intervalo temporal

% Visualiza resultados

subplot(211), stem(n(1:20),x(1:20))

title('Señal sinusoidal x[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

subplot(212), plot(t(1:20),x(1:20))

title('Señal sinusoidal x(t)'), grid

xlabel('Tiempo (seg)'), ylabel('Unidades Lineales')

 figure % Nueva figura

N=length(x); % Longitud de x[n]

% Espectro de x[n]

X=fftshift(fft(x,N))/N;

% Frecuencia discreta

f=linspace(-0.5,0.5-(1/N),N);

% Visualiza resultados

subplot(211), plot(f,abs(X))

title('Espectro discreto de x[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

subplot(212), plot(f*fs,abs(X))

title('Espectro continuo de x[n]'), grid

xlabel('Frecuencia (Hz)'), ylabel('Unidades Lineales')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

end

Ejercicio 2.3 SEÑAL CHIRP

c=1; % Control

if c==1

fs=8000; % Frecuencia de muestreo

T=1/fs; % Período de muestreo

f0=3000; % Frecuencia señal sinusoidal

beta=2000; % Beta

n=0:1*(7999); % Tiempo discreto

% Secuencia x[n]

x=sin(2*pi*f0*n*T+0.5*beta*(n*T).^2);

sound(x,fs) % Reproducción de x[n]

figure % Nueva figura

t=n*T; % Intervalo temporal

% Visualiza resultados

subplot(211), stem(n(1:100),x(1:100))

title('Señal chirp x[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

subplot(212), plot(t(1:100),x(1:100))

title('Señal chirp x(t)'), grid

xlabel('Tiempo (seg)'), ylabel('Unidades Lineales')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

figure % Nueva figura

N=length(x); % Longitud de x[n]

% Espectro de x[n]

X=fftshift(fft(x,N))/N;

% Frecuencia discreta

f=linspace(-0.5,0.5-(1/N),N);

% Visualiza resultados

subplot(211), plot(f,abs(X))

title('Espectro discreto de x[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

subplot(212), plot(f*fs,abs(X))

title('Espectro continuo de x[n]'), grid

xlabel('Frecuencia (Hz)'), ylabel('Unidades Lineales')

Para ver el gráfico seleccione la opción "Descargar" del menú superior

end

Ejercicio 2.4 CAMBIO DE LA FRECUENCIA DE MUESTREO

Ejercicio 2.4.1 Diezmado

c=1; % Control

if c==1

n=0:124; % Tiempo discreto

x1=0.4*(n-62); % Secuencia x1[n]

x1=(sinc(x1)).^2;

x2=0.2*(n-62); % Secuencia x2[n]

x2=(sinc(x2)).^2;

figure % Nueva figura

% Visualiza resultados

subplot(211), plot(n,x1)

title('Secuencia x1[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

subplot(212), plot(n,x2)

title('Secuencia x2[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

figure % Nueva figura

% Espectro de x1[n]

X1=fftshift(fft(x1,2048));

% Espectro de x2[n]

X2=fftshift(fft(x2,2048));

% Frecuencia discreta

f=linspace(-0.5,0.5-(1/2048),2048);

% Visualiza resultados

subplot(211), plot(f,abs(X1))

title('Espectro discreto de x1[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

subplot(212), plot(f,abs(X2))

title('Espectro discreto de x2[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

%--- Diezmado por un factor M

% Secuencia diezmada xd1[n]

xd1=x1(1:2:length(x1));

% Secuencia diezmada xd2[n]

xd2=x2(1:2:length(x2));

figure % Nueva figura

% Visualiza resultados

subplot(211), plot(xd1)

title('Secuencia diezmada xd1[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

subplot(212), plot(xd2)

title('Secuencia diezmada xd2[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

figure % Nueva figura

% Espectro de x1[n]

Xd1=fftshift(fft(xd1,2048));

% Espectro de x2[n]

Xd2=fftshift(fft(xd2,2048));

% Frecuencia discreta

f=linspace(-0.5,0.5-(1/2048),2048);

% Visualiza resultados

subplot(211), plot(f,abs(Xd1))

title('Espectro discreto de xd1[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

subplot(212), plot(f,abs(Xd2))

title('Espectro discreto de xd2[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

%--- Filtro antialiasing

fc=0.25; % Frecuencia de corte 0.5/M

% Respuesta al impulso del filtro

h=2*fc*sinc(2*fc*(-32:32)).*(hamming(65)');

figure % Nueva figura

freqz(h,1) % Visualiza resultados

%--- Filtrando x1[n]

xf1=conv(h,x1); % Secuencia filtrada xf2[n]

% Secuencia filtrada diezmada xfd2[n]

xfd1=xf1(1:2:length(xf1));

figure % Nueva figura

% Espectro de xfd2[n]

Xfd1=fftshift(fft(xfd1,2048));

% Frecuencia discreta

f=linspace(-0.5,0.5-(1/2048),2048);

% Visualiza resultados

plot(f,abs(Xfd1))

title('Espectro discreto de xfd1[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

end

Ejercicio 2.4.2 Interpolación

c=1; % Control

if c==1

n=0:124; % Tiempo discreto

x1=0.4*(n-62); % Secuencia x1[n]

x1=(sinc(x1)).^2;

x2=0.2*(n-62); % Secuencia x2[n]

x2=(sinc(x2)).^2;

%--- Interpolando por un factor L

figure % Nueva figura

% Visualiza resultados

stem(x1(52:72))

title('Secuencia x1[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

figure % Nueva figura

stem(x2(52:72))

title('Secuencia x2[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

% Secuencia interpolada xe1[n]

xe1=zeros(1,3*length(x1));

xe1(1:3:length(xe1))=x1;

% Secuencia interpolada xe2[n]

xe2=zeros(1,3*length(x2));

xe2(1:3:length(xe2))=x2;

figure % Nueva figura

% Visualiza resultados

stem(xe1(177:197))

title('Secuencia interpolada xe1[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

figure % Nueva figura

stem(xe2(177:197))

title('Secuencia interpolada xe2[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

figure % Nueva figura

% Espectro de xe1[n]

Xe1=fftshift(fft(xe1,2048));

% Espectro de xe2[n]

Xe2=fftshift(fft(xe2,2048));

% Frecuencia discreta

f=linspace(-0.5,0.5-(1/2048),2048);

% Visualiza resultados

subplot(211), plot(f,abs(Xe1))

title('Espectro discreto de xe1[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

subplot(212), plot(f,abs(Xe2))

title('Espectro discreto de xe2[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

%--- Filtro interpolador

fc=0.16; % Frecuencia de corte 0.5/L

% Respuesta al impulso del filtro

h=(2*fc*sinc(2*fc*(-32:32)).*(hamming(65)'))*3;

figure % Nueva figura

freqz(h,1) % Visualiza resultados

%--- Filtrando xe1[n] y xe2[n]

% Xfe1=conv(h,xe1); % Secuencia filtrada xfe1[n]

% xfe2=conv(h,xe2); % Secuencia filtrada xfe2[n]

xfe1=filter(h,1,xe1); % Secuencia filtrada xfe1[n]

xfe2=filter(h,1,xe2); % Secuencia filtrada xfe2[n]

figure % Nueva figura

% Visualiza resultados

stem(xfe1(209:229))

title('Secuencia interpolada filtrada xfe1[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

figure % Nueva figura

stem(xfe2(209:229))

title('Secuencia interpolada filtrada xfe2[n]'), grid

xlabel('n'), ylabel('Unidades Lineales')

figure % Nueva figura

% Espectro de xfe1[n]

Xfe1=fftshift(fft(xfe1,2048));

% Espectro de xfe2[n]

Xfe2=fftshift(fft(xfe2,2048));

% Frecuencia discreta

f=linspace(-0.5,0.5-(1/2048),2048);

% Visualiza resultados

subplot(211), plot(f,abs(Xfe1))

title('Espectro discreto de xfe1[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

subplot(212), plot(f,abs(Xfe2))

title('Espectro discreto de xfe2[n]'), grid

xlabel('fd=f/fs'), ylabel('Unidades Lineales')

end

Ejercicio 2.4.3 Cambio fs por un factor racional

c=1; % Control

if c==0

% Carga archivo de voz

[voz,fs]=wavread('ho1.wav');

% Versiones anteriores a MATLAB 5.3

% [voz fs]=loadw16('ho1.wav');

clc

disp('Pausa de voz... cualquier tecla para continuar')

sound(voz,fs) % Reproducción archivo de voz

pause

%--- Diezmado por un factor M

% Voz diezmada

vozd=voz(1:3:length(voz));

disp('Pausa de voz diezmada... cualquier tecla para continuar')

sound(vozd,fs/3) % Reproducción voz diezmada

pause

%--- Filtro antialiasing

fc=0.16; % Frecuencia de corte 0.5/M

% Respuesta al impulso del filtro

h=2*fc*sinc(2*fc*(-32:32)).*(hamming(65)');

%--- Filtrando voz

% Voz filtrada

vozf=filter(h,1,voz);

% Voz diezmada

vozfd=vozf(1:3:length(vozf));

disp('Pausa de voz filtrada diezmada... cualquier tecla para continuar')

sound(vozfd,fs/3) % Reproducción voz filtrada diezmada

pause

%--- Interpolando por un factor L

clc

disp('Pausa de voz... cualquier tecla para continuar')

sound(voz,fs) % Reproducción archivo de voz

pause

% Voz interpolada

voze=zeros(1,3*length(voz));

voze(1:3:length(voze))=voz;

disp('Pausa de voz interpolada... cualquier tecla para continuar')

sound(voze,3*fs) % Reproducción voz interpolada

pause

%--- Filtro interpolador

fc=0.16; % Frecuencia de corte 0.5/L

% Respuesta al impulso del filtro

h=(2*fc*sinc(2*fc*(-32:32)).*(hamming(65)'))*3;

%--- Filtrando voz

% Voz filtrada

vozfe=filter(h,1,voze);

disp('Pausa de voz filtrada expandida... cualquier tecla para continuar')

sound(vozfe,3*fs) % Reproducción voz filtrada expandida

pause

%--- Cambio de la frecuencia de muestreo fs'=fs*(L/M)

%--- Solución 1: L=10, M=8

%--- Solución 2: L=5, M=4

clc % Borra la pantalla

% Control del cambio de fs

L=5; % Factor de interpolación

M=4; % Factor de diezmado

%--- Interpolando por un factor L

disp('Pausa de voz... cualquier tecla para continuar')

sound(voz,fs) % Reproducción archivo de voz

pause

% Voz interpolada

voze=zeros(1,L*length(voz));

voze(1:L:length(voze))=voz;

%--- Filtro interpolador

fc=0.5/L; % Frecuencia de corte 0.5/L

% Respuesta al impulso del filtro

h=(2*fc*sinc(2*fc*(-32:32)).*(hamming(65)'))*L;

%--- Filtrando voz

% Voz filtrada

vozfe=filter(h,1,voze);

%--- Diezmado por un factor M

% Voz diezmada

vozfed=vozfe(1:M:length(vozfe));

% Reproducción voz filtrada, expandida y diezmada

disp('Pausa de voz filtrada, expandida y diezmada... cualquier tecla para continuar')

sound(vozfed,fs*(L/M))

pause

end

PRACTICA 3 La Transformada Z.

>> help polar

La función polar() permite hacer diagramas en forma polar, es decir, dado un numero complejo, polar() representa su módulo y fase.

POLAR(THETA, RHO) hace un diagrama usando los coordenadas polares del angulo de la THETA , en radianes, frente el radio RHO.

POLAR Polar coordinate plot.

POLAR(THETA, RHO) makes a plot using polar coordinates of

the angle THETA, in radians, versus the radius RHO.

POLAR(THETA,RHO,S) uses the linestyle specified in string S.

See PLOT for a description of legal linestyles.

See also PLOT, LOGLOG, SEMILOGX, SEMILOGY.

>> help roots

La función roots(C) encuentra las raíces polinómicas del vector C, donde los ceros y los polos son las raíces de los polinomios del numerador y del denominador de la función de transferencia, respectivemente.

ROOTS Find polynomial roots.

ROOTS(C) computes the roots of the polynomial whose coefficients

are the elements of the vector C. If C has N+1 components,

the polynomial is C(1)*X^N + ... + C(N)*X + C(N+1).

See also POLY, RESIDUE, FZERO.

>> help poly

La función poly() construye un polinomio a partir de sus raíces.

POLY Convert roots to polynomial.

POLY(A), when A is an N by N matrix, is a row vector with

N+1 elements which are the coefficients of the

characteristic polynomial, DET(lambda*EYE(SIZE(A)) - A) .

POLY(V), when V is a vector, is a vector whose elements are

the coefficients of the polynomial whose roots are the

elements of V . For vectors, ROOTS and POLY are inverse

functions of each other, up to ordering, scaling, and

roundoff error.

ROOTS(POLY(1:20)) generates Wilkinson's famous example.

See also ROOTS, CONV, RESIDUE, POLYVAL.

Ejercicio 3.2 MANEJO DE FUNCIONES DE TRANSFERENCIA

Ejercicio 3.2.1 REPRESENTACION GRAFICA DE POLOS Y CEROS

c=1; % Control

if c==1

%--- Punto 1

z=0.7*exp(j*pi/3); % Variable compleja z

figure % Nueva figura

% Visualiza resultados

polar(angle(z),abs(z),'o')

title('Variable compleja en forma polar'), grid

z=0.5+j*0.5; % Variable compleja z

Para ver el gráfico seleccione la opción "Descargar" del menú superior

figure % Nueva figura

% Visualiza resultados

polar(angle(z),abs(z),'o')

title('Variable compleja en forma cartesiana'), grid

 Para ver el gráfico seleccione la opción "Descargar" del menú superior

figure % Nueva figura

% Variables complejas

z=[0.5+j*0.5 0.7*exp(j*pi/3)];

% Visualiza resultados

polar(angle(z),abs(z),'x')

title('Variables complejas'), grid

Para ver el gráfico seleccione la opción "Descargar" del menú superior

%--- Punto 2

clc % Borra la pantalla

disp('Raices del polinomio z^2-z-2')

A=[1 -1 -2]; % Polinomio z^2-z-2

r=roots(A) % Raices del polinomio z^2-z-2

%--- Punto 3

disp('Polinomio asociado con las raices...')

A=poly([-1 2]) % Polinomio asociado con las raices...

%--- Punto 4

A=[1 -1 -2]; % Polinomio z^2-z-2

figure % Nueva figura

subplot(211)

zplane(A); % Diagrama de ceros y polos

title('Diagrama de ceros y polos para A,[]'), grid

subplot(212)

zplane(A'); % Diagrama de ceros y polos

title('Diagrama de ceros y polos para A`,[]'), grid

figure % Nueva figura

zplane(A); % Diagrama de ceros y polos

title('Diagrama de ceros y polos para A,[0]'), grid

%---

% Sistema con H1(z)

bH1=roots([1]); % Ceros de H1(z)

% Polos de H1(z)

aH1=roots([1 -0.85*sqrt(2) 0.7225]);

% Sistema con H2(z)

% Ceros de H2(z)

bH2=roots([1 -0.85*sqrt(2) 0.7225]);

aH2=roots([1]); % Polos de H2(z)

% Sistema con H3(z)

bH3=roots([1]) % Ceros de H3(z)

% Polos de H3(z)

aH3=roots([1 -0.95*sqrt(2) 0.9025])

% Sistema con H4(z)

% Ceros de H4(z)

bH4=roots([1 -0.95*sqrt(2) 0.9025]);

aH4=roots([1]) % Polos de H4(z)

figure % Nueva figura

subplot(223)

zplane(bH1,aH1); % Diagrama de ceros y polos

title('Sistema con H1(z)'), grid

subplot(221)

zplane(bH2,aH2); % Diagrama de ceros y polos

title('Sistema con H2(z)'), grid

subplot(224)

zplane(bH3,aH3); % Diagrama de ceros y polos

title('Sistema con H3(z)'), grid

subplot(222)

zplane(bH4,aH4); % Diagrama de ceros y polos

title('Sistema con H4(z)'), grid

end

Ejercicio 3.2.2 RESPUESTA EN FRECUENCIA A PARTIR DE LA TZ

c=1; % Control

if c==1

%--- Punto 1

clc % Borra la pantalla

disp('Evaluación del polinomio z^2-z-2 en z=-1')

A=[1 -1 -2]; % Polinomio z^2-z-2

z=-1; % Evaluación en z=-1

polyval(A,z) % Evaluación delpolinomio z^2-z-2 en z=-1

%---

disp('Evaluación [H(0), H(0.5)]...')

% Sistema con H1(z)

cH1=[1]; % Ceros de H1(z)

% Polos de H1(z)

pH1=[1 -0.85*sqrt(2) 0.7225];

z=1; % H(0)

H1=(polyval(cH1,z))/(polyval(pH1,z));

z=-1; % H(0.5)

H1=[H1 (polyval(cH1,z))/(polyval(pH1,z))];

% Sistema con H2(z)

% Ceros de H2(z)

cH2=[1 -0.85*sqrt(2) 0.7225];

pH2=[1]; % Polos de H2(z)

z=1; % H(0)

H2=(polyval(cH2,z))/(polyval(pH2,z));

z=-1; % H(0.5)

H2=[H2 (polyval(cH2,z))/(polyval(pH2,z))];

% Sistema con H3(z)

cH3=[1]; % Ceros de H3(z)

% Polos de H3(z)

pH3=[1 -0.95*sqrt(2) 0.9025];

z=1; % H(0)

H3=(polyval(cH3,z))/(polyval(pH3,z));

z=-1; % H(0.5)

H3=[H3 (polyval(cH3,z))/(polyval(pH3,z))];

% Sistema con H4(z)

% Ceros de H4(z)

cH4=[1 -0.95*sqrt(2) 0.9025];

pH4=[1]; % Polos de H4(z)

z=1; % H(0)

H4=(polyval(cH4,z))/(polyval(pH4,z));

z=-1; % H(0.5)

H4=[H4 (polyval(cH4,z))/(polyval(pH4,z))];

H1, H2, H3, H4 % Evaluación de H1, H2, H3 y H4

%--- Punto 2

% Respuesta en frecuencia

[H,W]=freqz([1 -1 -2],1,512,'whole');

figure % Nueva figura

subplot(211)

% Módulo respuesta en frecuencia

plot(W/(2*pi),abs(H))

title('Módulo de la respuesta en frecuencia'), grid

ylabel('Unidades Lineales')

subplot(212)

% Fase respuesta en frecuencia

plot(W/(2*pi),angle(H))

title('Fase de la respuesta en frecuencia'), grid

ylabel('Unidades Lineales')

%---

% Sistema con H1(z)

bH1=[1]; % Coeficientes b

% Coeficientes a

aH1=[1 -0.85*sqrt(2) 0.7225];

% Respuesta en frecuencia

[H1,W1]=freqz(bH1,aH1,512,'whole');

% Sistema con H2(z)

% Coeficientes b

bH2=[1 -0.85*sqrt(2) 0.7225];

aH2=[1]; % Coeficientes a

% Respuesta en frecuencia

[H2,W2]=freqz(bH2,aH2,512,'whole');

% Sistema con H3(z)

bH3=[1]; % Coeficientes b

% Coeficientes a

aH3=[1 -0.95*sqrt(2) 0.9025];

% Respuesta en frecuencia

[H3,W3]=freqz(bH3,aH3,512,'whole');

% Sistema con H4(z)

% Coeficientes b

bH4=[1 -0.95*sqrt(2) 0.9025];

aH4=[1]; % Coeficientes a

% Respuesta en frecuencia

[H4,W4]=freqz(bH4,aH4,512,'whole');

figure % Nueva figura

subplot(221)

% Módulo respuesta en frecuencia

plot(W1/(2*pi),abs(H1))

title('Módulo de H1'), grid

ylabel('Unidades Lineales')

subplot(223)

% Fase respuesta en frecuencia

plot(W1/(2*pi),angle(H1))

title('Fase de H1'), grid

ylabel('Unidades Lineales')

subplot(222)

% Módulo respuesta en frecuencia

plot(W2/(2*pi),abs(H2))

title('Módulo de H2'), grid

ylabel('Unidades Lineales')

subplot(224)

% Fase respuesta en frecuencia

plot(W2/(2*pi),angle(H2))

title('Fase de H2'), grid

ylabel('Unidades Lineales')

figure % Nueva figura

subplot(221)

% Módulo respuesta en frecuencia

plot(W3/(2*pi),abs(H3))

title('Módulo de H3'), grid

ylabel('Unidades Lineales')

subplot(223)

% Fase respuesta en frecuencia

plot(W3/(2*pi),angle(H3))

title('Fase de H3'), grid

ylabel('Unidades Lineales')

subplot(222)

% Módulo respuesta en frecuencia

plot(W4/(2*pi),abs(H4))

title('Módulo de H4'), grid

ylabel('Unidades Lineales')

subplot(224)

% Fase respuesta en frecuencia

plot(W4/(2*pi),angle(H4))

title('Fase de H4'), grid

ylabel('Unidades Lineales')

end

Ejercicio 3.2.3 RESPUESTA AL IMPULSO

Para ver el ejercicio y el gráfico seleccione la opción "Descargar" del menú superior

Ejercicio 3.2.4 DESCOMPOSICION EN FRACCIONES SIMPLES

c=1; % Control

if c==1

%--- Punto 1

a=[6 -5 1]; % Coeficientes a de y[n]

b=[1 -1]; % Coeficientes b de x[n]

[h,t]=impz(b,a,11); % h[n]

figure % Nueva figura

stem(t,h) % Respuesta al impulso del sistema

title('Respuesta al impulso del sistema'), grid

ylabel('Unidades Lineales')

%--- Punto 2

clc % Borra la pantalla

disp('Expansión en fracciones simples...')

% Expansión en fracciones simples

[r,p,k]=residue(b,a)

%--- Punto 3

n=0:10; % Tiempo discreto

% h[n]

h=-0.5*(0.5).^n+0.6667*(0.3333).^n;

figure % Nueva figura

stem(n,h) % Respuesta al impulso del sistema

title('Respuesta al impulso del sistema por fracciones simples'), grid

ylabel('Unidades Lineales')

end

Ejercicio 3.3 PROPIEDADES EN EL DOMINIO TRANSFORMADO

Ejercicio 3.3.1 FASE LINEAL

Para ver el ejercicio y el gráfico seleccione la opción "Descargar" del menú superior

PRACTICA 4. La DFT y la FFT. y PRACTCA 5.Diseño de filtros digitales.

>> help fft

FFT(X) calcula la Transformada de Fourier discreta (DFT) del vector X.

FFT Discrete Fourier transform.

FFT(X) is the discrete Fourier transform (DFT) of vector X. For

matrices, the FFT operation is applied to each column. For N-D

arrays, the FFT operation operates on the first non-singleton

dimension.

FFT(X,N) is the N-point FFT, padded with zeros if X has less

than N points and truncated if it has more.

FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the

dimension DIM.

For length N input vector x, the DFT is a length N vector X,

with elements

N

X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.

n=1

The inverse DFT (computed by IFFT) is given by

N

x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.

k=1

The relationship between the DFT and the Fourier coefficients a and b in

N/2

x(n) = a0 + sum a(k)*cos(2*pi*k*t(n)/(N*dt))+b(k)*sin(2*pi*k*t(n)/(N*dt))

k=1

is

a0 = X(1)/N, a(k) = 2*real(X(k+1))/N, b(k) = -2*imag(X(k+1))/N,

where x is a length N discrete signal sampled at times t with spacing dt.

See also IFFT, FFT2, IFFT2, FFTSHIFT.

>> help fftshift

Calcula la componente de la frecuencia cero en el centro del espectro.

FFTSHIFT Shift zero-frequency component to center of spectrum.

For vectors, FFTSHIFT(X) swaps the left and right halves of

X. For matrices, FFTSHIFT(X) swaps the first and third

quadrants and the second and fourth quadrants. For N-D

arrays, FFTSHIFT(X) swaps "half-spaces" of X along each

dimension.

FFTSHIFT(X,DIM) applies the FFTSHIFT operation along the

dimension DIM.

FFTSHIFT is useful for visualizing the Fourier transform with

the zero-frequency component in the middle of the spectrum.

See also IFFTSHIFT, FFT, FFT2, FFTN.

>> help angle

La función angle calcula la fase , en radianes, de una matriz compleja.

ANGLE Phase angle.

ANGLE(H) returns the phase angles, in radians, of a matrix with

complex elements.

See also ABS, UNWRAP.

>> help freqz

FREQZ Z calcula la transformada de la respuesta en frecuencia digital del filtro.

FREQZ Z-transform digital filter frequency response. When N is an

integer, [H,W] = FREQZ(B,A,N) returns the N-point frequency

vector W and the N-point complex frequency response vector G

of the filter B/A:

-1 -nb

jw B(z) b(1) + b(2)z + .... + b(nb+1)z

H(e) = ---- = ----------------------------

-1 -na

A(z) 1 + a(2)z + .... + a(na+1)z

given numerator and denominator coefficients in vectors B and

A. The frequency response is evaluated at N points equally

spaced around the upper half of the unit circle. To plot

magnitude and phase of a filter:

[h,w] = freqz(b,a,n);

mag = abs(h); phase = angle(h);

semilogy(w,mag), plot(w,phase)

FREQZ(B,A,N,'whole') uses N points around the whole unit circle.

FREQZ(B,A,W) returns the frequency response at frequencies

designated in vector W, normally between 0 and pi. (See

LOGSPACE to generate W). See also YULEWALK, FILTER, FFT,

INVFREQZ, and FREQS.

>> help fir1

B = FIR1(N, Wn) diseña un filtro digital de los "paso-bajo" de orden N y devuelve los coeficientes del filtro en el vector B de longitud N+1. La frecuencia de atajo Wn debe estar entre 0 < Wn < 1,0, donde 1,0 corresponde a la mitad del rango de la muestra. El filtro B es real y tiene fase lineal. El aumento normalizado del filtro en Wn es de -6 dB.

FIR1 FIR filter design using the window method.

B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter

and returns the filter coefficients in length N+1 vector B.

The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0

corresponding to half the sample rate. The filter B is real and

has linear phase. The normalized gain of the filter at Wn is

-6 dB.

>> help Kaiser

La función Kaiser() devuelve el número de elementos de la ventana Kaiser.

KAISER KAISER(N,beta) returns the BETA-valued N-point Kaiser window.

Ejercicio 3.3 PROPIEDADES EN EL DOMINIO TRANSFORMADO

Ejercicio 3.3.2 ESTABILDAD

Para ver el ejercicio y el gráfico seleccione la opción "Descargar" del menú superior

Ejercicio 3.3.3 SISTEMA PASO TODO

Para ver el ejercicio y el gráfico seleccione la opción "Descargar" del menú superior

Ejercicio 3.4.2 RESPUESTA EN FRECUENCIA DE POLOS

Para ver el ejercicio y el gráfico seleccione la opción "Descargar" del menú superior

Ejercicio 3.4.3 APLICACION: DISEÑO DE FILTROS EN HENDIDURA

Para ver el ejercicio y el gráfico seleccione la opción "Descargar" del menú superior

David Savall Climent

20 de mayo de 2003