Prof. Eduardo Parente Ribeiro

Introdução ao Octave e Matlab

Introdução
O matlab é um ambiente computacional para manipulação de matrizes e visualização dos dados. A versão atual é a 5.3, porem a versão que estamos utilizando é a 4.2. O Octave, atualmente chamado de GNU Octave, visto que é distribuído gratuitamente sob a licença GNU (acompanha os fontes em C++) funciona em ambiente unix enquanto o Matlab roda tanto em Windows quanto Unix. O Octave utiliza uma notação idêntica ao matlab, diferindo apenas em pequenos aspectos. De maneira breve, podemos dizer que o Octave incorpora alguns conceitos que o Matlab não utiliza, como estrutura de dados (na ultima versão do Matlab isto já foi incorporado), e o Matlab possui a parte gráfica mais desenvolvida (O Octave utiliza o Gnuplot para a geração dos gráficos). A versão do Octave que utilizamos é a 2.0.14.
 

Obtendo Ajuda
O primeiro comando a se aprender, e o mais importante é o comando ‘help’. Ele mostra uma lista dos tópicos é que existe ajuda disponível. Para obter ajuda sobre o tópico pode se usar ‘help tópico’, ou então ‘help nome-da-função’. Exemplos:
» help
» help signal
» help inv

No Octave, o manual online pode ser acessado com ‘help –i’, ou ‘help –i tópico’. O programa ‘info’ é invocado e para navegar nesta programa pode-se usar os seguintes comandos:

Espaço: Próxima Pag.
DEL: Pag. Anterior
b: inicio da seção (beggining)
e: final da seção (ending)
TAB: pula p/ proxima referencia
ENTER: Entra na referencia
u: volta pro nivel anterior (up)
n: Entra no proximo tópico (next)
p: Entra no tópico anterior (previous)

Formando um vetor:
» a = [ 1 2 3 7 6 5]

Resulta em:
a =

1 2 3 7 6 5

Formando uma Matriz
» b = [1 2 3; 4 5 6; 7 8 9]

Resulta em :
b =

1 2 3
4 5 6
7 8 9

Cada elemento do vetor ou matriz pode ser uma expressão:
» c = [3 –4 2.3 sqrt(2) (4+5)/2]

Resulta em:
c =

3 -4 2.3 1.41 4.5

Transpondo um vetor, ou uma matriz:
Usa-se o operador ‘ (apóstrofe que é a aspa simples. Não confundir com o acento agudo nem o grave).
» c = c’

Resulta em:
c =

3
-4
2.3
1.41
4.5

» b = b’

Resulta em:
b =

1 4 7
2 5 8
3 6 9

de novo.
» b = b’

resulta em:
b =

Formando um vetor de contagem.
Trata-se de um vetor como outro qualquer porem pode ser gerado com uma notação simplificada. Por ex.:
» i=1:10

Resulta em:
i =

1 2 3 4 5 6 7 8 9 10

ou então:

» x = 0:.1:1

Resulta em:
x =

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Acessando os elementos de um vetor ou matriz:

aproveitando o vetor x acima temos:
» x(2)

resulta em

0.1

que é o segundo elemento do vetor x. A numeração sempre começa de 1 e não de zero como na linguagem "C’.

De forma análoga, para a matriz b definida anteriormente:

» b(2,3)

Resulta em:
6

que é o elemento da segunda linha terceira coluna. Sempre (linha, coluna).

Tambem pode-se se acessar uma sub-matriz:

» b(1:2,1:2)

Resulta em:
1 2
4 5

Isto é equivalente a b([1 2],[1 2]), ou seja estamos usando vetores para indexar uma matriz (ou outro vetor). Outro exemplo:

» b([1 3],[1 3])

resulta em:
1 3
7 9

Selecionando uma coluna inteira (o mesmo pode ser usado para linha)

» b(:,3)

resulta em:
3
6
9

Modificando um elemento de uma matriz ou vetor.

Da mesma forma que acessamos para leitura um elemento, o mesmo pode ser modificado, colocando-o a esquerda do operador de atribuição "=".

» b(2,2)=0

resulta em:
1 2 3
4 0 6
7 8 9

Pode-se tambem operar em sub-matriz, mas deve-se sempre observar a seguinte regra: O tamanho da matriz do lado esquerdo do operador "=" deve ser igual ao tamanho da matriz do lado direito (mesmo número de linhas e colunas).

» b([1 3],[2 3])=[0.2 0.3; 0.8 0.9]

Resulta em:
1 0.2 0.3
4 0    6
7 0.8 0.9

Tratando uma matriz como vetor.
Às vezes é interessante operar com todos os elementos (ou alguns elementos) de uma matriz como se fosse um vetor. Por exemplo para achar o maior elemento de um vetor usa-se a função max().

» max(c)

resulta em:
4.5

Quando se usa a função max() em uma matriz, o que é retornado é um vetor contendo os maiores elementos de cada coluna (vide help max).

» max(b)

resulta em:
7 0.8 6

Se quiser obter o maior elemento de toda a matriz, basta tratar a matriz como um vetor:

» max(b(:))

resulta em:
7

Salvando as variáveis
Todas as variáveis que foram criadas permanecem definidas no ambiente até que se finalize o matlab. Para ve-las utilize o comando:

» who

ou então

» whos

Para salva-las em um arquivo:

» save nome_do_arquivo

Para carregar este ambiente posteriormente numa nova sessão do matlab use:
» load nome_do_arquivo

Outros simbolos

O simbolo % representa um comentario, todo restante da linha é ignorado.

A linha de comando pode opcionalmente ser terminada com ponto e virgula ; desta forma o matlab não emite eco. Ele processa o comando e passa para a próxima linha.

Os simbolo * representa a multiplicação matricial (linha por coluna) entre dois operandos. Se a matriz for de apenas um elemento esta operação é a propria multiplicacão entre os números. Para se obter uma multiplicação elemento a elemento entre duas matrizes deve se usar a operação .* (ponto asterisco). Para divisão o equivalente é /, e ./ .

A operação de potência é representada por ^ (circunflexo). Ex. 2^3 (2³). Lembrando que a potencia tambem pode ser fracionária: ex. 8^(1/3) (raiz cubica de oito).

Edição da linha de comando
Para modificar um comando anterior que já tenha sido entrado, seja para corrijir um erro ou para criar um novo comando pode-se se ‘navegar’ pelos comandos anteriores utilizando as setas para cima e para baixo. Depois as setas a direita e a esquerda posicionam o cursor no ponto desejado, basta fazer as alterações e clicar ‘enter’ em qualquer lugar da linha que esta será reprocessada.

Quando se desejar resgatar um linha que já foi digitada há muito, pode se utilizar a procura rápida bastando para tanto digitar no inicio da linha vazia as iniciais da linha desejada, e em seguida apertar seta para cima até a linha desejada aparecer.

Processamento de Sinais
Agora vamos estudar o toolbox de processamento de sinais. Os comandos contidos neste pacote podem ser listados com:

» help signal

Convolução
A função conv() calcula a convolução. Na verdade, este é um comando interno do matlab, para proporcionar uma alta eficiência no calculo.

» a=conv([1 1 1 1],[1 1 1 1 ])

Resulta em:

1 2 3 4 3 2 1

Vamos criar um sinal x=[1 1 1 1 1] e uma resposta impulsional h=[1 0.5 0.25]

A saída do sistema linear e invariante no tempo (LTI) é dada pela convulção da entrada com a resposta imulsional:

» y = conv(x,h)

Visualizando a exponencial complexa

» n=0:99; % define um vetor n indo de 0 a 99

» w=0.3; %define a frequência digital

» d=exp(j*w*n);

Para visualizar utilize o comando subplot:

» subplot(4,1,1) % divide a tela em 4 linhas x 1 coluna e seleciona o primeiro grafico

» plot(real(d))

» subplot(4,1,2) % seleciona o segundo grafico

» plot(imag(d))

» subplot(4,1,3)

» plot(abs(d))

» subplot(4,1,4)

» plot(angle(d)) % ou então plot(unwrap(angle(d)))

verifique agora outras frequências: w=0.6; w=0.3+2*pi;
 
 

Visualizando a FFT
defina um sinal x:

» t=0:31;

» w=2*pi/32;

» x=sin(w*t);

» X=fft(x);

» subplot(2,1,1); stem(t,x)

» subplot(2,1,2); stem(t,abs(X))

Agora em crie outra tela para plotar a frequência w=4*pi/32;

» figure(2)

» w=4*pi/32;

» x=sin(w*t);

» X=fft(x);

» subplot(2,1,1); stem(t,x)

» subplot(2,1,2); stem(t,abs(X))

Interprete o resultado. Porque o impulso que aparece em frequência caminhou para a direita?

Experiemente outras frequência: w=15*2*pi/32; na terceira tela. E w=4.5*pi/32 na quarta tela.

Por que o espectro da ultima senoide não é apenas um impulso, mas um impulso ‘contaminado’ ?

Efeito do Janelamento
Agora complete o sinal x com 64 zeros e plote a FFT :

» n=0:31;

» x=[sin(4*pi*n/32) zeros(1,64)];

» nn=0:95;

» X=fft(x);

» subplot(2,1,1); plot(nn,x)

» subplot(2,1,2); plot(nn,abs(X), nn,abs(X),’o’)

Interprete o resultado. Qual é a transformada de Fourier da janela retangular ? O que diz o teorema da convolução ?

Transformada das Janelas
Verifique a transformada de fourier de algumas janelas e plote-as no mesmo gráfico:

» j1=[zeros(1,32) hamming(64)’ zeros(1,32)];

» j2=[zeros(1,32) hanning(64)’ zeros(1,32)];

» j3=[zeros(1,32) ones(1,64) zeros(1,32);

» J1=fft(j1); J2=fft(j2); J3=fft(j3);

» n=0:127;

» subplot(2,1,1)

» plot(n,j1,n,j2,n,j3);

» subplot(2,1,2)

» plot(n,abs(J1),n,abs(J2),n,abs(J3));

Depois, dê um zoom, visualizando os pontos de 1 a 10:

» n=0:9;

» plot(n,abs(J1(1:10)), n,abs(J2(1:10)), n,abs(J3(1:10)));
 
 

Decimação e Interpolação
Vamos imaginar um sinal que foi amostrado com uma taxa bem superior a sua frequência máxima e que foram coletados Na=160 pontos. Depois deseja-se decimar o sinal (M=5) e em seguida interpolar (L=5).

» Na=160;

» na=0:(Na-1);

» xa=sin(2.3*pi*na/20);

» Nd=Na/5

» nd=0:(Nd-1);

» xd=xa(nd*5+1);

» Ni=Nd*5;

» ni=0:(Ni-1);

» xz(nd*5+1)=xd(nd+1);

» [b a]=ellip(4,1,40,0.2);

» xi=filter(b,a,xz);

» subplot(4,1,1)

» stem(xa)

» subplot(4,1,2)

» stem(xd)

» subplot(4,1,3)

» stem(xz)

» subplot(4,1,4)

» stem(xi)

Linearidade
Verifique a propriedade de linearidade da fft, somando duas senoides (32 pontos) e verificando o espectro. w1=2*pi/32 e w2=10*pi/32; Ex. x3=x2+x1; X3=fft(x3); Depois faça o mesmo, mas usando w2=10.5*pi/32.

Interprete o gráfico.
 
 

Filtragem Digital
veja a função filter():

» help filter

Considere o sinal x1 uma senoide de baixa frequência e o sinal x2 uma senoide de alta frequência e o sinal x sendo a soma destes dois:

» n=0:31;

» x1=sin(4*pi*n/32);

» x2=sin(12*pi*n/32);

» x=x1 + x2;

Visualize o sinal x e seu espectro:

» subplot(2,1,1); plot(x)

» X=fft(x);

» subplot(2,1,1); stem(abs(X))

Agora considere obtenha os coeficientes de um filtro elíptico de quarta ordem com frequência normalizada de corte igual a 0,25.

» [b a]=ellip(4,1,40,0.25)

Visualize a resposta em frequência do filtro e o espectro do sinal x no mesmo eixo:

[H w]=freqz(b,a);

ww=n*2*pi/32;

plot(w,abs(H),ww,abs(X)/16)

Vamos obter o sinal filtrado y de três maneiras distintas:

1) Usando a equação recursiva: função filter()

» y1=filter(b,a,x);

2) Fazendo a convolução com a resposta ao impulso

» d=[1 zeros(1,511)];

» h=filter(b,a,d);

» y2=conv(x,h);

3) A partir da multiplicação das transformadas

» H=fft(h);

» xx=[x zeros(1,512-32)];

» XX=fft(xx);

» Y3=XX .* H;

» y3=real(ifft(Y3));

Agora visualize as três respostas no mesmo eixo:

» plot(n,y1,n,y2(n+1),n,y3(n+1))

O resultado é coincidente e apenas uma curva aparece. Se quiser ter certeza, experimente:

» plot(n,y1,n,y2(n+1)+0.1,n,y3(n+1)+0.2)

Outras Funções importantes:
zplane():traça o diagrama de pólos e zeros

roots():calcula as raízes do polinômio

conv2(): convolução em 2 dimensões

fft2(): fft em 2 dimensões

filtdemo(): demonstração de proj. filtro
 
 

Exercícios
1) Calcule a Função de transferencia do filtro definido por b=[0.0186 0.0743 0.1114 0.0743 0.0186];

a=[1.0000 -1.5704 1.2756 -0.4844 0.0762]; por dois metodos diferentes:

i) usando a função freqz()

ii) aplicando um impulso e tirando a fft

Plote as duas respostas no mesmo eixo.
 

2) Plotar no mesmo eixo a resposta em frequência de quatro filtros de ordem 8 e freq. de corte=0.5. (butter, cheby1, cheby2, ellip). Visualizar na mesma tela (com subplots) os diagramas de pólos e zeros. Filtrar o sinal

» n=0:63;

» x=sin(0.1*pi*n)+cos(0.9*pi);
 
 

3) Levantar a resposta em freqûencia da função interp(). Basta aplicar um impulso e calcular a transformada da resposta impulsional.
 
 

4) Comparar a convolução circular e linear entre dois sinais: x1=[1 1 1 1] e x2=[5 4 3 2]. Sugestão: a função conv() calcula a convolução. Pelo teorema da convolução circular o produto das transformadas equivale a convolução circular no tempo. Agora preencha com zeros de tal maneira que ‘conv. circ.’ = ‘conv. linear’.
 
 

5) Utilizando o filtro dado pelos coeficientes b=[0.0186 0.0743 0.1114 0.0743 0.0186];

a=[1.0000 -1.5704 1.2756 -0.4844 0.0762]; realize a filtragem do sinal:

t=0:31;

w1=2*pi/32; w2=8*pi/32; w3=12.5*pi/32;

x=sin(w1*t)+sin(w2*t)+sin(w3*t)+0.5*rand(1,32);

Usando três métodos:

i) função filter()

ii) função conv()

iii) multiplicação por H(w)
 
 

Funções
Crie um função com o seu nome, que retorne a raiz quarta de um numero.
ex. crie com um editor de texto o arquivo ‘nome.m’ contendo:

Function a=nome(b)

a = b ^ (1/4);

depois execute:

» d = nome(16)

Resulta em:

2

Desta maneira os comandos da linguagem matlab podem ser expandido pelo usuário. Na verdade, muitas das funções que foram usadas até agora, estão escritas na linguagem do matlab. O matlab é composto de umas poucas funções internas e varias funções externas que são arquivos .m.
 

SCRIPTS
A diferença entre um script e uma função é que o segundo não possui o cabeçalho ‘function c=nome(a)’ . O script é apenas um arquivo que agrupa uma sequência de comandos para serem executados de uma vez.

Bibliografia

Matlab 4.0 User’s Guide (agosto 1992)
Octave 2.0.14, User's Documentation.