Processamento Digital de Sinais, Notas de aula de Engenharia Elétrica e Eletrônica. Faculdade de Tecnologia do Estado de São Paulo (FATEC)
Edilson.Camargo
Edilson.Camargo12 de Abril de 2016

Processamento Digital de Sinais, Notas de aula de Engenharia Elétrica e Eletrônica. Faculdade de Tecnologia do Estado de São Paulo (FATEC)

PDF (3.1 MB)
118 páginas
172Número de visitas
Descrição
Livro texto para ser usado em um curso de processamento digital de sinais
20pontos
Pontos de download necessários para baixar
este documento
baixar o documento
Pré-visualização3 páginas / 118
Esta é apenas uma pré-visualização
Consulte e baixe o documento completo
Esta é apenas uma pré-visualização
Consulte e baixe o documento completo
Pré-visualização finalizada
Consulte e baixe o documento completo
Esta é apenas uma pré-visualização
Consulte e baixe o documento completo
Esta é apenas uma pré-visualização
Consulte e baixe o documento completo
Pré-visualização finalizada
Consulte e baixe o documento completo

Processamento Digital de Sinais

Página - 1

Índice

Conceitos Básicos

1 - Introdução

2 - Matlab

3 - Processamento Digital de Sinal

4 - Código Matlab aplicado a Processamento de Sinais

Processamento Digital

5 - Processamento Digital de Imagens

6 - Processamento de Imagens com Matlab

7 - Processamento Digital de Áudio com Matlab

8 - Sistemas Discretos

Analise de Sinais

9 - Análise de Fourier

10 - Transformada Z

11 - Filtros Digitais

12 - Filtros IIR

Projetos de Filtros

13 - Exercícios de Filtros

14 - Projetos de Filtros FIR

15 - Projeto de Filtros por Amostragem de Frequência

16 - Processador Digital de Sinal - DSP

Página - 2

Conceitos Básicos

1. INTRODUÇÃO:

Nesse capítulo são apresentados conceitos básicos relativos ao

processamento digital de sinais. Inicialmente, é discutida a arquitetura de

sistemas de comunicação. Em seguida, o processamento de sinais e

definido, identificando-se seu objeto, seus agentes e suas ações, bem como.

é definida a sua arquitetura genérica. Posteriormente, é realizada uma

classificação de sinais. Por fim, é apresentada uma arquitetura de sistemas

de processamento digital.

1.1 Arquitetura de sistemas de comunicação

a) Arquitetura básica:

– Fonte de sinal.

– Transmissor.

– Canal (ou meio) de transmissão (ou de comunicação).

– Receptor.

– Destino.

b) Adequações em comunicação digital:

– Codificação de fonte.

– Codificação de canal.

1.2 Processamento de sinais

a) Objeto do processamento: sinal é definido como uma entidade que

carrega informação.

b) Agente do processamento: sistema.

 “Um sistema é um conjunto de elementos, que interagem entre si, com

o objetivo de realizar uma determinada função”.

 Arquitetura de um sistema: variáveis, elementos, topologia e função.

c) Domínio do processamento: domínio no qual a função do agente é

definida.

 Tempo/espaço (forma) × frequência (composição espectral).

Página - 3

d) Ação do processamento: função exercida pelo agente sobre o objeto.

 Conformação (tempo/espaço) × Alteração espectral (frequência).

e) Arquitetura genérica do processamento:

 Sinal de entrada (ou estímulo ou excitação).

 Condições iniciais (valores de todas as variáveis internas do sistema).

 Sistema.

 Sinal de saída (ou resposta).

f) Nomenclatura usual: “Sinal” (sinal desejado) × “Ruído” (sinal indesejado).

1.3 Tipos de sinais

a) Sinal: entidade que carrega informação.

b) Visão matemática de sinal: variável funcionalmente dependente de uma

ou mais variáveis independentes. Ex.: w = f(x, y).

c) Visão física de sinal: grandeza física.

d) Tipos de sinais de acordo com o número de variáveis independentes:

unidimensional, bidimensional, tridimensional, multidimensional.

Ex.: áudio = f(t), imagem = f(x, y), vídeo = f(x, y, t) e

tomografia/sismologia = f(v1, v2, · · · , vN, t).

e) Tipos de sinais de acordo com o tipo das variáveis:

 Sinal analógico: todas as variáveis são contínuas.

 Sinal amostrado: discretização das variáveis independentes

(amostragem).

 Sinal quantizado: discretização da variável dependente (quantização).

 Sinal digital: todas as variáveis são discretas (amostragem +

quantização).

f) Sinal amostrado e sinal digital: conjunto ordenado de valores numéricos

(sequencia numérica).

1.4 Arquitetura de sistemas de processamento digital

a) Sinal de entrada analógico: comumente, um sinal elétrico (tensão ou

corrente).

Página - 4

b) Pré-processamento analógico:

 Filtro anti-aliasing: com seletividade em frequência do tipo passa-

baixa.

 Amostragem e retenção (sample-and-hold ou S/H): responsável por

manter fixo o valor da amostra durante o tempo necessário para que

ela seja convertida em um número.

 Conversor Analógico-Digital (A/D).

c) Sinal de entrada digital: representação numérica computacional.

d) Processador de sinal digital (Digital Sinal Processor ou DSP).

e) Sinal de saída digital: representação numérica computacional.

f) Pós-processamento analógico:

 Conversor Digital-Analógico (D/A).

 Filtro de suavização (smoothing): com seletividade em frequência do

tipo passa-baixa.

g) Sinal de saída analógico: comumente, um sinal elétrico (tensão ou

corrente).

1.5 Implementações analógicas e digitais

1.5.1 Definição do problema

Dado que um mesmo sistema pode ser implementado por sinais/componentes

analógicos e digitais, naturalmente surge a questão de qual das opções é a melhor.

Diversas são as comparações encontradas na literatura técnica entre

implementações de sistemas que empregam sinais/componentes analógicos e

digitais. Adeptos de ambos os tipos de implementação facilmente definem

parâmetros que os favorecem e normalmente estabelecem comparações tentando

mostrar que a sua implementação de preferencia é, realmente, a melhor opção.

Antes de tecer meras comparações absolutas, muitas vezes passionalmente

polarizadas, por um tipo ou outro de implementação, deve-se levar lembrar que

estão sendo comparadas duas alternativas de implementação essencialmente

diferentes. Logo, os parâmetros de comparação normalmente empregados

fornecem apenas características individuais de cada tipo de implementação ao invés

de estabelecer bases reais de comparação.

Página - 5

Por um lado, uma implementação analógica mapeia sinais matemáticos em

sinais físicos e realiza as operações matemáticas através de componentes físicos

que possuem equações de definição capazes de realizar os cálculos necessários,

isolada ou conjuntamente com outros componentes.

Pode-se dizer que o problema matemático é transformado em um problema

físico e implementado por sinais/componentes físicos. Por sua vez, uma

implementação digital, apesar de obviamente empregar componentes físicos, não

utiliza os seus sinais e as suas equações diretamente.

Ao invés disso, é realizada uma codificação mais complexa, de tal forma que a

implementação simula diretamente as operações e os operandos matemáticos.

Nesse caso, pode-se dizer que o problema matemático não é mapeado em um

problema físico, sendo apenas simulado em uma implementação física, mas que é

virtualmente matemática.

A comparação direta torna-se ainda mais sem sentido a partir da constatação de

que sistemas digitais podem ser implementados por software, por hardware digital

ou por ambos simultaneamente.

Com tais conceitos em mente, não é difícil perceber que qualquer tentativa para

estabelecer parâmetros de comparação entre implementações analógicas e digitais

conduz apenas a definições de característica relativas de cada uma das opções

consideradas.

Também não é difícil perceber que, devido `as características próprias de

cada uma das implementações, cada uma delas poderá ser proposta como a melhor

solução para problemas específicos, os quais necessitem de tais características.

Por vezes, até mesmo uma solução que misture ambos os tipos de implementação

pode ser a melhor escolha.

Finalmente, é importante ressaltar que o levantamento de parâmetros

comparativos é importante e deve ser efetuado. NA o para que se defina qual das

duas opções de implementação é absolutamente a melhor, mas para que se possa

estabelecer uma base de dados que fundamente a decisão de projeto diante de

cada problema diferente.

Página - 6

1.5.2 Características das implementações analógicas e digitais

A seguir, levando-se em consideração uma implementação eletrônica, são

apresentadas algumas características da implementação de sistemas para os casos

analógico e digital.

a) Sinais;

 Analógico: valores contínuos de tensão e de corrente.

 Digital: sequencias de números, representados por uma determinada

codificação, contendo um número finito de símbolos. Normalmente, é

utilizada uma codificação binária.

b) Componentes básicos;

 Analógico: componentes passivos (resistor, capacitor e indutor) e

componentes ativos (transistor, Amp. Op.).

 Digital: atrasador unitário (registrador), multiplicador e somador.

c) Armazenamento de sinal por longo prazo;

 Analógico: por ser um processo que envolve valores contínuos de

grandezas físicas, sofre grande degradação.

 Digital: uma vez que, geralmente, envolve codificação binária, sofre

pequena degradação.

d) Ocupação de área;

 Analógico: de forma geral, menor ocupação.

 Digital: de forma geral, maior ocupação.

e) Repetitividade/reprodutibilidade;

 Analógico: uma vez que os parâmetros são físicos e contínuos,

necessita de um bom processo de fabricação.

 Digital: dado que os parâmetros são matemáticos, a fabricação é

repetitiva por construção.

f) Variabilidade na fabricação;

 Analógico: os componentes possuem um valor nominal e uma

incerteza associada ao processo de fabricação. Devem ser utilizadas

técnicas de projeto de sistemas que controlem a sensibilidade `a

variação dos valores dos componentes.

 Digital: valor matemático fixo, associado a quantidade símbolos

utilizados na codificação numérica.

Página - 7

g) Variabilidade na operação;

 Analógico: os componentes são influenciados por fatores ambientais

(temperatura, humidade, etc.), podem ser sujeitos a envelhecimento

(aging) e podem sofrer desgaste por uso. Devem ser utilizadas

técnicas de projeto de sistemas que controlem a sensibilidade `a

variação dos valores dos componentes.

 Digital: operação baseada em valores matemáticos fixos,

representados por uma quantidade finita de símbolos utilizados na

codificação numérica.

h) Variabilidade com as frequências envolvidas nos sinais;

 Analógico: as dimensões e a funcionalidade dos componentes podem

ser fortemente afetadas pela faixa de valores de frequência utilizada.

 Digital: na o afetado.

i) Fontes de erro;

 Analógico: além da variação provocada por fatores ambientais, os

componentes são diretamente afetados por ruídos dos mais variados

tipos, provocados pelos mais diversos mecanismos, sendo intrínsecos

aos próprios componentes e/ou induzidos por fontes externas.

 Digital: devido ao número finito de símbolos usado na codificação

numérica, surgem os seguintes problemas de aproximação numérica:

quantização dos valores das sequencias, quantização dos valores dos

coeficientes dos multiplicadores e aproximação dos valores finais das

operações (soma e multiplicação).

j) Programabilidade;

 Analógico: componentes podem ser fixos e/ou variáveis. No caso de

componentes fixos, devem ser desenvolvidas técnicas de projeto de

sistemas que permitam o controle da variação de um ou mais

parâmetros do sistema.

 Digital: naturalmente programável.

k) Complexidade funcional;

 Analógico: implementada com dificuldade.

 Digital: facilmente implementada.

l) Multiplexação temporal de sinais;

 Analógico: difícil implementação.

Página - 8

2. MatLab

O aplicativo Matlab é uma das ferramentas mais úteis para as áreas de

processamento de sinais, controle, e outras aplicações. O “Mat” do nome Matlab

está relacionado com “Matriz”, ou seja, é um laboratório que permite a manipulação

eficiente de matrizes. Pode-se dizer também que o Matlab permite também a

manipulação eficiente de vetores, já que um vetor de tamanho N pode ser visto

como uma matriz de Nx1. Uma das melhores formas de se familiarizar com o Matlab

é através do comando intro, que provê uma turnê pelas características básicas do

programa.

A aplicação predominante no presente curso será o processamento digital se sinais.

Para isto, preparamos algumas técnicas em Matlab que facilitarão o uso do

aplicativo para esta tarefa.

2.1Gerações de Sinais Amostrados

Como 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:

Para permitir o processamento de sinais por computador, o sinal analógico deve ser

Página - 9

digitalizado por um conversor analógico para digital (A/D). Neste processo o sinal é

amostrado em vários instantes, em intervalos de tempo constante. O intervalo de

tempo é o período de amostragem (T), e o recíproco de T é a frequência de

amostragem (fs= 1/T). A unidade de tempo do período de amostragem em geral é

segundos, e a de Frequência é em Hz. 100 Hz significa que 100 amostras do sinal

são tomadas em 1 segundo. Pode-se simular funções amostradas no Matlab. O

exemplo abaixo mostra como se pode simular a amostragem de uma senóide no

Matlab:

Exercício 1

Siga os passos descritos para simular a amostragem de uma onda senoidal no

seguinte exemplo. 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 e vírgula evita que a variável seja exibida na tela. Tente rodar o comando

acima e observe o que acontece. Caso você tenha errado, e o vetor inteiro comece

a aparecer na tela, você pode parar o processo usando CTRL-C.

O comando size(t) exibe o tamanho da matriz t, e o comando length(t) mostrará o

comprimento do vetor.

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 do vetor 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.

Página - 10

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

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 o eixo y mostrou simplesmente o número da amostra, mas não o

tempo correspondente a cada abcissa y.

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 leia

cuidadosamente 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:

Página - 11

plot(t,f,’b*’,t,f2,’c-.’)

E verifique os efeitos.

Agora tente os seguintes comandos:

stem(t,f)

stem(t,f2)

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 cores do traçado quanto a função plot.

2.2 Sinais 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

Matlab sempre interpretará o vetor como uma matriz, e o programa sempre

interpreta uma multiplicação como uma multiplicação de matrizes. Entretanto, há

uma opção que permite ao usuário fazer uma multiplicação ponto a ponto de dois

vetores ou matrizes (isto é o primeiro elemento do resultado é o produto dos

primeiros elementos dos fatores, o segundo elemento do resultado é o produto dos

segundos 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 2

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 em questão:

x=0:0.05:4;

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

plot(x,y)

Página - 12

Note o comando “.*”. Cada expressão entre parênteses é um novo vetor, e os

vetores são multiplicados ponto a ponto, resultando no efeito desejado.

Agora geraremos a função:

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

Para isto, criamos primeiro o vetor representando a variável x, e em seguida a

função em

questã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 da matriz à potência desejada. Tente digitar x^2, e veja que este comando

não funcionará. Tente descobrir 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 de intervalo 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 o intervalo, 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ça que 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? Use as seguintes escolhas de freqüência de amostragem: 5 Hz, 10 Hz, 20

Hz, 40 Hz, e 200 Hz. Use a funçã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 a função intro, e observe o uso da função bar. Em seguida digite help

bar, para analisar aquela função.

Página - 13

2.3 Exercícios de Aplicação

Nesta seção, veremos algumas funções extras do Matlab. Primeiramente,

consideremos o seguinte trecho a seguir. Analise atentamente o programa abaixo

(consulte o help para as funções que você não conhece):

% exemplo: geracão de uma forma de onda ruidosa

t=0:0.01:pi; %linha 1

y=2*sin(5*t); %linha 2

figure(GCF) %linha 3

plot(t,y) %linha 4

pause %linha 5

ruido=randn(1,length(t)); %linha 6

plot(t,ruido) %linha 7

pause %linha 8

y_ruido=y+ruido; %linha 9

plot(t,y_ruido); %linha 10

Execute o programa. A linha 1 cria o vetor do tempo de amostragem. A linha 2 cria

uma função senoidal (qual é a freqüência?). A linha três faz com que a janela com o

gráfico seja colocada à frente de todas as outras. A linha 4 plota uma senóide, como

mostrado abaixo:

O comando pause, na linha 5, causa uma pausa no programa. Para continuar o

programa, o usuário deve pressionar <Enter>. [Pergunta: o que acontece se este

comando for substituído por pause(2)? Dica: use o comando help].

Página - 14

A linha 6 utiliza o comando randn (utilize o help para aprender sobre este

comando). Este comando gera uma seqüência de números aleatórios com

distribuição normal, com uma linha e tantas colunas quanto for o comprimento do

vetor t. O comando lenght(t) fornece o número de elementos de t. Observe a

plotagem do sinal de ruído:

Em seguida, a linha 9 calcula uma nova variável y_noise, que é o sinal senoidal

mais o ruído. Observe o resultado:

Neste exemplo, você aprendeu como gerar um sinal com ruído. Agora, vejamos

outros comandos. Para isto considere o seguinte código:

% Exemplo 2

A=[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];

pause

size(A)

size(B)

Página - 15

pause

length(A)

length(B)

Após executar o código acima, realize manualmente, e em seguida usando o

Matlab, as seguintes 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 1

B=mean(A) %linha 2

C=std(A) %linha 3

Pause %linha 4

RUIDOSO=randn(1,5000); %linha 5

figure(gcf) %linha 6

plot(RUIDOSO) %linha 7

title('sinal ruidoso') %linha 8

pause %linha 9

MEDIA_RUIDOSO=mean(RUIDOSO) %linha 10

DP_RUIDOSO=std(RUIDOSO) %linha 11

Pause %linha 12

MAISRUIDO=5*randn(1,5000); %linha 13

Página - 16

MEDIA=mean(MAISRUIDO) %linha 14

DP=std(MAISRUIDO) %linha 15

Use o comando help para estudar os comandos mean, e std. Execute o programa.

Tome cuidado de só incluir o ponto-e-vírgula quando necessário. Estes comandos

calculam a média e o desvio padrão dos vetores em questão (relembre na literatura

o que estes termos significam). Na linha 2, o resultado, B, será a média dos

elementos do vetor A Na linha 3, o resultado C será o desvio padrão dos elementos

do vetor A. Calcule A e B manualmente, e cheque os resultados com os resultados

do Matlab. Na linha 5 um sinal de ruído normal é gerado, e é plotado na linha 6. A

média e o desvio padrão são calculados e mostrados nas linhas 10 e 11. Note que a

função randn gera ruído aleatório com média 0, e desvio padrão 1 (na verdade,

próximos de 0 e 1). Para obter um ruído com maior intensidade (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. Plote o sinal resultante. Determine os valores máximo e

mínimo deste sinal. Não se esqueça de incluir tí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 onda quadrada de 10 Hz, amostrada a uma freqüência de 1 KHz.

Página - 17

3. Processamento Digital de Sinais

Sinais estão presentes em diversas situações do dia-a-dia do ser humano. Um sinal

pode ser definido como uma função que carrega uma informação. A forma mais

comum para nós é a comunicação por sinal de voz. Nesse exemplo, temos o sinal

gerado pelo trato vocal e o sinal recebido pelo sistema auditivo. Apesar de ser o

mesmo sinal transmitido a forma como ele é processado é inerente ao receptor. O

processamento de sinais lida com a representação, transformação e manipulação

dos sinais e da informação que eles contêm. Até a década de 60, a tecnologia para

processamento de sinais era basicamente analógica. A evolução de computadores

e microprocessadores juntamente com diversos desenvolvimentos teóricos causou

um grande crescimento na tecnologia digital, surgindo o processamento digital de

sinais (PDS). Um aspecto fundamental do processamento digital de sinais é que ele

é baseado no processamento de sequências de amostras. Para tanto, o sinal

contínuo no tempo é convertido nessa sequência de amostras, i.e., convertido em

um sinal discreto no tempo. Após o processamento digital, a sequência de saída

pode ser convertida de volta a um sinal contínuo no tempo.

A maior parte do processamento de sinais envolve processar um sinal para obter

outro sinal. Normalmente, isso é conseguido por um processo conhecido como

filtragem.

Sinais podem ser classificados em quatro diferentes categorias dependendo de

características de tempo e dos tipos de valores que eles podem assumir. Sinais

contínuos no tempo (ou analógicos) são definidos para qualquer valor de tempo e

eles assumem valores no intervalo contínuo (a, b), onde a pode ser -∞ e b pode ser

+∞. Podem ser representados por uma função de variáveis contínuas. Sinais

discretos no tempo são definidos apenas para certos valores específicos de tempo.

Podem ser representados matematicamente por uma sequência de números reais

ou complexos, x. O n-ésimo número dessa sequência é denotado por x[n]. Assim, x

é formalmente escrito como:

onde n é um inteiro. Tais sequências são geradas a partir de um processo de

amostragem periódica de um sinal analógico. Assim, o valor numérico do nésimo

número da sequência é igual ao valor do sinal analógico xa(t) no tempo nT, i.e.:

Página - 18

Os valores de amplitude de sinais contínuos ou discretos no tempo podem ser

contínuos ou discretos. Se um sinal pode assumir qualquer valor dentro de um

espaço finito ou infinito, ele é dito um sinal contínuo em valores. Sinais digitais são

aqueles para os quais tanto o tempo quanto a amplitude são discretos. Ou seja, ele

é discreto no tempo e só pode assumir valores dentro de um conjunto finito de

possíveis valores (é discreto em valores).

Sinais também podem ser classificados em determinísticos ou aleatórios. Qualquer

sinal que podem ser unicamente descritos por uma expressão matemática, uma

tabela de dados ou uma regra bem definida é chamado determinístico. Esse termo é

usado para destacar que quaisquer valores passados, presentes e futuros do sinal

são conhecidos precisamente, sem incerteza. No entanto, em aplicações práticas,

os sinais não podem ser representados precisamente por equações matemáticas ou

suas descrições são muito complexas para uso. Isso indica que tais sinais têm

comportamentos imprevisíveis sendo chamados de sinais aleatórios.

3.1 Principais Tipos de Sinais

Em um estudo sobre processamento digital de sinais, alguns sinais são de mais

importância. Dentre eles, temos o impulso unitário, d[n], definido como:

Um dos mais importantes aspectos do impulso é que uma sequência arbitrária pode

ser representada como uma soma de impulsos escalonados e deslocados. Por

exemplo, a sequência p[n] abaixo:

Página - 19

pode ser representada como:

De forma mais geral, qualquer sequência x[n] pode ser representada como:

Outra sequência importante é o degrau unitário, u[n]:

O degrau relaciona-se com o impulso como:

Uma forma alternativa de representar o degrau em termos de impulso é obtida

interpretando o degrau em termos de uma soma de impulsos deslocados. Isso pode

ser expresso como:

Página - 20

Por outro lado, o impulso relaciona-se com o degrau unitário como:

Uma sequência exponencial é importante na análise de sistemas discretos e

invariantes no tempo. A forma geral de uma sequência exponencial é dada por:

3.2 Sistemas Lineares e Invariantes no Tempo

Uma classe importante de sistemas consiste naqueles que são lineares e

invariantes no tempo. Como dito acima, os sistemas lineares são aqueles que

obedecem ao princípio da superposição. Se a propriedade da linearidade é

combinada com a representação de uma sequência geral como uma combinação de

impulsos, então um sistema linear pode ser completamente caracterizado pela sua

resposta ao impulso. Seja hk[n] a resposta do sistema a δ[n – k]. Assim, como:

então

Página - 21

Pelo princípio da superposição, podemos escrever:

De acordo com essa equação, a resposta do sistema a qualquer entrada pode ser

expressa em termos da resposta a δ[n – k]. A propriedade da invariância no tempo

implica que, se h[n] é a resposta a δ[n], então a resposta a δ[n - k] é h[n – k]. Com

isso, podemos dizer que:

(Eq. 3.1)

Como consequência, um sistema linear invariante no tempo é completamente

descrito por sua resposta ao impulso. Essa equação é conhecida como soma de

convolução (convolution sum) que pode ser representada pela notação:

(Eq. 3.2)

Apesar da semelhança na notação, deve-se salientar que a soma de convolução

para sinais discretos não é uma aproximação da integral de convolução.

Propriedades da soma de convolução:

1) Comutatividade:

Isso pode ser facilmente justificável com uma mudança de variável na Eq. 3.1.

Especificamente, podemos fazer m = n – k.

2) Distributividade:

x[n]*(h1[n] + h2[n]) = x[n]*h1[n] + x[n]*h2[n]

Página - 22

3) Conexão em Cascata

4) Conexão em Paralelo

5) Causalidade

Como definido anteriormente, um sistema é dito causal se sua resposta não

depende de eventos futuros. Ou seja, para calcular a saída de y[n0], precisamos

apenas de x[n], n < n0. Isso implica na condição:

h[n] = 0, n < 0

Assim, para testar a causalidade basta testar se h[n] = 0 para n<0.

Página - 23

6) Estabilidade

A estabilidade é garantida se:

Para qualquer que seja a entrada x[n] de um sistema:

Assim, em geral, se um sistema linear invariante no tempo tem uma resposta ao

impulso h[n], então seu sistema inverso, se existir, tem resposta ao impulso hi[n]

definida pela relação:

Uma classe importante de sistemas lineares invariantes no tempo consiste daqueles

para os quais x[n] e y[n] se relacionam através de uma equação de diferenças de

coeficientes constantes lineares de n-ésima ordem da forma:

(Eq. 3.3)

Um exemplo de um tal sistema é um acumulador definido pela sequência cujo

diagrama de blocos pode ser visto na figura abaixo:

Página - 24

Tal sistema é representado pela equação de diferenças:

y[n] = y[n – 1] + x[n]

ou y[n] - y[n – 1] = x[n]

Pela Eq. 3.3, temos: N = 1, a0 = 1, a1 = -1, M = 0 e b0 = 1.

Assim, para cada valor de n a saída é dada pela entrada x[n] somada com o valor

anterior do acumulador, y[n – 1].

3.3 Operações entre sequencias

Vamos descrever algumas operações básicas em sequências.

a) Adição de sequências:

A adição de amostra por amostra é dada por:

{x1[n]} + {x2[n]} = {x1(n) + x2(n)}

Deve ser observado que o comprimento das sequências x1[n] e x2[n] deve ser o

mesmo. Se as sequências têm comprimentos diferentes, a menor deve ser

completada para que tenha o mesmo comprimento da maior. Normalmente, isso é

feito, acrescentando zeros à sequência (zero padding).

b) Multiplicação de sequências:

Novamente, é uma operação amostra por amostra e as questões de comprimento

das sequências devem ser consideradas:

{x1[n]}.{x2[n]} = {x1(n).x2(n)}

Página - 25

c) Mudança de escala:

Cada amostra de uma sequência é multiplicada por um escalar α:

α.{x[n]} = {α.x(n)}

d) Deslocamento:

Cada amostra x(n) é deslocada k posições:

y[n] = {x(n – k)}

Seja m = n – k, então n = m + k e a operação pode ser vista como:

y[m + k] = {x(m)}

e) Inversão:

A sequência é posta de trás para frente. Seja x[n] uma sequência de comprimento k.

Logo, y[n] será:

y[n] = {x(k – n)}

f) Soma de amostras:

Soma as amostras de uma sequência dentro de um intervalo:

g) Produto de amostras:

Similar ao anterior, mas com operação de produto em um intervalo.

Página - 26

h) Energia:

A energia de uma sequência x[n] é dada por:

i) Potência:

A potência média de uma sequência periódica x(n) pode ser calculada como:

Página - 27

4. Códigos MatLab aplicados a Processamento de sinais

A seguir são apesentados alguns códigos em Matlab que serão úteis em

processamento de sinais digitais.

4.1 Função Impulso

function [x, n] = impseq(n0, n1, n2) % Impulso

n = [n1:n2];

x = [(n-n0) == 0];

stem (x);

Exemplo 1:

>> impseq (5, 0, 10);

Exemplo 2.

x[n] = 2.d[n + 2] - d[n – 4], -5 £ n £ 5

>> n = [-5:5];

>> x = 2*impseq(-2, -5,5) - impseq(4, -5, 5);

>> stem (n, x); title ('Exemplo de Sequencia'); xlabel('n'); ylabel('x[n]');

Página - 28

4.2 Função Degrau

function [x, n] = stepseq(n0, n1, n2) % Degrau

n = [n1:n2];

x = [(n-n0) >= 0];

stem (x);

Exemplo 1.

>> stepseq (5, 0, 10);

Página - 29

Exemplo 2.

x[n] = n[u[n] – u[n – 10]] + 10e-0.3(n – 10)[u[n – 10] – u[n – 20]], 0 £ n £ 20

>> n = 0:20;

>> x1 = n.*(stepseq(0,0,20) - stepseq(10,0,20));

>> x2 = 10*exp(-0.3*(n-10)).*(stepseq(10,0,20) - stepseq(20,0,20));

>> x = x1 + x2;

>> stem(n,x); title('Sequencia de Degraus'); xlabel('n'); ylabel ('x[n]');

4.3 Senóide

function x = sinseq(n1,n2) % Senóide

n = [n1:0.1:n2];

x = 3*cos(0.1*pi*n + pi/3) + 2*sin(0.5*pi*n);

stem (x);

Exemplo:

>> sinseq (0, 10);

Página - 30

4.4 Operações em sequências

a) Adição de sinais

y[n] = x1[n] + x2[n]

function [y,n] = sigadd(x1,n1,x2,n2)

n = min(min(n1),min(n2)):max(max(n1),max(n2));

y1 = zeros(1, length(n));

y2 = y1;

y1(find((n>=min(n1))&(n<=max(n1))==1)) = x1;

y2(find((n>=min(n2))&(n<=max(n2))==1)) = x2;

y = y1 + y2;

b) Multiplicação de sinais

y[n] = x1[n].x2[n]

function [y,n] = sigmult(x1,n1,x2,n2)

n = min(min(n1),min(n2)):max(max(n1),max(n2));

y1 = zeros(1, length(n));

y2 = y1;

y1(find((n>=min(n1))&(n<=max(n1))==1)) = x1;

y2(find((n>=min(n2))&(n<=max(n2))==1)) = x2;

y = y1.*y2;

c) Deslocamento

y[n] = x[n – k]

Página - 31

function [y,n] = sigshift(x, m, n0)

n = m + n0;

y = x;

d) Inversão

y[n] = x[-n]

function [y,n] = sigfold(x,n)

y = fliplr(x);

n = -fliplr(n);

Exemplo: Seja x[n] = {1, 2, 3, 4, 5, 6, 7, 6, 5, 4, 3, 2, 1}. O valor em negrito

corresponde ao centro da sequência.

Sobre as sequências, temos que:

>> n = -2:10;

>> x = [1:7, 6:-1:1];

 Plote x1[n] = 2x(n – 5) – 3x[n + 4].

>> [x11, n11] = sigshift(x, n, 5);

>> [x12, n12] = sigshift(x, n, -4);

>> [x1, n1] = sigadd(2*x11,n11,-3*x12, n12);

>> stem (n1, x1); title(‘Sequencia’); xlabel (‘n’); ylabel (‘x1(n)’);

Página - 32

 Plote x2[n] = x[3 – n] + x[n].x[n – 2]

>> [x21, n21] = sigfold(x, n);

>> [x21, n21] = sigshift(x21, n21,3);

>> [x22, n22] = sigshift(x, n,2);

>> [x22, n22] = sigmult(x, n, x22, n22);

>> [x2, n2] = sigadd(x21, n21, x22, n22);

>> stem (n2, x2); title('Sequencia');

>> xlabel ('n'); ylabel ('x2(n)');

e) Convolução

Considere as sequências:

x = [3, 11, 7, 0, -1, 4, 2], -3 £ n £ 3

h = [2, 3, 0, -5, 2, 1]; -1 £ n £ 4

onde, novamente, os termos em negrito indicam a origem do eixo das abscissas.

As sequências podem ser vistas abaixo:

Página - 33

Podemos usar a função conv do MatLab diretamente:

>> x = [3, 11, 7, 0, -1, 4, 2];

>> h = [2, 3, 0, -5, 2, 1];

>> y = conv (x, h);

y = 6 31 47 6 -51 -5 41 18 -22 -3 8 2

O problema do uso da função conv é que não sabemos, na resposta, onde está

a origem da sequência. Para tanto, vamos criar uma nova função:

function [y, ny] = conv_m (x, nx, h, nh)

nyb = nx(1) + nh(1);

nye = nx(length(x)) + nh(length(h));

ny = [nyb:nye];

y = conv(h, x);

>> x = [3, 11, 7, 0, -1, 4, 2];

>> nx = [-3:3];

>> h = [2, 3, 0, -5, 2, 1];

>> nh = [-1:4];

>> [y, ny] = conv_m (x, nx, h, nh)

y = 6 31 47 6 -51 -5 41 18 -22 -3 8 2

ny = -4 -3 -2 -1 0 1 2 3 4 5 6 7

A amplitude -51 está no ponto de origem (ny = 0).

Página - 34

4.5 Equações de Diferenças e Resposta ao Impulso

Exemplo: Dada a seguinte equação de diferenças:

y[n] – y[n – 1] + 0.9y[n – 2] = x[n], para todo n

a) Calcule e plote sua resposta ao impulso h[n] para n = -20,.., 120.

Como vimos anteriormente, uma equação de diferenças é da forma:

De acordo com a equação dada, temos:

a = [1, -1, 0.9] e b = [1]

No MatLab, fazemos:

>> x = impseq(0, -20, 120);

>> n = [-20:120];

>> h = filter(b, a, x);

>> stem(n, h); title('Resposta ao impulso'); xlabel('n'); ylabel('h[n]');

Página - 35

b) Calcule e plote sua resposta ao degrau s[n] para n = -20,.., 120.

No MatLab, fazemos:

>> x = stepseq(0, -20, 120);

>> n = [-20:120];

>> h = filter(b, a, x);

>> stem(n, h); title('Resposta ao degrau'); xlabel('n'); ylabel('s[n]');

Página - 36

c) O sistema é estável? Como vimos, um sistema é estável se:

Assim, no MatLab, basta fazermos:

>> sum(abs(h))

Ans = 14.8785

Logo, o sistema é estável.

Página - 37

5. Processamento digital de imagens

Ao longo dos anos, o interesse em técnicas de processamento digital de

imagens (PDI) vem aumentando. Diversas são as aplicações que precisam utilizar

algum algoritmo associado ao processamento digital de imagens. Na prática, os

principais objetivos relacionados com PDI são: i) melhoria na qualidade da imagem

para observação humana ou para um mais reconhecimento automático por

máquinas; e ii) compressão para transmissão ou armazenamento de imagens. Cada

um desses objetivos leva a um grande leque de técnicas e algoritmos.

No primeiro caso, podemos destacar aplicações relacionadas com imagens via

satélite, onde, em geral, as imagens chegam com distorções e ruídos que precisam

ser filtrados para uma análise humana. Da mesma forma, a quantidade de estudos

voltados para a visão computacional vêm crescendo mundialmente e, na maioria

dos casos, o primeiro passo para um sistema de visão computacional é a aplicação

de algoritmos de processamento de imagens para, por exemplo, extrair a borda de

objetos, permitindo uma melhor análise da sua forma.

No segundo caso, em termos de compressão e transmissão, podemos observar

aplicações voltadas para a Internet que, cada vez mais, precisam usar imagens para

transmitir informações. Algumas dessas aplicações exigem imagens de alta

resolução, como o sensoriamento remoto. Nesses casos, é preciso usar algoritmos

que provoquem uma grande redução na quantidade de dados, sem perder

informação. Esse é o caso também de imagens médicas, mas, nesse tipo de

aplicação, a preocupação maior é o armazenamento das imagens. Por exemplo,

nos Estados Unidos, mamografia digitais precisam ser armazenadas por 5 anos

após o exame. Cada exame gera diversas lâminas que são imagens de alta

resolução ocupando GigaBytes de espaço de armazenamento. Nesse caso

também, é essencial que a compressão seja feita sem perdas.

Uma imagem digital pode ser entendida como uma função f(x,y), onde x e y são as

coordenadas espaciais e a função f é o nível de cinza (ou de brilho) naquele ponto.

Quando x, y e f estão numa escala finita e discreta, dizemos que temos uma

imagem digital.

A imagem digital pode ser gerada diretamente através de um dispositivo digital ou

pode ser digitalizada a partir de uma imagem real. Esse segundo caso é o mais

clássico, onde dispositivos, como scanners, por exemplo, faziam a transposição da

Página - 38

imagem do meio real para o digital. Para tanto, alguns processos são fundamentais

como veremos posteriormente.

Em termos de pesquisa e desenvolvimento, processamento digital de imagens está

inserido na grande área de processamento gráfico. Nessa área, encontramos ainda

elementos de computação gráfica e de visão computacional. [Gomes e Velho, 1995]

relacionam essas áreas de acordo com a representação gráfica da Fig. 9.1. Essa

representação e essa relação estão associadas ao tipo de dado que entra e que sai

de um sistema específico. Notadamente, essa representação está relacionada com

a percepção humana para os dados inseridos e retornados pelos sistemas. Para os

dispositivos, como computadores, tudo é apenas processamento de dados.

Página - 39

6. Processamento de Imagem com Matlab

Nesta atividade exploraremos como imagens são mostradas e representadas como

matrizes no Matlab. A informação de cores é codificada em uma colormap (tabela

de cores a ser usada). Imagens podem ser armazenadas em arquivos

imagen_name.mat e carregadas usando o comando load image_name.

O Matlab tem uma série de imagens padrões disponível.

1. Digite load clown e verifique as variáveis carregadas digitando whos

>> whos

Name Size Bytes Class

X 200x320 512000 double array

caption 2x1 4 char array

map 81x3 1944 double array

Foram criadas três variáveis: caption (guarda as informações de cabeçalho da

imagem), map

(informação de cores), x (guarda a informações de intensidade da imagem).

2. Olhe alguns elementos do vetor x (por exemplo, digite x(65:75, 100:110)). Os

valores na matriz devem ser inteiros. Quantos bits estão sendo usados na

codificação?

3. Para mostrar a imagem, digite image(X). Você deve ver o palhaço numa janela,

apesar de sua cor parecer não-natural.

Página - 40

4. Para obter as cores corretas, digite colormap(map). Agora o palhaço deve

aparecer nas cores corretas.

5. Para aumentar o brilho da imagem, pode-se utilizar o comando brighten(beta)

sendo beta um número entre -1 e 1. Números maiores que zero resultam numa

imagem mais brilhante e menores do que zero numa imagem mais escura. Por

exemplo,

>> brighten(.7)

6. Podem-se utilizar operações matriciais normalmente para trabalhar com imagens.

Por exemplo, transpor a matriz implica na transposição da imagem. Para colocar a

imagem de ponta cabeça, basta inverter a ordem das linhas da matriz.

X1 = X';

X2 = X(end:-1:1,:);

figure(1);

subplot(221); image(X1);

subplot(222); image(X2);

colormap(map);image(X(end:-1:1, :))

colormap(map);

Página - 41

7. Gere a imagem a seguir:

8. Podemos selecionar um pedaço da imagem, pegando algumas linhas e colunas

da matriz. Por exemplo, o olho esquerdo do palhaço pode ser obtido usando:

Xolho = X(50:100, 150:250);

image(Xolho);

colormap(map);

9. Usando o comando filter e a Equação 1 podemos borrar a imagem na horizontal e

na vertical usando os seguintes comandos:

N = 25;

Yvert = filter(ones(1,N)/N,1, X);

subplot(221);image(Yvert);

Página - 42

colormap(map);

Yhori = filter(ones(1,N)/N,1, X');

Yhori = Yhori';

subplot(222);image(Yhori);colormap(map);

10. A imagem pode ser aproximadamente recuperada utilizando um filtro inverso.

Por exemplo,

para desfazer o borrão vertical, usamos:

Yvolta = filter(1, ones(1,N)/N, Yvert);

subplot(223);image(Yvolta);colormap(map);

B. Imagem no fomato jpeg

O formato jpeg é um dos mais utilizados para codificação de imagens. O comando

[X] = imread( FILENAME,’jpeg’) pode ser utilizado para importar uma imagem neste

formato para o Matlab. É gerado um vetor X tridimensional com os componentes

RGB da imagem. A matriz X(:,:,1) contém a informação do vermelho (“red”), a matriz

X(:,:,2) a informação do verde e a matriz X(:,:,3) a informação do azul (“blue”).

O comando IMWRITE(A,FILENAME,’jpeg’) grava a matriz A no arquivo FILENAME.

jpg .

Página - 43

Muitos outros formatos são possíveis. Leia o help dos comandos acima para mais

informações. A seguinte seqüência de comandos lê a figura moinho.jpg e mostra

suas componentes:

X = imread('moinho.jpg','jpeg');

Xred = X(:,:,1);

Xgreen = X(:,:,2);

Xblue = X(:,:,3);

subplot(221);image(X); title('moinho.jpg');

subplot(222);image(Xred); title('Componente red');

subplot(223);image(Xgreen); title('Componente green');

subplot(224);image(Xblue); title('Componente blue');

Num arquivo jpeg, as cores podem ser trabalhadas separadamente. Por exemplo,

os seguintes comandos aumentam a intensidade do verde da imagem. Tente outras

configurações de cores.

X = imread('moinho.jpg','jpeg');

Xred = X(:,:,1);

Xgreen = X(:,:,2);

Xblue = X(:,:,3);

Xgreennovo = double(Xgreen)*3;

Xrednovo = double(Xred);

Xbluenovo = double(X(:,:,3));

Xnovo(:,:,1) = uint8(Xrednovo);

Xnovo(:,:,2) = uint8(Xgreennovo);

Página - 44

Xnovo(:,:,3) = uint8(Xbluenovo);

subplot(221);image(X); title('moinho.jpg');

subplot(222);image(Xnovo); title('Cores modificadas');

C. O que está escrito?

A seguinte foto foi tirada de uma placa com o carro em movimento. Ela esta

armazenada no arquivo: testeplacaverm.mat

Página - 45

7. Processamento de Áudio com MatLab

Um som pode ser gravado com a função wavrecord que gera um arquivo do tipo

wav:

>> som = wavrecord (16000, 8000, 1, ‘double’);

Esse comando grava 16000 amostras com uma taxa de amostragem de 8 kHz e

o armazena no vetor som do tipo double. O parâmetro ‘1’ indica que a gravação

é feita em apenas um canal (mono) e não em estéreo. Para tocar o som, basta

usar o comando soundsc:

>> soundsc (som);

Às vezes, é importante informar a frequência de amostragem:

>> soundsc (som, 8000);

O arquivo de som pode ser plotado em um gráfico como um vetor comum:

>> plot (som);

É comum precisarmos trabalhar com o arquivo de som normalizado. Para tanto,

usamos:

>> som = som/max(abs(som));

Isso não muda o sinal apenas o apresenta com amplitude entre -1 e 1.

A leitura de um arquivo wav pode ser feita com o comando wavread:

>> som = wavread (‘teste.wav’);

ou

>> [som, fs] = wavread (‘teste.wav’);

Nesse segundo caso, a frequência de amostragem do sinal é salva na variável

fs. O comando wavread não permite que o parâmetro de entrada seja uma

variável com o nome do arquivo. A entrada precisa ser especificamente o nome

do arquivo de som.

O processamento do sinal de áudio pode ser feito através de filtros como o

processamento de sinais, usando a função filter:

>> y = filter (b, a, som);

onde a e b são os coeficientes do filtro a ser aplicado (para filtros FIR, a = 1).

Exemplo 1:

>> som = wavread(‘a_casa.wav’);

Página - 46

>> plot (som);

Calculando a transformada de Fourier:

>> som_spec = fft (som, 256);

>> plot (abs(som_spec));

Um gráfico em formato mais padrão pode ser obtido com as baixas frequências

centralizadas:

>> plot (abs(fftshift(som_spec)));

Suponha um filtro IIR com função de transferência: H(z) = 1 – 0,9375z-1:

>> h = [1 -0.9375];

>> y = filter(h, 1, som);

>> soundsc(y, 22000);

Página - 47

Ouviremos um som mais nasal. Vamos fazer o mesmo processamento, mas

dividindo o som em amostras de 240 frames cada que serão filtradas

individualmente e re-agrupadas depois:

>> w = 240;

>> n = floor(length(som)/w);

>> for k=1:n

seg = som(1+(k - 1)*w:k*w);

segf = filter(h, 1, seg);

outsp(1+(k-1)*w:k*w) = segf;

end

>> soundsc(outsp, 22000);

Observem, a seguir (Fig. 10.13), a plotagem do sinal original (som) e filtrado

(outsp), assim como suas transformadas de Fourier:

Fig. 10.13. (coluna da esquerda) Sinal original e sua transformada e (coluna da

direita) sinal filtrado e sua transformada.

Página - 48

Exemplo 2: Como um exemplo da dificuldade de analisar um sinal que muda

constantemente, vamos construir um sinal “ágil na frequência” (um que muda

rapidamente suas características de frequência):

>> y = chirp([0:0.001:5],0,5,500);

>> soundsc (y); % Escute o som para entende-lo!!

>> z = [y, y(length(y):-1:1), y];

>> f = abs(fft(z, 8192));

>> plot(f(1:4096));

Mas esse gráfico representa mesmo o sinal que criamos? Vamos observar

melhor o sinal quebrando ele em janelas e plotando-as como uma “queda

d’água” (waterfall) para ver como as frequências mudam pelo tempo.

>> s = spectrogram(z, 1024);

>> waterfall(abs(s)');

Página - 49

A plotagem em “queda d’água” mostra cerca de 30 pedaços no tempo, cada um

correspondendo a uma FFT-512 e indicando claramente que diferentes

componentes de frequência estão presentes durante cada período de tempo.

Veja o resultado para diferentes FFTs de N-Pontos:

>> s = spectrogram(z, 256);

>> waterfall(abs(s)');

>> s = spectrogram(z, 8192);

>> figure, waterfall(abs(s)');

Página - 50

Exemplo 3: Criando música no MatLab:

Primeiro, vamos crier uma onda senoidal de amplitude A = 1, com uma

frequência de 523,25 Hz (correspondente a um pitch C em um piano; uma oitava

acima do C médio):

cnote = sin(2*pi*523.25*(0:0.000125:0.5));

Esse vetor cnote contém amostras da onda senoidal de t = 0s a t = 0.5s, as

amostras são separadas de 0.000125s (que é o intervalo de amostragem Ts).

Note que esse intervalo de amostragem corresponde à frequência de

amostragem de 8 kHz (1/Ts = fs).

Podemos graver esse som com o comando wavwrite:

wavwrite(cnote, ‘c.wav’);

E temos a primeira nota.

O seguinte site apresenta a frequência de diferentes notas (apresentadas na

tabela a seguir): http://www.dolmetsch.com/musictheory27.htm.

Página - 51

Usando essa informação, podemos criar diferentes notas no MatLab. Observe

que existem diferentes oitavas da mesma forma que existem diferentes teclas

em um piano. Aqui estão algumas no MatLab:

f = sin(2*pi*174.61*(0:0.000125:0.5));

g = sin(2*pi*195.99*(0:0.000125:0.5));

a = sin(2*pi*220*(0:0.000125:0.5));

b = sin(2*pi*246.94*(0:0.000125:0.5));

Pas de música faça:

line1 = [a,b,c,d,e,f];

line2 = [a,b,c,d,e,f];

As letras correspodem às notas que você cria no MatLab de acordo com a

tabela anterior. Coloque as notas na ordem que você quiser tocá-las. Para crier

uma música use:

song = [line1, line2];

e toque com sound() ou soundsc().

Página - 52

8. Sistemas Discretos

Para um sistema linear e invariante (SLI) caracterizado por uma resposta

impulsional h[n], a saída y[n], para uma entrada x[n], é dada pelo somatório da

convolução:

Se o sistema também for causal, então:

A expressão acima pode ser facilmente computada usando-se a função conv(h,x).

Exemplo 8.1: Um SLIC possui uma resposta impulsional dada por h[n] = (0,85)n

u[n]. Para uma entrada x[n] = u[n] - u[n-15], determine as 40 primeiras amostras da

saída y[n].

/Programa Exemplo 8.1/

% Programa para computar a convolução de duas sequências

% h[n]=(0.85)^n e x[n]=u[n] - u[n-15]

N=40;

n=0:N-1;

h=(.85).^n;

x=[ones(1,15),zeros(1,25)];

y=conv(x,h);

subplot(3,1,1)

stem(n,x),xlabel('n'),ylabel('x[n]'),title('Entrada')

subplot(3,1,2)

stem(n,h)

xlabel('n'),ylabel('h[n]'),title('Resposta Impulsional')

subplot(3,1,3)

stem(n,y(1:N)),xlabel('n'),ylabel('y[n]'),title('Saída')

Página - 53

A entrada x[n], a resposta impulsional do sistema h[n], e a saída do sistema y[n] são

mostradas na figura abaixo.

Figura 8.1: Resposta de um SLI a uma entrada tipo pulso.

A saída de um SLI discreto também pode ser obtida através da solução da

equação diferença que descreve este sistema.

O Matlab possui uma função filter(a,b,x), em que a = [1, a1,a2,...aN], b =

[b0,b1,....bM], e x é o vetor com as amostras da entrada, para resolver

numericamente a equação diferença acima.

Exemplo 9.2: Determinar as respostas impulsional e ao degrau do sistema descrito

pela equação diferença y[n] - 1,2728y[n-1] + 0,81y[n-2] = 0,5372x[n].

Página - 54

/Programa Exemplo 8.2/

% Programa para computar as respostas impulsional e ao degrau

%do sistema descrito por y[n]-1.2728y[n-1]+0.81y[n-2] =0,532x[n].

% número de amostras

N = 50;

% coeficientes da ED

b = 0.5372; a = [1 -1.2728 0.81];

% gera o impulso

imp = [1 zeros(1,N-1)];

% gera o degrau unitário

deg = [ones(1,N)];

% cômputo da resposta impulsional

h = filter(b,a,imp);

% cômputo da resposta ao degrau

y = filter(b,a,deg);

n = 0:N-1;

subplot(2,1,1)

stem(n,h)

xlabel('n'),ylabel('h[n]'),title('Resposta Impulsional')

subplot(2,1,2)

stem(n,y)

xlabel('n'),ylabel('y[n]'),title('Resposta ao Degrau')

A Figura 8.2 mostra as respostas impulsional e ao degrau do sistema do exemplo

acima

Página - 55

Figura 8.2: Respostas impulsional e ao degrau do sistema do Exemplo 8.4.

No Exemplo 8.2, o sistema é dito recursivo ou IIR (Infinite Impulse Response),

pois a resposta impulsional é de duração infinita. Há também os sistemas não-

recursivos ou FIR (Finite Impulse Response), cuja resposta impulsional é de

duração finita. Um exemplo destes sistemas é o filtro de média móvel descrito pela

seguinte equação diferença:

Exemplo 8.3: Um sinal s[n] = 3 + 4(0,95)nsen(π/8 n) é contaminado por um ruído

com distribuição uniforme entre -0,5 e 0,5. Será utilizado um filtro de média móvel

para diminuir os efeitos do ruído aditivo.

Página - 56

Figura 8.3: Sinal s[n] e ruído r[n] com distribuição uniforme

As sequências s[n] e r[n] são somadas, resultando na seqüência x[n], que nada

mais é do que o sinal original s[n] contaminado pelo ruído aditivo r[n]. Para se

atenuar o efeito indesejável deste ruído, pode-se passar o sinal x[n] por um filtro de

médias móveis.

No presente exemplo, usou-se um filtro com 3 atrasos (M = 3), a Figura 9.4a

mostra o sinal original s[n], o ruído r[n], e a soma dos dois x[n]. A Figura 9.4b mostra

as curvas da entrada x[n] e da saída y[n] do filtro de médias móveis. Deve-se

observar que a saída y[n] é muito próxima do sinal original s[n], exceto por um

atraso de uma amostra, o que é consequência do processo de filtragem.

Figura 8.4: a) sinal s[n], ruído r[n] e x[n]=s[n]+r[n]; b) sinal s[n] e saída do filtro y[n]

Página - 57

/Programa Exemplo 8.3/

% Programa que exemplifica o uso de filtros de média móvel

N = 50;

n = 0:N-1;

s = (4*sin(n*pi/8).*(.95).^n) + 3*ones(1,N);

% ruído com distribuição uniforme

r = rand(1,N) - 0.5;

figure(1)

subplot(2,1,1)

stem(n,s),xlabel('n'),ylabel('s[n]'),title('Sinal sem Ruido')

subplot(2,1,2)

stem(n,r),xlabel('n'),ylabel('r[n]'),title('Ruido')

x = s + r;

M = 3; % número de atrasos do filtro

b = ones(1,M)/M;

y = filter(b,1,x);

figure(2)

subplot(2,1,1)

plot(n,r,'g-',n,s,'y--',n,x,'b:')

xlabel('n'),ylabel('Amplitudes')

legend('g-','r[n]','y--','s[n]','b:','x[n]')

axis([0 50 -2 8]);

subplot(2,1,2)

plot(n,s,'y-',n,y,'b:'),xlabel('n'),ylabel('Amplitudes')

legend('y-','s[n]','b:','y[n]')

axis([0 50 -2 8]);

Página - 58

9. Analise de Fourier

A transformada de Fourier de uma seqüência discreta é definida como:

em que X(ejw) é uma função contínua e complexa.

Exemplo 9.1: Calcular a transformada de Fourier da seqüência x[n] = (0,5)n u[n].

Observe que na solução acima, o somatório nada mais é do que a soma de uma

progressão geométrica cuja razão é (0,5e-jw).

No Exemplo 9.1, como x[n] é uma seqüência de duração infinita não se pode usar o

Matlab para computar X(ejw) diretamente. Entretanto, pode-se computar X(ejw)

usando-se a expressão acima, no intervalo [0,π] e então traçar gráficos de

magnitude e fase, ou das partes real e imaginária.

Exemplo 9.2: Computar e plotar os espectros de magnitude e fase de X(ejw), bem

como as suas partes real e imaginária.

/Programa Exemplo 9.2/

% Programa para computar e plotar os espectros de amplitude e

% fase da DTFT de x[n]=(0,5)^n u[n], a partir da expressão

% X(w)= exp(jw)/[exp(jw) - 0,5].

N = 256;

% eixo das frequências dividido em N pontos

w = (0:N-1)*pi/N;

ex = exp(j*w);

X = ex./(ex - .5*ones(1,N));

% cômputo do módulo de X(jw)

ampl = abs(X);

Página - 59

% cômputo da fase de X(jw)

fase = angle(X);

rex = real(X);

imx = imag(X);

% normalização do eixo das frequências

wnorm = w/pi;

subplot(2,2,1)

plot(wnorm,ampl)

xlabel('Frequência normalizada') ,ylabel('Magnitude'), title('Espectro de Amplitude')

subplot(2,2,3)

plot(wnorm,fase),xlabel('frequência Normalizada') ylabel('Radianos'),title('Espectro

de Fase')

subplot(2,2,2)

plot(wnorm,rex),xlabel('Freqüência Normalizada'),

ylabel('Real[X(jw)]'),title('Parte Real')

subplot(2,2,4)

plot(wnorm,imx),xlabel('Freqüência Normalizada') ylabel('Imag[X(jw)]'),title('Parte

Imaginária')

Figura 9.1: Gráficos da transformada de Fourier: magnitude e fase; partes real e imaginária.

Página - 60

Nos gráficos da Figura 9.1 deve-se observar que o eixo das frequências está

normalizado, ou seja, a frequência π corresponde ao valor 1. Muitos autores

chamam este ponto frequência de Nyquist, pois é a frequência que corresponde à

metade da frequência de amostragem.

Um outro ponto que deve ser observado é o que se refere à forma de se

plotar uma função complexa X(ejw). Isto pode ser feito através de gráficos que

mostrem as partes real e imaginária, ou de gráficos de magnitude e fase. Esta

última forma é a preferida quando a função complexa estiver relacionada a sistemas

lineares discretos, como por exemplo filtros digitais.

No caso em que a sequência discreta x[n] for de duração finita, então será

possível se utilizar o Matlab para o cômputo da transformada de Fourier. Na

realidade o que se faz é computar a Transformada Discreta de Fourier (DFT), que

é uma sequência discreta, definida como:

Pode-se mostrar que os X[k] são amostras de X(ejw) igualmente espaçadas de 2π/N

no círculo unitário, ou seja X(ejw) é a envoltória das amostras representadas por

X[k].

Exemplo 9.3:

Repetir o Exemplo 9.1, utilizando a DFT para computar X(ejw). Como x[n]

deve ser finita, o que se pode fazer é obter uma versão truncada de x[n] num

intervalo finito adequado.

No programa abaixo foi utilizada a função fft(x,M) para computar a DFT da

sequência x de comprimento N, com M pontos. Se M for omitido o cômputo da DFT

se dará com o número de pontos de x, neste caso N. Se M>N, então serão

acrescentados (M-N) zeros ao final da sequência x. Este procedimento não altera a

forma da envoltória de X[k], apenas a define melhor.

A FFT (Fast Fourier Transform) é apenas um algoritmo que computa de

forma eficiente a DFT e esses algoritmos são mais rápidos quando M for uma

potência inteira de 2 (M= 2p).

Página - 61

/Programa Exemplo 9.3/

% Programa para computar a transformada de Fourier usando o

% algoritmo de FFT

N = 256;

% geração da sequência x[n]

x = (.5).^(0:N-1);

% cômputo da DFT de x[n]

Xc = fft(x);

% frequências no intervalo [0,pi]

X = Xc(1:N/2+1);

% cômputo do módulo de X(jw)

ampl = abs(X);

% cômputo da fase de X(jw)

fase = angle(X);

rex = real(X);

imx = imag(X);

% normalização do eixo das frequências

wnorm = 0:2/N:1;

subplot(2,2,1)

plot(wnorm,ampl)

xlabel('Freqüência Normalizada'),ylabel('Magnitude')

title('Espectro de Amplitude')

subplot(2,2,3)

plot(wnorm,fase)

xlabel('Freqüência Normalizada'),ylabel('Radianos')

title('Espectro de Fase')

subplot(2,2,2)

plot(wnorm,rex)

xlabel('Freqüência Normalizada'),ylabel('Real[X(jw)]')

title('Parte Real')

subplot(2,2,4)

plot(wnorm,imx)

Página - 62

xlabel('Freqüência Normalizada'),ylabel('Imag[X(jw)]')

title('Parte Imaginária')

Comparando-se os resultados mostrados nas Figuras 9.1 com os da Figura 9.2,

pode-se ver que a aproximação obtida utilizando a DFT é muito boa.

Figura 9.2: Gráficos da transformada de Fourier obtidos através da DFT.

Exemplo 9.4:

Computar e plotar a DFT da sequência x[n] = u[n] -u[n-8], com 8, 16, 32, e 64

pontos.

/Programa Exemplo 9.4/

% Programa que ilustra o efeito de se apendar zeros a uma

% seqüência finita

% x[n]=u[n]-u[n-8]

x = ones(1,8);

% fft com 8 pontos

XO = fft(x);

% fft com 16 pontos, 8 zeros acrescentados a x[n]

X1 = fft(x,16);

% fft com 32 pontos, 24 zeros acrescentados a x[n]

X2 = fft(x,32);

Página - 63

% fft com 64 pontos, 56 zeros acrescentados a x[n]

X3 = fft(x,64);

subplot(2,2,1)

stem([(0:4)/4],abs(XO(1:5)))

xlabel('Freqüência Normalizada'),ylabel('Magnitude')

title('N = 8')

subplot(2,2,3)

stem([(0:8)/8],abs(X1(1:9)))

xlabel('Freqüência Normalizada'),ylabel('Magnitude')

title('N = 16 (8 zeros)')

subplot(2,2,2)

stem([(0:16)/16],abs(X2(1:17)))

xlabel('Freqüência Normalizada'),ylabel('Magnitude')

title('N = 32 (24 zeros)')

subplot(2,2,4)

stem([(0:32)/32],abs(X3(1:33)))

xlabel('Freqüência Normalizada'),ylabel('Magnitude')

title('N = 64 (56 zeros)')

O programa acima computa 4 DFTs com diferentes números de pontos, a

primeira com 8 pontos, a segunda com 16 pontos, sendo que destes 8 são zeros

que foram acrescentados ao final de x[n]. As outras duas são computadas com 32 e

64 pontos, sendo que nestes casos foram acrescentados 24 e 56 zeros

respectivamente.

Os resultados mostrados na Figura 8.3 indicam que ao se acrescentar zeros

a uma sequência finita, obtém-se uma amostragem mais fina da envoltória X(ejw),

sem que isto altere a sua forma.

Página - 64

Figura 9.3: O efeito de se acrescentar zeros ao final de uma sequência finita

Exercício 9.5: Para a sequência x[n] = u[n] - 2u[n-8] + u[n-16] compute e plote |X[k]|,

para N=16, e N=64 (acrescentando-se 48 zeros). Comente os resultados.

Página - 65

10. Transformada Z

A transformada z unilateral é definida como:

No caso de SLI a transformada z de h[n] é H(z), que é uma função racional

da forma H(z)=N(z)/D(z), em que N(z) e D(z) são polinômios em z. As raízes de N(z)

= 0 são os zeros de H(z), e as raízes de D(z) = 0 correspondem aos pólos de H(z).

H(z) é chamada de função sistema ou função de transferência, servindo

juntamente com h[n], para caracterizar os SLIs. O Matlab possui funções que

permitem decompor H(z) em seus pólos e zeros ([z,p,k]=tf2zp(num,den)), ou fazer

um gráfico dos pólos e zeros no plano z ( zplane(num,den) ).

Há também uma função que permite obter H(z) na forma racional, a partir de

seus pólos e zeros ( [num,den] = zp2tf(z,p,k) ).

Exemplo 10.1:Expressar a FT, dada abaixo, na forma fatorada, fazer um gráfico de

pólos e zeros e indicar a região de convergência.

/Programa Exemplo 10.1/

% Programa para determinar na forma fatorada uma função

% racional e fazer um diagrama de pólos e zeros

num = [1 -1.2071 -2.6464 2.1213];

den = [1 0.0272 -0.4446 0.5439 0.3240];

[zeros polos ganho] = tf2zp(num,den)

zplane(num,den)

/Resultados/

Página - 66

zeros =

2.0000

-1.5000

0.7071

polos =

0.6364 + 0.6364i

0.6364 - 0.6364i

-0.8000

-0.5000

ganho =

1

A figura 10.1 mostra o diagrama de pólos e zeros de H(z). A região de

convergência é o exterior do círculo de raio | 0,6364 + 0,6364i | = 0,9, que

corresponde ao módulo dos pólos mais afastados da origem, para uma solução

causal.

A transformada inversa de H(z) é a resposta impulsional h[n]. Esta

transformada pode ser computada de diversas maneiras, porém no caso em que

H(z) for racional, é mais conveniente utilizar o método da expansão em frações

parciais. Com este método é possível se obter uma expressão para h[n] na forma

fechada. O Matlab possui a função [r,p,k]=residuez(num,den), que permite computar

os pólos de H(z), os coeficientes da expansão em frações parciais.

Página - 67

Figura 10.1: Diagrama de pólos e zeros de H(z) do exemplo acima.

Exemplo 10.2: Determinar h[n] da FT dada abaixo, pelo método da expansão em

frações parciais.

/Programa Exemplo 10.2/

% Programa para determinar os pólos e os coeficientes da

% expansão em frações parciais de uma função racional em z.

num = [0 1 -1];

den = [1 -1.1 .14 .08];

[coeficientes polos constante] = residuez(num,den)

Página - 68

Resultados

coeficientes =

-0.6667

2.3810

-1.7143

polos =

0.8000

0.5000

-0.2000

constante =

[]

A expansão em frações parciais será da forma:

O que resulta numa transformada inversa causal dada por:

Página - 69

11. Filtros Digitais

Um dos maiores campos de aplicação em processamento digital de sinal é o

dimensionamento de filtros. Em geral, estamos interessados em manipular o sinal.

Por exemplo, podemos querer retirar algum ruído de um sinal, como no caso de um

sinal de voz, onde o ruído deve ser diferenciado da voz propriamente dita. Para isso,

filtros são utilizados.

Filtros estão envolvidos em diversas partes de um sistema de processamento

digital de sinal. Eles podem ser implementados tanto em hardware quanto em

software e atuam em sinais digitais de diversas naturezas, como sons, voz, imagem

ou vídeo. Em cada caso, os filtros assumem particularidades diferentes. Vamos

entender um pouco como se dá o processo em sinais e, em seguida, particularizar

para o caso de imagens digitais.

Filtros digitais são formados por poucos componentes. Basicamente são

apenas multiplicadores, somadores e elementos de retardo (delay). Desses,

multiplicadores e somadores implementam essas operações aritméticas em

sequências discretas. Existem basicamente duas principais classes de filtros

digitais:

12.1 Filtros do tipo FIR (Finite Inpulse Response)

A estrutura de um filtro FIR é bastante regular e, uma vez definidos os coeficientes,

o filtro pode ser completamente especificado. Esses são os coeficientes do filtro. Na

abaixo podemos ver uma estrutura simples de um filtro FIR. Podemos observar que

a passagem pelos componentes do filtro se dá sempre da esquerda para a direita.

Por isso, esse filtro é chamado também de feed-forward.

Página - 70

Suponha, baseado na figura acima, que o sistema tem uma entrada x[n] = [1, 0].

Sendo um sistema causal, a entrada para n < 0 é igual a zero. Assim, para n = 0,

temos:

E, para n = 1:

Logo, a saída seria y[n]=[0.5, 0.5]. A equação para cada termo é:

y[0] = 0.5.x[0] + 0.5.x[-1]

y[1] = 0.5.x[1] + 0.5.x[0]

Ou, de forma geral:

Os filtros FIR são expressos como:

com função de transferência:

Página - 71

A resposta ao impulso h[n] é dada por:

e a representação em equação de diferenças é:

Um filtro FIR pode ser representado em forma direta como:

Por simplicidade, pode-se representar um filtro FIR apenas com seus coeficientes.

Por exemplo, seja um filtro FIR com coeficientes [1, -1], ele pode ser representado

como:

a) Sistemas FIR com fase linear do Tipo I

Um sistema do tipo I tem resposta ao impulso simétrica

com M um inteiro par (observe que isso gera um número ímpar de amostras). O

atraso M/2 é um inteiro. A resposta em frequência é:

Considerando a condição de similaridade, essa resposta em frequência pode ser

expressa como:

Página - 72

onde a[0] = h[M/2],

a[k] = 2h[(M/2) – k], k = 1, 2, 3, ..., M/2

Assim, H(ejw) tem a forma de AP(ejw)e-jwM/2.

 Tipo I, resposta simétrica, M par

h[n] = 1, 0 £ n £ 4, e 0, caso contrário

h[n] = h[M – n]

A resposta em frequência é:

Magnitude e fase:

/Programa Filtro FIR I/

Exemplo 1: O seguinte código implementa no MatLab um filtro FIR tipo I:

h = [3 4 5 6 5 4 3]/30; % Veja que vai de zero a seis

M = length(h);

N = (M-1)/2;

L = 512;

H = fft([h zeros(1,L-M)]);

k = 0:L-1;

Página - 73

W = exp(j*2*pi/L);

A = H.* W.^(N*k);

A = real(A);

figure(1)

w = [0:L-1]*2*pi/(L-1);

subplot(2,1,1)

plot(w/pi,abs(H))

ylabel('|H(\omega)| = |A(\omega)|')

xlabel('\omega/\pi')

subplot(2,1,2)

plot(w/pi,A)

ylabel('A(\omega)')

xlabel('\omega/\pi')

Resultado:

Página - 74

b) Sistemas FIR com fase linear do Tipo II

Um sistema do tipo II tem uma resposta ao impulso simétrica como

Mas com M um inteiro ímpar. Nesse caso, H(ejw) pode ser expresso como:

Onde:

b[k] = 2.h[(M + 1)/2 - k], k = 1, 2, ..., (M + 1)/2

Novamente H(ejw) tem a forma de AP(ejw)e-jwM/2.

 Tipo II, resposta simétrica, M ímpar

h[n] = 1, 0 £ n £ 5, e 0, caso contrário

h[n] = h[M – n]

A resposta em frequência é:

Magnitude e fase:

Página - 75

/Programa Filtro FIR II/

Exemplo 2: O seguinte código implementa no MatLab um filtro FIR tipo II:

h = [3 5 6 7 7 6 5 3]/42;

M = length(h);

N = (M-1)/2;

L = 512;

H = fft([h zeros(1,L-M)]);

k = 0:L-1;

W = exp(j*2*pi/L);

A = H.* W.^(N*k);

A = real(A);

figure(2)

w = [0:L-1]*2*pi/(L-1);

subplot(2,1,1)

plot(w/pi,abs(H))

ylabel('|H(\omega)| = |A(\omega)|')

xlabel('\omega/\pi')

subplot(2,1,2)

plot(w/pi,A)

ylabel('A(\omega)')

xlabel('\omega/\pi')

Página - 76

Resultado

c) Sistemas FIR com fase linear do Tipo III

Se o sistema tem uma resposta ao impulso assimétrica:

h[n] = -h[M - n], 0 £ n £ M

com M um inteiro par, então H(ejw) tem a forma:

onde

c[k] = 2h[(M/2) - k], k = 1, 2, ..., M/2

Nesse caso, H(ejw) tem a forma:

Página - 77

 Tipo III, resposta assimétrica, M par

Se h[n] = δ[n] - δ[n – 2]

Então H(ejw) = 1 – e-2jw = j[2.sin(w/2)]e-jw

Magnitude e fase:

d) Sistemas FIR com fase linear do Tipo IV

Se o sistema tem uma resposta ao impulso assimétrica:

h[n] = -h[M - n], 0 £ n £ M

com M um inteiro ímpar, então H(ejw) tem a forma:

onde

d[k] = 2h[(M+1)/2 - k], k = 1, 2, ..., (M + 1)/2

e H(ejw) tem a forma:

Página - 78

 Tipo IV, resposta assimétrica, M ímpar Se h[n] = δ[n] - δ[n – 1]

Então H(ejw) = 1 – e-jw = j[2.sin(w/2)]e-jw/2

Magnitude e fase:

Página - 79

12. Filtros IIR

Os filtros FIR usam apenas cálculos feed forward. Se feedback é permitido, a

resposta de um filtro ao impulso não será necessariamente finita. Assim, filtros com

feedback são chamados de filtros com resposta ao Impulso Infinita (IIR – Infinite

Impulse Response). Por exemplo, o filtro abaixo:

cuja equação que descreve sua saída é dada por:

Se um impulso passa por esse filtro teremos como resposta as seguintes saídas,

considerando o filtro causal:

Ou seja, mesmo quando a entrada se anula, o filtro continua apresentando uma

saída. Essa saída diminui, mas não torna-se zero. Claro que, na prática, a resposta

chega a zero em algum momento. Considere então o filtro a seguir:

Página - 80

Esse filtro tem saídas:

Esse é um filtro IIR. A saída será sempre 0,8 mesmo a entrada permanecendo

0. Nesse próximo exemplo, a saída cresce mesmo com entrada zero.

De uma maneira geral, os filtros IIR são expressos como:

Página - 81

com a seguinte função de sistema:

Suas formas Direta I e Direta II são mostradas nas figuras abaixo

a) Forma direta I de um filtro IIR

Página - 82

b) Forma direta II de um filtro IIR

Exemplo 1: Considere a função de sistema:

A forma direta I e II podem ser desenhadas como:

Forma direta I:

Página - 83

e

Forma direta II:

Exemplo 2: Conexão em cascata

Considere a função de sistema:

Como os pólos e zeros são reais, uma estrutura em cascata tem seções com

coeficientes reais. Duas estruturas em cascata equivalentes podem ser criadas para

essa função:

Página - 84

Exemplo 3: Conexão em paralelo:

Considere a função de sistema (observe que é a mesma função anterior):

A forma paralela para esse sistema com uma seção de segunda ordem é:

Como os pólos são reais, podemos obter ainda uma forma paralela alternativa

expandindo H(z) como:

que gera o diagrama abaixo apenas com seções de primeira ordem:

Página - 85

13. Exercícios de Filtros

1. Um filtro IIR é definido por:

Desenhe sua estrutura na Forma Direta I ou na Forma Direta II.

2. Um filtro causal e invariante no tempo é definido por:

Desenhe sua estrutura na Forma Direta I ou na Forma Direta II. Esse filtro é FIR ou

IIR?

3. Um filtro IIR é definido por:

Desenhe sua estrutura na Forma Direta I ou na Forma Direta II

4. Modifique os códigos das dos filtros FIR tipo I e II para implementar exemplos de

filtros FIR dos tipos III e IV.

Página - 86

14. Projeto de Filtros FIR

Tanto a aproximação quanto a implementação podem ser realizadas de diversas

maneiras diferentes, com o resultado de que não existe uma solução única para o

problema de projeto de filtros com um conjunto prescrito de especificações.

Todavia, podemos mencionar três diferentes abordagens para o projeto de filtros

analógicos e digitais:

Abordagem analógica, a qual se aplica à classe de filtros analógicos.

Abordagem de analógico para digital, em que a motivação é projetar um

filtro digital lançando mão de um projeto de filtro analógico.

Abordagem digital direta a qual se aplica à classe de filtros digitais.

Para o projeto de filtros FIR, as técnicas são divididas nas seguintes categorias:

• Projeto usando janelas

• Método da amostragem em frequência

• Projeto equirriple ótimo

• Projeto de mínimos quadrados

14.1 Projeto usando janelas

A ideia básica de um projeto por janelas é selecionar um filtro seletor de frequências

ideal apropriado (que sempre é não-causal e de resposta ao impulso infinita) e

então truncar sua resposta ao impulso em uma janela para obter um filtro FIR causal

e de fase linear. Assim, o foco está na escolha de uma função de janelamento e um

filtro ideal apropriados.

Seja Hd(ejw) um filtro seletivo de frequência ideal que tem magnitude unitária e

características de fase linear sobre sua banda de passagem, e resposta zero na

banda de corte. Um filtro passa-baixa (FPB) ideal de largura de banda wc < π é

dado por:

Página - 87

onde wc é chamado de frequência de corte (cut-off) e a é chamado de atraso de

amostra (sample delay). A resposta ao impulso desse filtro é de duração infinita e é

dada por:

Para obter um filtro FIR a partir de hd[n], precisamos truncar hd[n] em ambos os

lados. Para obter um filtro FIR causal de fase linear h[n] de comprimento M,

devemos ter:

e

Essa operação é chamada de janelamento. Em geral, h[n] pode ser pensado como

sendo formado pelo produto de hd[n] e uma janela w[n] tal que:

h[n] = hd[n].w[n]

onde w[n] é alguma função simétrica com respeito a a no intervalo 0 £ n £ M – 1

e 0 fora desse intervalo. Dependendo de como obtivermos w[n] acima, temos

diferentes projetos de filtros.

Por exemplo:

é uma janela retangular.

Página - 88

No domínio da frequência, a resposta H(ejw) do filtro FIR causal é dada pela

convolução de Hd(ejw) e a resposta da janela W(ejw):

Podemos ver essa convolução na Fig. 5.3 para uma janela comum.

Fig. 15.1. Operação de janelamento no domínio da frequência.

Observações:

1. Como a janela w[n] tem comprimento finito igual a M, sua resposta em frequência

tem uma região de pico central (lóbulo principal) cuja largura é proporcional a 1/M e

tem lóbulos laterais com pesos menores.

2. A convolução gera uma versão da resposta ideal Hd(ejw), mas com algumas

distorções (ondulações).

3. A largura da banda de transição é proporcional a 1/M.

4. Os lóbulos laterais produzem ondulações que têm forma similar tanto na banda

de passagem quanto na de corte.

Projeto usando janelas: Para uma dada especificação de filtro, escolha um filtro de

comprimento M e uma função janela w[n] para a mais estreita largura do lóbulo

principal e a menor atenuação nos lóbulos laterais possível.

Da observação 4 acima, podemos notar que a tolerância d1 da banda de passagem

e a tolerância d2 da banda de corte não podem ser especificadas de forma

independente. Geralmente, toma-se d1 = d2.

Página - 89

Vamos descrever alguns tipos comuns de funções de janelamento.

a) Janela Retangular

Essa é a janela mais simples, mas que provê o pior desempenho em termos de

atenuação da banda de corte. Ela é definida como:

sendo sua resposta em frequência:

A magnitude da função sen[w(M + 1)/2]/sen(w/2) é mostrada na Fig. 5.4 para o caso

de M = 7. Note que W(ejw) tem fase linear generalizada. À medida que M aumenta,

a largura do lóbulo principal diminui.

Fig. 15.2. Magnitude da transformada de Fourier de uma janela retangular (M =7).

Página - 90

Pode ser provado que a largura do lóbulo central é Dwm = 4p/(M + 1) para uma

janela retangular. O primeiro zero de W(ejw) ocorre quando:

sen (w(M + 1)/2) = 0 ⇒ w(M + 1)/2 = p ⇒ 2p/(M + 1)

Assim, a largura do lóbulo central é o dobro desse valor (já que envolve os valores

negativos e positivos):

2w = 4p/(M + 1)

Observa-se também que a magnitude do primeiro lóbulo lateral é aproximadamente

em w = 3p/(M + 1) e é dada por:

À medida que M cresce, a largura de cada lóbulo lateral diminui, mas a área sobre

cada um permanece constante. Assim, as amplitudes relativas dos picos laterais

vão permanecer constantes e a atenuação da banda de passagem permanece em

cerca de 21 dB. Isso significa que as ondulações vão sofrer um pico perto das

bordas das bandas. Isso é conhecido como fenômeno de Gibbs (Fig. 5.5). Esse

fenômeno ocorre por causa da transição brusca de 0 para 1 (e de 1 para 0) da

janela retangular.

Página - 91

Fig. 15.3. Fenômeno de Gibbs: Pico das ondulações nas fronteiras entre as bandas.

b) Janela Triangular ou de Bartlett

Bartlett sugeriu uma transição mais suave para evitar o fenômeno de Gibbs. Isso

seria conseguido através de uma janela triangular da forma:

Essa janela e sua resposta em frequência podem ser vistas na Fig. 5.6

Página - 92

Fig. 15.4. Janela triangular.

c) Janela de Hanning

(homenagem a Julius von Hann, meteorologist austríaco)

Fig. 15.5. Janela de Hanning

d) Janela de Hamming (Richard Hamming, matemático americano)

Tem uma quantidade menor de descontinuidades em relação à janela de Hanning.

Página - 93

Fig. 15.6. Janela de Hamming.

e) Janela de Blackman

Também é similar às duas anteriores, mas tem um segundo harmônico o que faz

com que ela se aproxime de zero com mais suavidade.

Fig.15.7. Janela de Blackman.

Tanto a janela de Bartlett, quanto Hamming, Hanning e Blackman têm lóbulos

laterais menores do que os da janela retangular. No entanto, para o mesmo valor de

M, a largura do lóbulo principal também é mais larga para essas janelas se

comparadas à janela retangular. Consequentemente, essas janelas conseguem

Página - 94

uma convolução no domínio da frequência mais suave e, como resultado, a região

de transição na resposta do filtro FIR é mais larga.

Para reduzir a largura da região de transição podemos aumentar o comprimento da

janela, o que resulta em filtros mais largos.

A tabela a seguir resume algumas características no domínio da frequência dessas

janelas.

f) Janela de Kaiser (James F. Kaiser)

Esta é a melhor janela. Ela e considerada ótima porque provê um lóbulo principal

largo para a dada atenuação da banda de corte, o que implica a mais brusca banda

de transição. A função foi definida por Kaiser e é dada por:

I0(.) é a função de Bessel modificada de ordem zero:

Página - 95

Fig. 15.8. Variadas formas da Janela de Kaiser.

Na expressão de w[n], existem dois parâmetros:

1. O comprimento M

2. O parâmetro b

Variando b e M, é possível ajustar a amplitude dos lóbulos laterais. Kaiser encontrou

duas fórmulas que permitem achar M e b de modo a atender às especificações do

filtro. Assim, dado que d1 é fixo (especificado), a frequência de corte wP da banda

de passagem do filtro passa-baixa é a maior frequência tal que:

|H(ejw)| ³ 1 - d1

A frequência da banda de corte tem tolerância d2, satisfazendo:

|H(ejw)| £ d2

A largura da banda de transição é:

Dw = wS - wP

Dado:

A = -20log10 d,

considerando d1 = d2. Kaiser mostrou que:

Página - 96

Além disso, dados Dw e A, M é aproximadamente:

O procedimento para projetar um filtro passa-baixa digital FIR usando a janela de

Kaiser consiste nos seguintes passos:

i) Estabelecer as especificações wP, wS e d.

ii) Estabelecer a frequência de corte wc do filtro passa-baixa ideal ao qual se

aplicará a janela (wc = (wP + wS)/2).

iii) Calcular A = 20log10 d e Dw = wP - wS e usar as fórmulas de Kaiser para

encontrar os valores de M e b.

iv) Encontra a resposta ao impulso do filtro através de h[n]=hd[n]w[n], onde w[n] é a

janela de Kaiser e hd[n] = Á-1[Hd(ejw)].

Devido à complexidade de cálculos com funções de Bessel, o projeto dessas

janelas não é fácil. A equação de w[n] definida por Kaiser tem valores encontrados

empiricamente e são definidos sem prova.

Exemplo: Projetar, usando janelas de Kaiser, um filtro passa-baixa com as

seguintes especificações: wP = 0,4p, wS = 0,6p e d = 0,001.

wc = (wS + wP)/2 = 0,5p

Dw = wS - wP = 0,2p

A = -20log10 d = 60 dB

Como A > 50, pela Eq. 5.1:

b = 0,1102(A – 8,7) @ 5,633

M = (A - 8)/(2,285Dw) @ 36,219 ⇒ M = 37 (M inteiro)

A resposta ao impulso é:

com w[n] dado pela definição da janela de Kaiser.

Página - 97

Implementações no MatLab

O MatLab tem diversas funções para implementar janelas:

1. w = rectwin(M): Janela retangular

2. w = bartlett(M): Janela de Bartlett

3. w = hanning(M): Janela de Hanning

4. w = hamming(M): Janela de Hamming

5. w = blackman(M): Janela de Blackman

6. w = kaiser(M, beta): Janela de Kaiser

Antes de projetarmos alguns exemplos, vamos implementar duas funções base

importantes para os exemplos a seguir. Uma implementa uma resposta ao impulso

ideal de um filtro passa-baixa hd[n]. A outra função implementa a plotagem no

domínio da frequência, apresentando também a resposta em magnitude absoluta e

em escala dB (é uma variação da função freqz do MatLab).

Função 1:

function hd = ideal_lp(wc, M)

% Ideal low pass filter

% wc = cutoff frequency

% M = length of the ideal filter

alpha = (M - 1)/2

n = [0:(M-1)];

m = n - alpha + eps;

hd = sin(wc*m)./(pi*m);

Função 2:

function [db, mag, pha, w] = freqz_m(b, a)

% Versao modificada da funcao freqz

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

H = (H(1:501))';

w = (w(1:501))';

mag = abs(H);

db = 20*log10((mag + eps)/(max(mag)));

pha = angle(H);

Página - 98

Exemplo 1:

Projetar um filtro passa-baixa FIR com as seguintes especificações wP = 0,2p, RP =

0,25 dB, wS = 0,3p e AS = 50 dB.

Tanto a janela de Hamming quanto a de Blackman provêem atenuação de mais de

50 dB. Vamos escolher a janela de Hamming que provê a menor banda de transição

e assim tem a menor ordem.

wp = 0.2*pi; ws = 0.3*pi;

tr_width = ws - wp;

M = ceil(6.6*pi/tr_width) + 1

n = [0:M-1];

wc = (ws + wp)/2;

hd = ideal_lp (wc, M);

w_ham = (hamming(M))';

h = hd.*w_ham;

[db, mag, pha, w] = freqz_m(h, [1]);

delta_w = 2*pi/1000;

Rp = -(min(db(1:wp/delta_w+1)))

As = -round(max(db(ws/delta_w+1:501)))

subplot(1, 1, 1)

subplot (2, 2, 1); stem(n, hd); title('Resposta ao Impulso Ideal');

axis([0 M-1 -0.1 0.3]);xlabel('n');ylabel('hd[n]');

subplot (2, 2, 2); stem(n, w_ham); title('Janela de Hamming');

axis([0 M-1 0 1.1]);xlabel('n');ylabel('w[n]');

subplot (2, 2, 3); stem(n, h); title('Resposta ao Impulso Atual');

axis([0 M-1 -0.1 0.3]);xlabel('n');ylabel('h[n]');

subplot (2, 2, 4); plot(w/pi, db); title('Magnitude em dB');grid

axis([0 1 -100 10]);xlabel('frequencia em pi unidades');ylabel('Decibeis');

M = 67

alpha = 33

Rp = 0,0394

As = 52

Página - 99

Exemplo 2:

Resolva o exemplo anterior usando uma janela de Kaiser.

wp = 0.2*pi; ws = 0.3*pi; As = 50;

tr_width = ws - wp;

M = ceil((As - 7.95)/(14.36*tr_width/(2*pi))+1) + 1

n = [0:M-1];

beta = 0.1102*(As - 8.7)

wc = (ws + wp)/2;

hd = ideal_lp (wc, M);

w_kai = (kaiser(M, beta))';

h = hd.*w_kai;

[db, mag, pha, w] = freqz_m(h, [1]);

delta_w = 2*pi/1000;

As = -round(max(db(ws/delta_w+1:501)))

subplot(1, 1, 1)

subplot (2, 2, 1); stem(n, hd); title('Resposta ao Impulso Ideal');

axis([0 M-1 -0.1 0.3]);xlabel('n');ylabel('hd[n]');

subplot (2, 2, 2); stem(n, w_kai); title('Janela de Kaiser');

axis([0 M-1 0 1.1]);xlabel('n');ylabel('w[n]');

subplot (2, 2, 3); stem(n, h); title('Resposta ao Impulso Atual');

comentários (0)
Até o momento nenhum comentário
Seja o primeiro a comentar!
Esta é apenas uma pré-visualização
Consulte e baixe o documento completo
Docsity is not optimized for the browser you're using. In order to have a better experience we suggest you to use Internet Explorer 9+, Chrome, Firefox or Safari! Download Google Chrome