20

Click here to load reader

Tutorial MatLab

Embed Size (px)

Citation preview

Page 1: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

1

MATLAB1) Introdução

O aplicativo Matlab é uma das ferramentas mais úteis para as áreas de processamento desinais, controle, e outras aplicações. O “Mat” do nome Matlab está relacionado com “Matriz”, ouseja, é um laboratório que permite a manipulação eficiente de matrizes. Pode-se dizer também que oMatlab permite também a manipulação eficiente de vetores, já que um vetor de tamanho N pode servisto como uma matriz de Nx1.

Uma das melhores formas de se familiarizar com o Matlab é através do comando intro, queprovê uma turnê pelas características básicas do programa.

Exercício 1

Use o comando intro e, com atenção, procure entender todos os exemplos fornecidos.

A aplicação predominante no presente curso será o processamento digital se sinais. Paraisto, preparamos algumas técnicas em Matlab que facilitarão o uso do aplicativo para esta tarefa.

2) Geração de Sinais AmostradosComo primeiro exemplo, geraremos a seguinte função senoidal no Matlab:

f(t)=sin(2*pi*fo*t)

Onde fo=1 Hz.

A função é mostrada na seguinte figura:

0 0.2 0.4 0.6 0.8 1-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Tempo (seg.)

f(t)

Senoide de 1 Hz

Page 2: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

2

Para permitir o processamento de sinais por computador, o sinal analógico deve serdigitalizado por um conversor analógico para digital (A/D). Neste processo o sinal é amostrado emvários instantes, em intervalos de tempo constante. O intervalo de tempo é o período de amostragem(T), e o recíproco de T é a freqüência de amostragem (fs= 1/T). A unidade de tempo do período deamostragem em geral é segundos, e a de Freqüência é em Hz. 100 Hz significa que 100 amostras dosinal são tomadas em 1 segundo.

Pode-se simular funções amostradas no Matlab. O exemplo abaixo mostra como se podesimular a amostragem de uma senóide no Matlab:

Exercício 2

Siga os passos descritos para simular a amostragem de uma onda senoidal no seguinteexemplo. Digite os comandos na Matlab, e observe os efeitos.

1) Crie o eixo do tempo:t=0:0.05:1;

O comando acima cria uma variável (vetor) com os elementos 0, 0.05, 0.1, ..., 1. O ponto evírgula evita que a variável seja exibida na tela. Tente rodar o comando acima e observe o queacontece. Caso você tenha errado, e o vetor inteiro comece a aparecer na tela, você pode parar oprocesso usando CTRL-C.

O comando size(t) exibe o tamanho da matriz t, e o comando length(t) mostrará o comprimento dovetor. Experimente estes comandos.

2) Crie a função:f=sin(2*pi*t);

O comando acima cria um vetor cujas componentes são a função senoidal calculada a cada valor dovetor t. Assim, no Matlab, você pode calcular o seno de um vetor, ou mesmo de uma matriz.Execute este comando.

3) Plotar a função:plot(t,f)

Este comando faz a plotagem do sinal, tendo f no eixo y, e o correspondente t no eixo x.Execute este comando.

4) Coloque o Titulo e as variáveis dos eixos x e y:title(‘Funcao Senoidal’)xlabel(‘Tempo (segundos)’)ylabel(‘Amplitude (Volts)’)

Page 3: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

3

Experimente os comandos acima.

Tente também os seguintes comandos:plot(f,t)

plot(f)

plot(t)

O que aconteceu no primeiro comando? No segundo comando, a senóide foi plotada, mas oeixo y mostrou simplesmente o número da amostra, mas não o tempo correspondente a cada abcissay. Explique o que acontece no terceiro exemplo.

Tente também as seguintes opções:plot(t,f,’*’)plot(t,f,’b*’)plot(t,f,’.’)plot(t,f,’c.’)plot(t,f,’-’)plot(t,f,’y-’)

Observe as cores e características dos gráficos. Agora digite help plot, e leiacuidadosamente as descrições da função plot.

Agora digite a seguinte função:grid

Note o efeito da função, e digite help grid, para ver as características desta função.

Agora gere uma segunda função:f2=sin(2*pi*2*t)

E plote as funções f e f2 no mesmo gráfico:plot(t,f,t,f2)

Verifique as cores dos gráficos. Agora digite:plot(t,f,’b*’,t,f2,’c-.’)

E verifique os efeitos.

Agora tente os seguintes comandos:stem(t,f)

stem(t,f2)

Page 4: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

4

O que estes comandos fazem? Digite help stem, e estude o funcionamento desta função.Tente mudar as cores dos traçados desta função. Note que a função não é tão flexível quanto a coresdo traçado quanto a função plot.

3) Geração de sinais mais complexos

Digite os seguintes vetores:A=[0 1 2 3 4 5 6];

B=[0 2 4 6 8 10 12];

Tente multiplicar A*B. Note que o Matlab não realizará tal operação. O fato é que o Matlabsempre interpretará o vetor como uma matriz, e o programa sempre interpreta uma multiplicaçãocomo uma multiplicação de matrizes. Entretanto, há uma opção que permite ao usuário fazer umamultiplicação ponto a ponto de dois vetores ou matrizes (isto é o primeiro elemento do resultado é oproduto dos primeiros elementos dos fatores, o segundo elemento do resultado é o produto dossegundos elementos dos fatores, e assim por diante). Tente digitar:

C=A.*B

Note que o vetor C é o resultado da multiplicação ponto a ponto do dos vetores A e B.

Exercício 3

Queremos plotar a seguinte função no intervalo 0<x<4:

y(x)=(x-1)(x-2)

Para isto, criamos primeiro o vetor representando a variável x, e em seguida a função emquestão:

x=0:0.05:4;

y=(x-1).*(x-2);

plot(x,y)

Note o comando “.*”. Cada expressão entre parênteses é um novo vetor, e os vetores sãomultiplicados ponto a ponto, resultando no efeito desejado.

Agora geraremos a função:

y(x)=x3-6x2+11x-6

Page 5: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

5

Para isto, criamos primeiro o vetor representando a variável x, e em seguida a função emquestão:

x=0:0.05:4;

y=x.^3-6*x.^2+11*x-6;

plot(x,y)

Note que o operador “.^” é também um operador ponto a ponto, que eleve cada elemento damatriz à potência desejada. Tente digitar x^2, e veja que este comando não funcionará. Tentedescobrir porque (lembre-se que o Matlab sempre assume operações matriciais).

PROBLEMAS

1) Gere e plot a função f(t)=exp(-t)sin(2*5*pi*t), no intervalo 0<x<6. Use várias escolhas deintervalo para o eixo x: 0.1 sec., 0.05 sec., etc., e escolha o que tem melhor aparência.

2) Ache graficamente as raízes da função y=x3+2x2-x+1. (Dica: procure, por tentativa e erro ointervalo, a o período de amostragem do eixo x, realizando sucessivas plotagens das tentativas).

3) Plote a seguinte função: y=sin(x)*cos(2x)/(2+sin(x)). Escolha um eixo x apropriado (Não se esqueçaque as multiplicações e divisões devem ser operações ponto a ponto).

4) Gere a função f(t)=sin(2*pi*10*t), no intervalo 0<x<10. Qual é a freqüência desta senóide? Useas seguintes escolhas de freqüência de amostragem: 5 Hz, 10 Hz, 20 Hz, 40 Hz, e 200 Hz. Use afunção stem. Observe bem o efeito, pois será analisado na próxima aula.

5) Gere a função f(t)=exp(-t)sin(2*5*pi*t), e depois plote um histograma da mesma. Dica: use afunção intro, e observe o uso da função bar. Em seguida digite help bar, para analisar aquelafunção.

Page 6: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

6

EXERCÍCIOS COM O MATLAB

Nesta seção, veremos algumas funções extras do Matlab. Primeiramente, consideremos oseguinte trecho a seguir. Analise atentamente o programa abaixo (consulte o help para as funçõesque você não conhece):

% exemplo: geracão de uma forma de onda ruidosat=0:0.01:pi; %linha 1y=2*sin(5*t); %linha 2figure(GCF) %linha 3plot(t,y) %linha 4pause %linha 5ruido=randn(1,length(t)); %linha 6plot(t,ruido) %linha 7pause %linha 8y_ruido=y+ruido; %linha 9plot(t,y_ruido); %linha 10

Execute o programa. A linha 1 cria o vetor do tempo de amostragem. A linha 2 cria umafunção senoidal (qual é a freqüência?). A linha três faz com que a janela com o gráfico sejacolocada à frente de todas as outras. A linha 4 plota uma senóide, como mostrado abaixo:

0 0.5 1 1.5 2 2.5 3 3.5-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

O comando pause, na linha 5, causa uma pausa no programa. Para continuar o programa, ousuário deve pressionar <Enter>. [Pergunta: o que acontece se este comando for substituído porpause(2)? Dica: use o comando help].

A linha 6 utiliza o comando randn (utilize o help para aprender sobre este comando). Estecomando gera uma seqüência de números aleatórios com distribuição normal, com uma linha etantas colunas quanto for o comprimento do vetor t. O comando lenght(t) fornece o número deelementos de t. Observe a plotagem do sinal de ruído:

Page 7: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

7

0 0.5 1 1.5 2 2.5 3 3.5-3

-2

-1

0

1

2

3

Em seguida, a linha 9 calcula uma nova variável y_noise, que é o sinal senoidal mais oruído. Observe o resultado:

0 0.5 1 1.5 2 2.5 3 3.5-5

-4

-3

-2

-1

0

1

2

3

4

Neste exemplo, você aprendeu como gerar um sinal com ruído.Agora, vejamos outros comandos. Para isto considere o seguinte código:

% Exemplo 2A=[1 2 3; 4 5 6]B=[3 4 ; 3 6 ; 7 8]C=[1 2 3 ; 4 5 6 ; 7 8 9];D=[5 6 ; 7 8 ; 9 10];pausesize(A)size(B)pause

Page 8: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

8

length(A)length(B)

Após executar o código acima, realize manualmente, e em seguida usando o Matlab, asseguintes operações (indique quando a operação não é possível):

A*B=

B*A=

A*C=

A^2=

A.^2=

C^2=

C.^2=

A’=

C*inv(C)=

[B B B]=

[A;A;A]

Veremos agora alguns comandos estatísticos. Considere o seguinte trecho de código:

% exemplo 3

A=[ 1 2 3 4 5 6 7] %linha 1B=mean(A) %linha 2C=std(A) %linha 3Pause %linha 4RUIDOSO=randn(1,5000); %linha 5figure(gcf) %linha 6plot(RUIDOSO) %linha 7title('sinal ruidoso') %linha 8pause %linha 9MEDIA_RUIDOSO=mean(RUIDOSO) %linha 10DP_RUIDOSO=std(RUIDOSO) %linha 11Pause %linha 12MAISRUIDO=5*randn(1,5000); %linha 13MEDIA=mean(MAISRUIDO) %linha 14DP=std(MAISRUIDO) %linha 15

Page 9: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

9

Use o comando help para estudar os comandos mean, e std. Execute o programa. Tomecuidado de só incluir o ponto-e-vírgula quando necessário. Estes comandos calculam a média e odesvio padrão dos vetores em questão (relembre na literatura o que estes termos significam). Nalinha 2, o resultado, B, será a média dos elementos do vetor A Na linha 3, o resultado C será odesvio padrão dos elementos do vetor A. Calcule A e B manualmente, e cheque os resultados comos resultados do Matlab.

Na linha 5 um sinal de ruído normal é gerado, e é plotado na linha 6. A média e o desviopadrão são calculados e mostrados nas linhas 10 e 11. Note que a função randn gera ruído aleatóriocom média 0, e desvio padrão 1 (na verdade, próximos de 0 e 1). Para obter um ruído com maiorintensidade (maior desvio padrão), basta multiplicar o ruído pelo desvio padrão desejado, como éfeito nas linhas 13, 14 e 15. Analise com cuidado estas três últimas linhas.

EXERCÍCIOS

1) Usando o comando help, aprenda a lidar com os comandos max, min, mean, median, sort.2) Usando o comando help, aprenda a usar os comandos fix, floor, ceil, round, mod,

rem, sign.3) Gere uma senóide com freqüência de 50 Hz e amplitude 3V, amostrada a uma freqüência de

amostragem de 1 kHz. Adicione a esta senóide ruído gaussiano com desvio padrão de 2V. Ploteo sinal resultante. Determine os valores máximo e mínimo deste sinal. Não se esqueça de incluirtítulo e unidades. Procure plotar o sinal de forma que suas características fiquem claras ao leitor(deve ficar claro que o sinal é uma senóide com ruído).

4) Usando o comando help square, aprenda como gerar uma onda quadrada. Gere uma ondaquadrada de 10 Hz, amostrada a uma freqüência de 1 KHz.

Page 10: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

10

REVISÃO: NÚMEROS COMPLEXOS

1) Representação dos Números Complexos

Números complexos são uma estrutura criada como uma extensão dos números reais. Àprimeira vista, esta estrutura parece ser artificial e inútil. Entretanto, na prática suas aplicações sãoinúmeras. Este capítulo apresentará uma revisão resumida do tópico, devida à sua importância nasséries de Fourier.

Um número complexo é um par ordenado de números reais (R,I), e é muito usualmenterepresentado na seguinte forma:

C=R+iI, ou C=R+jI

Sendo que matemáticos preferem o uso do i, enquanto os engenheiros preferem o j. R échamado parte real, e I é chamado parte imaginária do número complexo. Note que se I é zero, onúmero se transforma em puramente complexo, e se R é zero, o número é puramente imaginário.

Os números complexos seguem certas regras algébricas, que a princípio parecem estranhas,mas que permitem a extensão de várias técnicas matemáticas, bem como possibilitam um grandenúmeros de aplicações. É muito comum se representar os números complexos em um planoordenado x,y, da forma ilustrada na Figura 1.

(x,y)

X=M sin(α)

M cos(α)

M

Real

imaginario

a=tan-1(y/x)

Figura 1 - Representação gráfica de um número complexo

O módulo ou magnitude M de um número complexo é representado geometricamente pelareta ligando a origem à coordenado do par ordenado, e é definido como:

M=sqrt(x2+y2)

O ângulo ou fase α, do número complexo é o ângulo entre a reta mencionada e o eixo x, e édefinido como:

α=tan-1(I/R)

É fácil mostrar, observando a geometria da Figura 1 que:

Page 11: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

11

R=M cos(α)

I=M sin(α)

Muitas vezes o número complexo é representado por seu módulo e sua fase, da seguintemaneira:

M ∠ α

Exercício 1:

Considere os seguintes números:

A=1+j2

B=-2+j3

C=4-j1

D=-2-j2

1) Plote manualmente cada um do números, e enumere em que quadrante o número está (1o., 2o.,3o., ou 4o.) e depois plote os mesmo usando o Matlab (plot(A) plota o número complexo, caso omesmo realmente seja complexo).

2) Para cada número, calcule manualmente o módulo M, e depois faça os cálculos com o Matlab(comando abs(A)). Confira os resultados com os cálculos manuais.

3) Calcule manualmente os ângulos α (atente para os quadrantes) e em seguida repita os cálculosem Matlab (comando atan(x/y)). Tente também a função atan2(x/y), e observe que a mesmaleva em consideração o quadrante em que o número está.

4) A partir dos valores de M e α calculados, determine x e y. Faça o mesmo em Matlab (Exemplo,x=M*cos(alfa), y=M*sin(alfa)).

5) Digite A=1+j2, e em seguida, digite:

real(A) imag(A)

O que acontece? Qual é a função dos comandos acima?1

6) Digite o seguinte código:

A=[ 1+j*2 ; 2+j*3 ; 4-j*1 ; 2-j*2 ]plot(A)M=abs(A)gridalfa1=atan(real(A)./imag(A))

Page 12: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

12

alfa2=atan2(real(A),imag(A))Areal=M.*cos(alfa2)Aimag=M.*sin(alfa2)

Tente entender todos os comandos. Note que os exercícios de 1 a 4 foram implementadosnesta linha. Note o caráter matricial (ou vetorial) das variáveis (preste atenção nos pontos do .*, e ./.

6) Converta os seguintes números para a forma R+jI:

2 π/3

3 60o

4 π/2

2) Álgebra dos Números Complexos

Os números complexos são sujeitos a regras de álgebra, algumas delas semelhantes àálgebra dos números reais, e algumas diferentes. Quando os números complexos são puramentereais, a álgebra simplesmente se transforma na álgebra dos números reais.

Suponha os números complexos A e B:

A=Ar+ j Ai

B=Br+ j Bi

Estes números complexos são sujeitos às seguintes regras:

Soma: A+B=(Ar + Br) + j ( Ai + Bi)

Valor de j: j=i=sqrt(-1)

Quadrado de i: i2=-1

Multiplicação: A . B=( Ar+ j Ai).( Br+ j Bi)= Ar Br +j Ai Br+j BrAi+j2AiBi

= (ArBr - AiBi)+j (ArBi+AiBr)

Raiz n-ésima (fórmula de de Moivre):

a+2kpA1/n= n |A| ∠ n

Complexo conjugado de R+jI : R-jI

Multiplicação de um número complexo por seu conjugado:

(A+jB)(A-jB)=A2+B2 (note a eliminação da parte complexa)

Page 13: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

13

Eliminação de partes complexas no denominador:

A+jB = A+jB C-jD = (AC-BD) + j(AD-BC)C+jD C+jD C-jD C2+D2

Exercício 2

Suponha:

M= 1+j2

N= 3+j4

Calcule primeiro manualmente, e depois em Matlab. Dê os resultados na forma R+jI :

1) M+N

2) MN

3) M/N

4) M2

5) sqrt(M)

6) M1/3

7) M1/5

8) Determine todas as raízes cúbicas de 8 (Fórmula de de Moivre ou Matlab)

9) Eleve todas as raízes encontradas em na questão 8 ao cubo.

3) Fórmula de Euler:

Duas fórmulas fazem os números complexos muito úteis:

ejω = cos ω + j sin ω

e-jω = cos ω - j sin ω

Ou equivalentemente:

cos ω = ejω + e-jω

2sin ω = ejω − e-jω

j2Note que

cos ω = Re (ejω)

sin ω = Im (ejω)

Page 14: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

14

Exercício 3

1) Defina uma variável de tempo com 0<t<5, amostrada em intervalos de 0.01 segundos, e calcule afunção

f(t)=e2jπt

Usando seus conhecimentos, tente prever o comportamento desta função. Note que a funçãoé formada por uma parte real senoidal, e uma parte imaginária cossenoidal. Plote a parte real e aparte imaginária na Função no Matlab:

ParteReal=real(f);ParteImag=imag(f);subplot(2,1,1)plot(t,ParteReal)subplot(2,1,2)plot(t,ParteImag)

2) A função também pode ser dividida em uma parte função módulo e fase. Tente prever a formadas funções módulo e fase da função seus conhecimentos. Use o Matlab para plotar estasfunções:

Módulo=abs(f);Fase=atan2(imag(f),real(f));subplot(2,1,1)plot(t,Módulo)subplot(2,1,2)plot(t,Fase)

3) Multiplique a função complexa f, calculada no exercício anterior, pelo número complexo 1+j1, echame a nova função de f2. Plote as partes reais de f1 e f2 no mesmo gráfico. Observe o efeito damultiplicação por um número complexo. Por que este efeito ocorre? Preste bastante atenção aesta questão, pois é central para o entendimento posterior dos conceitos de resposta emfreqüência, e FFT.

4) Gere a função f(t)=sin(2πt). Multiplique esta função por 1+j1, gerando f2. Plote no mesmográfico as partes imaginárias de f e f2. Note os efeitos, e explique-os.

A FFT e a DFT

O Matlab permite o cálculo fácil da DFT e da FFT. Se tivermos um vetor A, de n elementos,a DFT do mesmo é:

Page 15: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

15

FFTdeA=fft(A);

E para a FFT inversa, é usada a seguinte forma:IFFTdeA=ifft(FFTdeA);

Caso o número de elementos de A seja da forma 2k, o Matlab usa o algoritmo da FFT (maisrápido). Caso contrário, é usado o algoritmo da DFT. Considere o seguinte exemplo:t=(0:.01:.99);f=sin(2*pi*2*t);F=fft(f);subplot(2,1,1)plot(abs(F),'ko')subplot(2,1,2)plot(angle(F),'ko')

O resultado é mostrado a seguir:

0 10 20 30 40 50 60 70 80 90 1000

10

20

30

40

50

0 10 20 30 40 50 60 70 80 90 100-4

-2

0

2

4

Note que foram plotados o ângulo e a fase da FFT. Isto porque o resultado é um vetor comnúmeros complexos. Tente plotar o vetor sem a função abs. Explique a razão para o resultadoaparentemente absurdo.

Neste exemplo, é calculada a DFT de 1 segundo de amostragem de uma senóide com 2 Hz.

Page 16: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

16

Note que o resultado da fase foi estranho. A razão é que a fase é o resultado de divisões denúmeros muito pequenos por números muito pequenos, resultando num valor não significativo.

Substitua agora a primeira linha port=(0:.01:1.3);

e execute o programa. O resultado é mostrado a seguir:

0 20 40 60 80 100 120 1400

10

20

30

40

50

0 20 40 60 80 100 120 140-4

-2

0

2

4

Nota-se que neste último caso ocorreu o fenômeno do “leakage”. Este fenômeno nãoocorreu no primeiro caso porque as amostras completavam um ciclo completo (a menos de umaamostra). Para que não ocorra leakage, é necessário que esta precaução seja tomada.

Voltemos agora ao primeiro exemplo. Nota-se que o eixo dos x não corresponde ãfreqüência e sim ao número da amostra. Para que seja mostrada a amostra, é necessário que amesma seja calculada. Para tal, suponhamos que a freqüência de amostragem seja fs e que o sinalcontenha N amostras. Neste caso, a primeira amostra corresponderia à freqüência zero, e a amostraN (que não existe, já que as amostras vão de 0 a N-1) corresponderia à freqüência de amostragem.Assim, a freqüência de cada amostra, na DFT resultante corresponderia a fs/N. Por exemplo, a

Page 17: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

17

primeira amostra corresponderia à freqüência zero, a Segunda a fsN, a terceira a 2fs/N, e a n-ésima a(n-1)fs/N.

Outra observação importante é com respeito à amplitude. A amplitude da senóide noexemplo é de 1 V, e a amplitude da DFT foi de 50. Isto se deve ao fato de que o resultado daFFT é sempre multiplicado pelo número de amostras. Assim, para se Ter uma estimativaconfiável das componentes espectrais, é necessário que se divida o resultado por N.

Tendo essas considerações em mente, execute o seguinte código:

clear % Limpa variáveisclf % Limpa gráficost=(0:.01:.99); % Gera instantes das amostrasfs=100; % Freqüência de amostragemN=length(t); % Determina o número de amostrasf=1+sin(2*pi*10*t);% Amostra o sinal formado por

% uma cte + senóideF=fft(f); % Calcula a DFTFesc=F/N; % Escalona para o valor corretofreq=(0:N-1)*fs/N; % Calcula eixo das freqüênciasplot(freq,abs(Fesc),’ko’) % Plota o módulo da DFTxlabel(‘Freqüência ’) % Label do eixo xylabel(‘Amplitude’) % Label do eixo y

O resultado é mostrado a seguir:

Page 18: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

18

0 10 20 30 40 50 60 70 80 90 1000

0.2

0.4

0.6

0.8

1

1.2

1.4

Freqüência

Am

plitu

de

FILTROS DIGITAIS

O Matlab apresenta muitas funções que facilitam a implementação de filtros digitais.Observe o help da função filter.

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

When X is a matrix, FILTER operates on the columns of X. When X

Page 19: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

19

is an N-D array, FILTER operates along the first non-singleton dimension.

[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, FILTFILT (in the Signal Processing Toolbox).

Assim, a função tem, em geral, três argumentos: os vetores A e B, correspondentes aoscoeficientes dos filtros, e o vetor de entrada. O resultado da função é o vetor de saída.

Antes de fazermos um exemplo com a função filter, vamos observar que o Matlab possuivárias ferramentas que permitem a determinação dos coeficientes A e B. Um exemplo é a funçãobutter. Observe o help desta função:

BUTTER Butterworth digital and analog filter design. [B,A] = BUTTER(N,Wn) designs an N'th order lowpass digital Butterworth filter and returns the filter coefficients in length N+1 vectors B and A. The cut-off frequency Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample rate.

If Wn is a two-element vector, Wn = [W1 W2], BUTTER returns an order 2N bandpass filter with passband W1 < W < W2. [B,A] = BUTTER(N,Wn,'high') designs a highpass filter. [B,A] = BUTTER(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].

When used with three left-hand arguments, as in [Z,P,K] = BUTTER(...), the zeros and poles are returned in length N column vectors Z and P, and the gain in scalar K.

When used with four left-hand arguments, as in [A,B,C,D] = BUTTER(...), state-space matrices are returned.

BUTTER(N,Wn,'s'), BUTTER(N,Wn,'high','s') and BUTTER(N,Wn,'stop','s') design analog Butterworth filters. In this case, Wn can be bigger than 1.0.

See also BUTTORD, BESSELF, CHEBY1, CHEBY2, ELLIP, FREQZ, FILTER.

Page 20: Tutorial MatLab

Tópicos em Processamento de SinaisDepartamento de Engenharia Elétrica

Faculdade de TecnologiaUniversidade de Brasília

20

Vê-se que esta função determina os coeficientes de um filtro de Butterworth de ordem n, com umafreqüência de corte fc. A freqüência de corte é indicada de forma que o valor 1 corresponda àmetade da freqüência de amostragem fs. Por exemplo, se fs=100Hz, 0.1 corresponderia a umafreqüência de corte de 5 Hz.

Como exemplo, será filtrada uma senóide de 15 Hz, amostrada a 100 Hz, por um filtro comfc=10 Hz.

clearclft=(0:.01:1);f=sin(2*pi*5*t); % amostragem do sinalsubplot(2,1,1)plot(t,f) % Sinal de entrada[B,A]=butter(1,0.1); % Calcula coeficientes do filtro de 2ª ordem,

% fs/2=50Hz e fo=5Hz, entao Wn=5/50=0.1saida=filter(B,A,f); % Filtra o sinal. Resultando no vetor de

% saída.subplot(2,1,2)plot(t,saida)

Execute este código, e explique o resultado.

EXERCÍCIOS

1) Gerar um sinal f(x)=1+sin(2100t), amostrado a 1KHz.

2) Passar o sinal por um filtro passa-baixas Butterworth com fc=50 Hz.

3) Passar o sinal por um filtro passa-baixas Butterworth com fc=50 Hz.

4) O Matlab possui a função freqz, que permite plotar a resposta em freqüência de filtros digitais.Examine o help desta função, e use-a para plotar a resposta em freqüência do filtro do últimoexemplo.