Baixe Introdução ao Processamento Digital de Sinais no Octave e Matlab e outras Notas de aula em PDF para Sinais e Sistemas, somente na Docsity!
Processamento Digital de Sinais – Laboratório
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)
Manipulações Básicas
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. 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
» b = b’ Resulta em: b = 1 4 7 2 5 8 3 6 9 de novo. » b = b’ resulta em: b = 1 2 3 4 5 6 7 8 9 O vetor nada mais é do que um caso particular da matriz. O vetor linha é uma matriz de uma linha e várias colunas enquanto o vetor coluna é uma matriz de uma coluna e várias linhas. 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: Resulta em: i = 1 2 3 4 5 6 7 8 9 10 ou então: » x = 0:.1: Resulta em: x = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
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. 4 0 6 7 0.8 0. 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:
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^3 ). 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(jwn); 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)
» 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=152pi/32; na terceira tela. E w=4.5pi/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(4pi*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.3pina/20); » Nd=Na/ » nd=0:(Nd-1); » xd=xa(nd5+1); » Ni=Nd5; » 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)
plot(w,abs(H),ww,abs(X)/16) Vamos obter o sinal filtrado y de três maneiras distintas:
- Usando a equação recursiva: função filter() » y1=filter(b,a,x);
- Fazendo a convolução com a resposta ao impulso » d=[1 zeros(1,511)]; » h=filter(b,a,d); » y2=conv(x,h);
- 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
- 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.
- 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.1pin)+cos(0.9*pi);
- Levantar a resposta em freqûencia da função interp(). Basta aplicar um impulso e calcular a transformada da resposta impulsional.
- 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’.
- 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=2pi/32; w2=8pi/32; w3=12.5pi/32; x=sin(w1t)+sin(w2t)+sin(w3t)+0.5*rand(1,32); Usando três métodos: i) função filter() ii) função conv() iii) multiplicação por H(w)
