Matlab-Geração-Sinais-Interação-com-Sistemas

Parte 1 Aprendendo a trabalhar com sinais e a interagir com sistemas Esta prática tem por propósito ilustrar a geração de sinais básicos e a trabalhar com eles com as principais funções do Matlab Ponto 1 geração de sinais básicos Amostra unitária δn esta função é gerada usando o código impseqn0 inicio fim conforme exemplo sinall n1 impseq0 3 3 sinal1 0 0 0 1 0 0 0 n1 3 2 1 0 1 2 3 sinal2 n2 impseq3 1 5 sinal2 0 0 1 0 0 n2 1 2 3 4 5 Degrau unitário un esta função é gerada usando o código stepseqn0 inicio fim conforme exemplo sinall n1 stepseq0 3 3 sinal1 0 0 0 1 1 1 1 n1 3 2 1 0 1 2 3 sinal2 n2 stepseq 3 1 5 sinal2 0 0 1 1 1 n2 1 2 3 4 5 Sequência exponencial complexa para entender melhor o comportamento de sinais experimente gerar um e visualizálo conforme exemplo abaixo da função n 530 aexp023jn subplot221 stemn reala titleParte real subplot222 stemn imaga titleParte imaginária subplot223 stemn absa titleMagnitude subplot224 stemn anglea titleFase bexp02n cexp3jn figure subplot211 stemn b titleMódulo subplot212 stemn c titleFase Ponto 2 Amostrando sinais Para simular o processo de conversão de tempo contínuo para tempo discreto utilize o código abaixo inicializacoes F 60 freq analogica Fs 400 freq aquisicao ts 1Fs tempo aquisicao numciclos 3 equivale a 3 ciclos t 01e6numciclos1F tempo simulacao n 0roundnumciclosFsF numero amostras parte analogica fanalogico sin2piFt plottfanalogico parte digital fdigital sin2piFFsn hold on stemntsfdigital Ponto 3 Réplicas espectrais A ideia deste ponto é mostrar a réplica espectral Para isto considere um sistema cuja taxa de aquisição é Fs6000 amostras por segundo Se considerarmos a frequência de 1kHz as frequências de 7k 1kFs 13k1k2Fs 19k1k3Fs e assim por diante passam pelos mesmos pontos que o sinal de 1kHz se confundindo com eles e gerando o indesejado fenômeno de aliasing inicializacoes F1 1000 freq analogica Fs 6000 freq aquisicao ts 1Fs equivale a 3 ciclos numciclos 1 t 01e6numciclos1F1 tempo simulacao n 0roundnumciclosFsF1 numero amostras parte analogica fanalogico1 sin2piF1t fanalogico2 sin2pi7000t fanalogico3 sin2pi13000t plottfanalogico1r hold on plottfanalogico2b plottfanalogico3k legend1k7k13k parte digital fdigital sin2piF1Fsn stemntsfdigital Ponto 4 operação com sinais Soma a soma de dois sinais em função de n pode ser feita usando a função sigaddsinal1 ndosinal1 sinal2 ndosinal2 conforme exemplo abaixo x1 1 2 3 2 1 n1 2 1 0 1 2 x2 2 3 4 n2 0 1 2 x3 n3 sigaddx1n1x2n2 stemn3x3 Multiplicação a multiplicação de dois sinais em função de n pode ser feita usando a função sigmultsinal1 ndosinal1 sinal2 ndosinal2 conforme exemplo abaixo x1 1 2 3 2 1 n1 2 1 0 1 2 x2 2 3 4 n2 1 0 1 x3 n3 sigmultx1n1x2n2 Universidade Federal de Uberlândia Faculdade de Engenharia Elétrica Prof Alan Petrônio Pinheiro wwwalanengbr alaneletricaufubr Curso de EngEletrônica e Telecomunicações campus Patos de Minas Disciplina de Processamento Digital de Sinais Versão documento 13 3 6 stemn3x3 Deslocamento o deslocamento avanço ou atraso de um sinal pode ser gerado através da função sigshift conforme mostra o exemplo 1 2 3 4 5 x1 1 2 3 2 1 n1 0 1 2 3 4 x2 n2 sigshiftx1n12 mesmo que x1n2 subplot211 stemn1x1 titleOriginal subplot212 stemn2x2 titleAtrasada 2 unid Reflexão a reflexão em torno do eixo vertical é dada por ynxn e é implementada conforme mostra o exemplo 1 2 3 4 x1 1 2 3 4 5 n1 0 1 2 3 4 x2 n2 sigfoldx1n1 stemn2 x2 Convolução embora o Matlab tenha uma função própria para convolução a conv ela leva em conta que os sinais são sempre casuais ou seja xn0 para n0 Para evitar esta limitação é aconselhável usar a função convm conforme ilustra exemplo abaixo 1 2 3 4 5 6 x1 3 11 7 0 1 4 2 n1 33 h 2 3 0 5 2 1 nh 14 y ny convmx1n1hnh stemny y Parte 2 DFT Nesta prática pretendese demonstrar o uso da transformada discreta de Fourier DFT e a principal função do Matlab para cálculo da transformada discreta de Fourier a FFT Ponto 1 calcular manualmente a transformada Para calcular a DFT de um sinal manualmente utilize o código abaixo Vale destacar que neste caso consideramos o sinal como causal em n0 Note que neste exemplo não se nota a preocupação em otimizar o código e se pode usar tanto a forma de Euler linha 10 ou pela sua forma expandida linha 11 Ambos casos conduzem a um mesmo resultado Código 31 Cálculo manual da DFT clear clc sinal 1 2 3 4 5 sinal Fs 1000 taxa de aquisição N sizesinal2 num de amostras y zeros1N calcular a DFT for m0N1 acumulador 0 for n0N1 acumulador acumulador sinaln1cos2pinmNjsin2pinmN acumulador acumulador sinaln1expj2pinmN end ym1acumulador end calcular eixo frequencias f m 0N1 f mFsN plotar sinais magX absy angX angley realX realy imagX imagy subplot221 stemfmagX titlemagnitude subplot223 stemfangX titlefase subplot222 stemfrealX titlereal subplot224 stemfimagX titleimaginario Ponto 2 fft O código ilustrado anteriormente é capaz de calcular a DFT mas não tem eficiência computacional pois não faz uso das propriedades de simetria de cálculo da DFT Uma função computacionalmente mais eficiente para cálculo da DFT é a fft fast Fourier transform O código 32 ilustra seu uso Código 32 Cálculo manual da DFT usando a função FFT do Matlab clc clear gera um sinal sintetico com 3 componentes de frequencia Fs 1000 taxa de aquisicao do sinal ts 1Fs N 1000 numero amostras de amostras analisadas t 0N1ts sinal 06sin2pi50t 2sin2pi120t 03sin2pi400t ruido 04randnsizet calcula a DFT usando a funcao fft y fftsinalruido y yN se desejar podese dividir por N para acomodar os calculos calcular o eixo das frequencias m 0N1 f mFsN y y1N2 pegando so a primeira metade já que a outra é cópia 16 f f1N2 pegando so a primeira metade já que a outra é cópia 17 calcular a potência do espectro em db 18 magnitude absy 19 fase angley 20 fref maxmagnitude escolhe o maior valor para ser a referencia para normalizacao 21 ydb 20log10magnitudefref lembre que por exemplo 0db 1 unidade 22 plotar 23 subplot 311 plotf magnitude titlemagnitude espectro xlabelfreqHz 24 ylabelAmplitude 25 subplot 312 plotf ydb titlepotencia espectro 26 subplot 313 plotf fase titlefase Código 33 Cálculo manual da DFT usando a função FFT do Matlab 1 N 100 2 janela1 windowblackmanharrisN 3 janela2 windowhammingN 4 janela3 windowgausswinN25 5 wvtooljanela1 janela2 janela3 bartbamnwin bartlett blackman blackmanharris bohmanwin chebwin flattopwin gausswin hamming hann kaiser nuttallwin parzenwin rectwin taylorwin triang tukeywin Para ilustrar o uso de uma janela de Hanning podemos proceder como no código 34 que é uma versão adaptada do código 32 Veja a diferença da DFT de um sinal janelado e um puro não janelado Código 34 Uso de janelas para cálculo da DFT visando diminuir o leakage 1 clc clear 2 gera um sinal sintetico 3 Fs 200 4 ts 1Fs 5 N 300 6 t 0N1ts 7 sinal 06sin2pi13t 8 cria um segundo sinal janelado 9 janela windowhannN 10 sinaljan sinaljanela 11 calcula a DFT usando a funcao fft 12 y fftsinal 13 yjan fftsinaljan 14 calcular o eixo das frequencias 15 m 0N1 16 f mFsN 17 y y1N2 18 yjan yjan1N2 19 f f1N2 20 calcular a potência do espectro em db 21 magnitude absy 22 magnitudejan absyjan 23 fref maxmagnitude 24 frefjan maxmagnitudejan 25 ydb 20log10magnitudefref 26 ydbjan 20log10magnitudejanfrefjan 27 plota 28 subplot 321 plott sinal titlesinal nao janelado 29 subplot 323 stemf magnitude titlemagnitude espec sem jan ylim0 maxmagnitude 30 subplot 325 plotf ydb titlepotencia espectro sem janela 31 subplot 322 plott sinaljan titlesinal janelado 32 subplot 324 stemf magnitudejan titlemagnitude espec com jan ylim0 maxmagnitude 33 subplot 326 plotf ydbjan titlepotencia espectro com janela Ponto 4 diferença entre resolução espectral e densidade espectral Para ilustrar a diferença entre resolução espectral e densidade espectral utilizemos um exemplo Considere o sinal xncos048πncos052πn Calcula a DFT do sinal considerando A um número de amostras N 10 e preencha as outras 90 amostras com zero chegando a N100 B um número de amostras N 100 sem preenchimento com zero Código 35 Código para solução dos dois itens anteriores Observe que é usada a função DFT da biblioteca empregada neste apostila 1 clc clear 2 n 099 3 Fs 100 4 x cos048pin cos052pin 5 calcula DFT com N10 e plota 6 n1 09 7 sinal1 x110 8 Y1 dftsinal1 10 9 magY1 absY116 10 N1 5 11 m1 05 12 f1 m1FsN1 13 subplot231 stemsinal1 xlabeln titleN10 14 subplot234 stemf1 magY1 xlabelfreqHz 15 calc DFT com N1090 zeros e plota melhor resol espectro 16 sinal2 x110 zeros190 preencheu sinal que essencialmente é o mesmo 17 Y2 dftsinal2100 18 magY2 absY2150 19 N2 50 20 m2 150 21 f2 m2FsN2 22 subplot232 stemsinal2 xlabeln titleN1090 zeros com boa resol espec 23 subplot235 stemf2 magY2 xlabelfreqHz 24 calc DFT com N100 25 sinal3 x1100 26 Y3 dftsinal3100 27 magY3 absY3150 28 N3 50 29 m3 150 30 f3 m3FsN3 31 subplot233 stemsinal3 xlabeln titleN100 com boa densid espec 32 subplot236 stemf3 magY3 xlabelfreqHz Parte 3 Projeto de filtros FIR Nesta prática pretendese demonstrar i o projeto de filtros de resposta finita FIR finite impulse response através do cálculo de seus coeficientes e ii a implementação do filtro Pra isto usamos a técnica de janelas Resumo pontos abordados nesta parte Ponto 1 Execução de uma filtragem usando função conv Ponto 2 Calculando coeficientes FIR fase a fase Ponto 3 Calculando coeficientes FIR usando funções fir1 e fir2 Ponto 4 Calculando coeficientes FIR usando ParksMcClellan Ponto 1 Implementação de filtros FIR Para implementar um filtro FIR é necessário que se tenha em mãos os coeficientes do filtro Por isto antes é necessário aprender a projetar um filtro que é o que faremos nos próximos tópicos Por enquanto vamos imaginar que já temos os coeficientes prontos e nossa intenção se resume a implementar o filtro considerando que o projeto que é a designação dos coeficientes do filtro já esteja pronta O Código 41 mostra como é feita a implementação de um filtro FIR passa baixas com coeficientes já dados liha 2 e 3 Note que a execução do filtro é feita por uma convolução da entrada com a resposta impulsiva hn do filtro na linha 17 Código 41 Implementação de um filtro FIR dados seus coeficientes 1 clc clear 2 h 04375 03142 00625 00935 00625 00418 00625 00124 00625 00124 3 00625 00418 00625 00935 00625 03142 4 nh 015 5 PARTE 1 Gera um sinal analógico entradaanalogica e sua versão discretizada entradadiscretizada 6 Fs 20000 taxa de amostragem 20mil amostrasseg 7 t 01Fs002 amostragem de 002seg 8 N maxt1Fs 9 n 0N quantidade de amostras da entrada para filtrar 10 entradaanalogica sin2pi 200t sin2pi25t sin2pi 5000t sin2pi 11 9000t 12 entradadiscretizada sin2pi200nFs sin2pi25nFs sin2pi5000nFs sin2pi9000nFs 13 14 PARTE 2 convolui entrada com h 17 saida conventradadiscretizada h 18 19 PARTE 3 PLOTA RESULTADOS 20 calcula espectro entrada 21 fftsinalentrada fftentradadiscretizadaN 22 entrada nFsN 23 calcula espectro saida 24 N3 sizesaida2 25 fftsinalsaida fftsaidaN3 26 n3 0sizefftsinalsaida21 27 fsaida n3FsN3 28 calcula espectro filtro 29 fftrespfiltro ffth 30 n4 0sizefftrespfiltro21 31 N4 sizen42 32 fh n4FsN4 33 plota espectros 34 subplot311 stemnhh xlabeln titleCoef filtro 35 subplot312 plotfentrada1N2 absfftsinalentrada1N2 36 hold on 37 plotfsaida1N32 absfftsinalsaida1N32r 38 plotfh1N42 absfftrespfiltro1N42g 39 legendentradasaídaresp impulsiva 40 xlabelFreq Hz 41 titleEspectros dos sinais e filtro 42 plota sinais 43 subplot313 plotn entradadiscretizada 44 hold on 45 plotn3 saidar 46 legendentradasaída 47 xlabeln 48 titleSinais discretos do filtro Uma outra forma muito usual de se executar uma filtragem no Matlab é usando a função filter Para isto basta substituir a linha código saida conventradadiscretizada h do exemplo anterior por saida filterh1entradadiscretizada Ponto 2 Projetar um filtro FIR pelo método de janelas Projetar um filtro FIR usando o método de janelas é uma tarefa relativamente fácil No exemplo do código 42 dividimos o código em oito partes para facilitar seu entendimento Na parte 1 linhas 3 a 7 especificamos as características do filtro São elas o Fs indica a taxa de aquisição máxima do sinal de entrada do filtro o Fc indica a frequência de corte em Hertz o N número de coeficientes do filtro que pretendemos projetar o NumCoeffiltro dos N coeficientes calculados devemos usar apenas parte deles que será designado por esta variável NumCoeffiltro Na parte 2 linhas 10 a 16 especificamos o espectro ideal para o filtro baseado nos valores coletados na parte 1 A partir deste espectro idealizado encontramos na parte 3 linha 19 quais devem ser os coeficientes ideais do filtro para implementação deste espectro idealizado Note que se o filtro é passa alta linha 20 devemos multiplicar os coeficientes encontrados por cos2πnFs2Fscosπn1111 Isto transforma o filtro passabaixas em passaaltas Já na linha 25 o código testa se foi escolhido um filtro passafaixas para transformar o filtro passabaixas em um passafaixas através do deslocamento de frequência feito pela multiplicação dos coeficientes do filtro passabaixa por cos2πnfcFs onde fcc indica a frequência de centro da banda de passagem Na parte 4 cuidados para gerar um gráfico de coeficientes simétrico em relação ao eixo vertical pois a ifft só estima metade da função sync A outra metade deve ser preenchida através de espelhamento Na parte 5 temos a opção de usar uma janela aos coeficientes para atenuar o ripple da banda passante Aplicar uma janela é importante quando NumCoeffiltro N pois nesta situação é gerado um ripple que pode ser atenuado por algumas funções de janelamento Finalmente nas partes 6 a 8 geramos um sinal sintético com 4 componentes de frequência para testar o filtro e calculamos os espectros do sinal sintético de entrada o espectro da resposta impulsiva do filtro e o espectro do sinal de saída Faça testes para ver os efeitos de filtragem Código 42 Projeto de um filtro FIR passabaixas passalatas e faixas 1 clc clear Universidade Federal de Uberlândia Faculdade de Engenharia Elétrica Prof Alan Petrônio Pinheiro wwwalanengbr alaneletricaufubr Curso de EngEletrônica e Telecomunicações campus Patos de Minas Disciplina de Processamento Digital de Sinais Versão documento 13 9 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 PARTE 1 dados do filtro Fs 20000 taxa de amostragem fc 2000 frequencia de corte do filtro N 64 Quantidade total de coeficientes filtro NumCoefFiltro 64 Qtos coeficientes vou usar do total tipofiltro input Tipo filtro1passabaixas 2passaaltas 3passafaixa PARTE 2 gera o espectro ideal do filtro H zeros1N resposta em freq ideal do filtro fresol FsN calcula resolução frequencia mcorte roundfcfresol1 estima indice m de fc H ones1 mcorte 1 zeros1 N2mcorte zeros1 N2mcorte ones1 mcorte 1 gera espectro ideal do filtro m 0N1 f mfresol define o eixo de frequencias do grafico 1 PARTE 3 gera todos coef ideais do filtro hideal ifftH if tipofiltro2 se passaalta multiplicar coef por 1111 n 0N1 deslocamentof cos2pinFs2Fs hideal hidealdeslocamentof end if tipofiltro3 se passafaixa multiplicar coef por modulacao n 0N1 fcc input Digite em Hz frequência central deslocamentof cos2pifccnFs hideal hidealdeslocamentof end PARTE 4 GERA SIMETRIA PAR NA FUNCAO SYNC DOS COEFICIENTES hidealN2N hideal1N21 for i2N2 hideali1 hidealNi1 end inicio N2 NumCoefFiltro2 1 fim N2 NumCoefFiltro2 h realhidealiniciofim pega so parte dos coeficientes ideias do filtro PARTE 5 opcional aplica janela resposta input Deseja aplicar janela aos coef 1sim 2nao if resposta 1 jan windowblackmanfiminicio1 h hjan end PARTE 6 testa a implementação filtro com sinal sintetico Nsinalsintetico 400 quantidade de amostras do sinal sintetico n1 0Nsinalsintetico1 entradadiscretizada sin2pi500n1Fs sin2pi2500n1Fs sin2pi5000n1Fs sin2pi8000n1Fs saida conventradadiscretizada h Nsinalsaida sizesaida2 Nrespimpulsitva NumCoefFiltro PARTE 7 calcula espectros sinal entrada saida e hn fftsinalentrada fftentradadiscretizadaNsinalsintetico fftsinalsaida fftsaidaNsinalsaida fftrespfiltro ffth fentrada n1FsNsinalsintetico n3 0sizefftsinalsaida21 fsaida n3FsNsinalsaida n2 0sizefftrespfiltro21 fh n2FsNrespimpulsitva PARTE 8 plota subplot221 stemfH titleHf idealizado xlabelfHz subplot222 stemrealhideal xlabeln hold on steminiciofimhr titleCoeficientes hn do filtro ylabelAmplitude xlabeln legendTodosSelecionados subplot223 plotfentrada1Nsinalsintetico2 absfftsinalentrada1Nsinalsintetico2 hold on plotfsaida1Nsinalsaida2 absfftsinalsaida1Nsinalsaida2r plotfh1Nrespimpulsitva2 absfftrespfiltro1Nrespimpulsitva2g 79 legendentradasaídaresp impulsiva 80 xlabelFreq Hz 81 titleEspectros dos sinais e filtro 82 subplot224 plotn1 entradadiscretizada 83 hold on 84 plotn3 saidar 85 legendentradasaída 86 xlabeln 87 titleSinais discretos do filtro Uma forma muito interessante de se avaliar a resposta impulsiva de um filtro é usando a função freqz do Matlab O argumento desta função são os coeficientes hn do filtro e o resultado é a magnitude da resposta em frequência do filtro designado como mag no código 43 e as correspondentes frequências normalizadas por π Assim a maior frequência tem valor 1 e indica a frequência de Nyquest que em nosso exemplo é Fs220000210000Hz Código 43 Visualização da resposta em frequência de um filtro 1 magf freqzh 2 plotfpiabsmag Ponto 3 Projetar um filtro FIR usando a função do Matlab FIR1 e FIR2 A função FIR1 do Matlab tem um dos seguintes formatos coeficientes fir1nWn coeficientes fir1nWnftype coeficientes fir1nWnwindow coeficientes fir1nWnftypewindow Onde Wn indica a faixa de frequências que determinam o filtro n a ordem ftype o tipo que pode ser high ou stop e window os coeficientes da janela que pode ser opcional Observe os exemplos do código 44 A função FIR1 usa o método de projeto conhecido como projeto baseado por resposta ao impulso Código 44 Cálculos coeficientes filtro FIR usando função FIR1 1 ordem 50 2 resolplotfreq 512 3 4 Exemplo 1a passabaixas frequência normalizada 5 coef fir1ordem 03 6 freqzcoef1resolplotfreq 7 8 Exemplo 1b passabaixas frequência nominal considerando Fs20k 9 fc 3000 10 Fs 20e3 11 coef fir1ordem fcFs2 12 H freq freqzcoef1resolplotfreq 20e3 13 plotfreqabsH 14 15 16 Exemplo 2 passaaltas 17 coef fir1ordem 04high 18 freqzcoef1resolplotfreq 19 20 21 Exemplo 3 passaaltas com janela de Hanning 22 coef fir1ordem 06 high hannordem1 freqzcoef1resolplotfreq Exemplo 4 passabanda coef fir1ordem 03 05 freqzcoef1resolplotfreq Exemplo 5 rejeitabanda coef fir1ordem 03 07 stop freqzcoef1resolplotfreq No código 44 devese escolher frequências normalizaas Isto quer dizer que a maior frequência que pode ser representada pelo sistema que equivale a Fs2 segundo critério de Nyquest é denominada de 1 Assim se desejamos projetar um filtro para cortar em 3000Hz usando um sistema com Fs20000 sabese que a máxima frequência de 10000Hz equivale a 1 Logo a frequência de 3000Hz equivale 300010000 03 A função freqz analisa a resposta em frequência dos coeficientes Um outro tipo de projeto dos coeficientes do filtro FIR é conhecido como amostragem baseada em frequência e é implementada pela função FIR2 conforme ilustra o código 45 A função FIR2 tem geralmente o formato fir2nfmwindow onde n indica o ordem do filtro f as frequências do filtro m as respectivas amplitudes que devem ter as frequências especificadas pelo vetor f e window o tipo de janela usado para os n1 coeficientes do filtro o uso de janelas é opcional e pode ser deixado em branco conforme ilustra alguns dos exemplos do código 45 Código 45 Cálculos coeficientes filtro FIR usando função FIR2 1 ordem 50 2 resolplotfreq 512 3 4 Exemplo 1 passabaixas 5 f 0 03 04 1 m 1 1 0 0 6 coef fir2ordem f m 7 freqzcoef1resolplotfreq 8 9 Exemplo 2 passaaltas 10 f 0 05 065 1 m 0 0 1 1 11 coef fir2ordem f m 12 freqzcoef1resolplotfreq 13 14 Exemplo 3 passaaltas com janela 15 f 0 05 065 1 m 0 0 1 1 16 coef fir2ordem f m hannordem1 17 freqzcoef1resolplotfreq 18 19 Exemplo 4 passabanda 20 f 0 04 05 07 08 1 m 0 0 1 1 0 0 21 coef fir2ordem f m 22 freqzcoef1resolplotfreq 23 24 Exemplo 5 rejeitabanda 25 f 0 04 05 07 08 1 m 1 1 0 0 1 1 26 coef fir2ordem f m 27 freqzcoef1resolplotfreq Ponto 4 Projetar um filtro FIR pelo método de ParksMcClellan É uma das técnicas mais usadas atualmente para se encontrar coeficientes pois é uma técnica otimizada Seu algoritmo é complexo por isto vamos nos ater somente a ensinar como proceder para usar as funções do Matlab para encontrar os valores dos coeficientes do filtro FIR Para isto devemos lembrar os principais parâmetros de um filtro como ilustra a figura abaixo Universidade Federal de Uberlândia Faculdade de Engenharia Elétrica Prof Alan Petrônio Pinheiro wwwalanengbr alaneletricaufubr Curso de EngEletrônica e Telecomunicações campus Patos de Minas Disciplina de Processamento Digital de Sinais Versão documento 13 12 Baseados no que vemos na figura anterior devemos especificar estes principais parâmetros para se determinar os coeficientes Para isto podemos usar a função firpm do Matlab como ilustra o Código 44 Para projetar um filtro passafaixas devemos usar as especificações das linhas 2 e 3 e um filtro passabaixas as linhas 6 e 7 Utilize a função freqz para ver os resultados e entender como funciona a especificação destes filtros Lembrese que a frequência deve ser normalizada Ou seja a frequência 1 indica a frequência de Nyquest Código 46 Ilustração do uso da função firpm que executa o algoritmo de ParksMcClellan para encontrar coeficientes de filtro FIR 1 2 3 4 5 6 7 8 9 10 coeficientes para filtro passafaixa f 0 04 044 056 06 1 a 0 0 1 1 0 0 coeficientes para filtro passabaixas f 0 04 05 1 m 1 1 0 0 coef firpm32fm freqzcoef Note no exemplo anterior que escolhemos 32 coeficientes Mas seria este valor correto Geralmente o projeto de um filtro envolve as especificações da atenuação ripple e banda de transição Assim devemos inserir estes valores para descobrir de quantos coeficientes precisamos para conseguir estes valores Assim o modo mais eficiente de se usar a função firpm é antes determinar os parâmetros ótimos do filtro usando a função firpmord como exemplifica o código 45 Para isto considere segundo figura anterior wc1500Hz ws2000Hz δ1001 e δ201 com Fs8000 amostras por segundo Note que usando estes argumentos usamos a função firpmord para encontrar os melhores valores para a ordem do filtro e os pontos chaves de frequência do filtro Nas linhas 4 e 5 plotamos a resposta em frequência do filtro Código 47 Uso otimizado do filtro FIR usando ParksMcClellan Exemplo de um filtro passabaixas 1 2 3 Fs 8000 nfoaow firpmord1500 20001 0001 01Fs coef firpmnfoaow freqzcoef1512 No código da sequência é ilustrado o emprego da mesma função para um filtro passabanda Observe que diferentemente das funções fir1 e fir2 a função firpmord emprega um modo diferente para representar os ganhos e atenuações nas bandas de passagem ou rejeição No caso do exemplo abaixo Universidade Federal de Uberlândia Faculdade de Engenharia Elétrica Prof Alan Petrônio Pinheiro wwwalanengbr alaneletricaufubr Curso de EngEletrônica e Telecomunicações campus Patos de Minas Disciplina de Processamento Digital de Sinais Versão documento 13 13 temos três bandas A primeira delas é de rejeição pois o ganho é próximo de 0 conforme especificado no código e vai de 0 a 1500Hz A segunda é de passagem pois tenho ganho 1 e vai de 2000 a 3000 Hz E a terceira banda é de rejeição e vai de 3500 Hz até Fs2 80002 4000 Hz Código 48 Exemplo de um filtro passabanda ParksMcClellan 1 2 3 4 5 6 7 8 9 10 11 12 13 Fs 8000 k1 01 atenuacao primeira banda de rejeicao k2 005 ripple da banda de passagem k3 01 atenuacao segunda banda de rejeicao n fo ao w firpmord1500 2000 3000 3500 0 1 0 k1 k2 k3 Fs b firpmn fo ao w calcula o eixo da freq espectro w freqzb1512 retorna o espectro em espectro e o eixo de frequen em w freq wpi normaliza as frequencias de 0 ate 1 onde 1Fs2 w está de 0 até pi freq freqFs2 pego o vetor freq normalizado e multiplico pela máx freq que é Fs2 plota subplot211 plotfreq 20log10absespectro titleescala dB grid on subplot212 plotfreq absespectro titleescala linear grid on No exemplo anterior podese perceber que no vetor de frequências deve ser designado para cada transição as duas frequências que designam a banda de transição Por exemplo para projetar o filtro da figura abaixo devese observar que ele tem 5 bandas e 4 bandas de transição Como cada banda de transição é especificada por um par de frequências teremos que especificar 8 frequências para a função firpmord O código necessário para implementar este filtro é mostrado logo abaixo Código 49 Exemplo do uso da função firpmord para o filtro ilustrado na figura 1 2 3 4 5 6 7 8 9 10 11 12 13 Fs 40000 k1 005 atenuacao nas bandas de rejeicao k2 001 ripple nas bandas de passagem n fo ao w firpmord5000 5500 9500 10000 15000 15500 17500 18000 0 1 0 1 0 k1 k2 k1 k2 k1 Fs b firpmn fo ao w calcula o eixo da freq espectro w freqzb1512 freq wpi freq freqFs2 plota subplot211 plotfreq 20log10absespectro titleescala dB grid on subplot212 plotfreq absespectro titleescala linear grid on Parte 4 Transformada Z em Matlab O objetivo desta prática é aprender as funções do Matlab que podem ser utilizadas para analisar sinais e sistemas no espaço Z Vale lembrar que geralmente um sinal ou um sistema pode ser descrito como uma relação de dois polinômios definidos na forma Sz b0 b1 z1bM zM a0 a1 z1aN zN Assim os coeficientes bM expressam o numerador desta razão e aN os coeficientes do polinômio de denominador Lembrese desta convenção pois a inversão desta ordem pode originar um resultado errado Os próximos pontos descrevem os principais recursos para se trabalhar com transformadas z Ponto 1 Determinação da transformada z O Matlab não tem uma função mágica que se possa inserir uma equação de sistema ou sinal e automaticamente calcular sua transformada Z Muitas vezes é preciso trabalhar matematicamente na expressão para fazer um arranjo e poder usar as funções do Matlab O próximo exemplo ilustra isto Exemplo usando as propriedades da transformada Z e a tabelada Z determine a transformada Z do sinal xn n205n2 cosπ3n2 un2 Solução a análise do sinal revela que ele foi atrasado em 2 unidades Aplicando a tabela de propriedades ver propriedade de atraso do sinal podemos reescrever Xz z2 n05n cosπ3 n un Uma nova análise na tabela de propriedades revela que um sinal multiplicado por n como é este o caso deve ser derivado no espaço Z tal como Xz z2 z ddz 05n cosπ3 n un Aplicando a transformada Z ver tabelada Z nos termos contido entre chaves temos Z 05n cosπ3 n un 1 05 cosπ3 z1 1 2 05 cosπ3 z1 025 z2 1 025 z1 1 05 z1 025 z2 Inserindo o resultado de 22 em 21 temos Xz z1 ddz 1 025 z1 1 05 z1 025 z2 Aplicando a derivada e rearranjando os termos temos finalmente Xz 025 z3 05 z4 00625 z5 1 z1 075 z2 025 z3 00625 z4 Agora que calculamos a transformada do sinal e a colocamos no formato desejado ilustrado em 20 usaremos o Matlab da seguinte forma 1 clear clc 2 teste 1 usando na forma de eq diferencias 3 n 07 4 sinal n212n2cospin230 0 1 1 1 1 1 1 5 6 teste 2 usando transf z 7 b0 0 0 025 05 00625 coef b do numerador 8 a1 1 075 025 00625 coef a do denominador 9 entrada 1 0 0 0 0 0 0 0 10 sinal filterb a entrada Ponto 2 trabalhando com a equação de diferenças Como visto na teoria de PDS uma equação de diferenças pode representar o formato matemático de um determinado sinal de saída ou até mesmo expressar matematicamente como um sistema age em um sinal de entrada produzindo uma saída Seu formato básico é dado por yn Σm0M bm xnm Σk1N ak ynk Σk0N ak ynk Σm0M bm xnm Assim ak são os coeficientes que multiplicam y e bm os que multiplicam x A função filter do Matlab pode ser usada para resolver numericamente equações diferença como mostra o exemplo abaixo Exemplo dada a equação diferença abaixo calcule yn yn1 09 yn2 xn a a resposta impulsiva hn do sistema entre o intervalo n20 a 100 b a saída do sistema considerando uma entrada de degrau unitário entre o intervalo n20 a 100 O sistema é estável Solução a reescrevendo a equação para o modelo padrão temos yn yn1 09 yn2 xn Que em Matlab é escrito como 1 b1 2 a1 1 09 3 x impseq0 20 120 4 n 20120 5 h filter b a x 6 subplot211 stemn h titleresposta impulsiva hn b Para achar a saída façamos 1 x stepseq0 20 120 2 y filter b a x 3 subplot212 stemn y titleresposta ao degrau Para verificar se o sistema é estável usaremos o critério de que uma entrada limitada ou somável como é este o caso deve também produzir uma saída limitada somável Assim basta fazer sumabsy que resultará em um valor próximo de 148 concluindo que o sistema é estável Outra forma de fazer esta análise é aplicando a transformada Z na equação do sistema que será visto mais a frente como Z ynxn Hz 1 1 1 z1 09 z2 Analisando as raízes de denominador da relação anterior que chamaremos mais à frente de polos vemos pelo código abaixo que os polos que são valores complexos conjugados tem valores positivos em sua parte real e menores que 1 o que indica que o sistema é estável quando se aplica esta determinada entrada à equação deste sistema 1 z rootsa 2 z 3 05000 08062i 05000 08062i Ponto 3 Decomposição resíduos de uma função Para estimar a função inversa de uma transformada Z geralmente desejase fracionar a razão de polinômios que é representada na forma de 20 Para isto considere o exemplo abaixo Exemplo Considere a função abaixo Faça sua decomposição em resíduos da forma mais simples possível Xz z 3 z2 4 z 1 Solução para proceder com a solução do sistema é necessário reescrever a função de tal forma que só tenhamos expoentes de z negativo Assim fazemos Xz z 3 z2 4 z 1 z2 z2 z1 3 4 z1 z2 Usando o Matlab temos 1 b 0 1 2 a 3 4 1 3 R p C residuezba 4 R 5 05000 6 05000 7 p 8 10000 9 03333 10 C 11 Este resultado quer dizer que Xz z 3 z2 4 z 1 05 1 10 z1 05 1 0333 z1 A forma simplificada ou residual permite uma fácil comparação com a tabelada Z e assim a consequente transformação do sinal do espaço Z para sua forma original No código anterior R indica os resíduos coeficiente do numerador p indica os pólos raízes do polinômio do denominador e C as constantes que devem ser somadas aos termos encontrados É também possível pegar uma expressão já na forma residual e convertêla na forma racional como mostra o código abaixo Que equivale a forma original mostrada no início desta solução Xz 0 03334z¹10 1333z¹ 03334z² Ponto 4 Computação da transformada inversa usando frações parciais É também comum a necessidade de se calcular a transformada inversa através da análise de frações parciais conforme ilustra o exemplo abaixo Exemplo dada a função abaixo determine a transformada inversa da equação abaixo Xz 1 109z¹²109z¹ ROC z09 Solução para iniciar a solução vamos usar a função poly para calcular o polinômio do denominador Em seguida usando o que foi aprendido no ponto anterior vamos usar a função residuez para decompor a expressão polinomial Veja o código Isto gera a função Xz 0251 09z¹ 05 1 09z¹² 025 1 09z¹ 0251 09z¹ 05z09 091 09z¹² 0251 09z¹ Que pela análise da tabelada Z ver anexo podemos inferir que xn vale xn 02509ⁿun 59n 109ⁿ¹un 1 02509ⁿun Ponto 5 Computação da resposta em frequência fase e plano Z de uma função A variável z é complexa e por isto possui um módulo e uma fase Como é sabido uma função racional no formato da equação 20 possui polos p e zeros z em um padrão similar a Xz z z0z z1z zM z p0z p1z pN É portanto importante conhecer o comportamento da localização dos polos e zeros além de como a magnitude e fase de z podem afetar o comportamento de Xz O exemplo abaixo ilustra isto Exemplo Para o sistema yn09yn1xn determine a encontre a transformada z de yn e desenhe seus polos e zeros b plote a magnitude e fase de para quando zejw ou seja r1 lembrando que zrejw c determine a resposta impulsiva hn do sistema a partir de sua equação de diferenças Solução a a transformada Z pode ser calculada como yn 09yn 1 xn Hz 1 1 09z¹ Para plotar polos contidos no vetor a e zeros contidos no vetor b temos b Para responder a pergunta b usaremos a função freqz tal como No código anterior a constante 100 linha 3 indica o número de componentes de frequências analisadas c Para calcular hn basta encontrar a transformada inversa de Hz que conhecemos Analisando a tabelada Z podemos chegar ao resultado Hz 1 1 09z¹ inversa hn 09ⁿun Ponto 6 plotando a superfície tridimensional de uma função transferência no espaço Z Às vezes por questões didáticas é interessante a visualização da transformada Z no espaço tridimensional Para fazer isto utilize o código abaixo Universidade Federal de Uberlândia Faculdade de Engenharia Elétrica Prof Alan Petrônio Pinheiro wwwalanengbr alaneletricaufubr Curso de EngEletrônica e Telecomunicações campus Patos de Minas Disciplina de Processamento Digital de Sinais Versão documento 13 19 Anexos da prática A Tabela de algumas funções da transformada Z B Tabela de propriedades da transformada Z Parte 5 Projeto de filtros IIR Nesta prática pretendese demonstrar o projeto de filtros IIR a partir de transformações de filtros analógicos Este projeto é realizado por etapas sequenciais Deste modo os pontos subsequentes representam as etapas necessárias para um projeto de um filtro IIR e sua implementação em cada um dos pontos abordados aqui Destacase que este material é complementar as aulas teóricas e por isto não tem como objetivo ensinar todo o conteúdo de filtros IIR Seu propósito é reforçar os conteúdos ensinados em sala e exemplificar como projetar um filtro em Matlab usando funções específicas Ponto 1 Implementação de um filtro analógico passabaixas Butterworth Este ponto se propõem a descrever o projeto básico de um filtro Butterworth analógico passabaixas usando o Matlab Um filtro analógico é especificado no domínio de Laplace especificado pela variável complexa s O primeiro passo é descobrir a ordem de um filtro que deve ser usado para satisfazer a determinados requerimentos Estas especificações são novamente mostradas no gráfico abaixo Um filtro analógico é especificado no domínio de Laplace representado pela variável complexa s e há várias maneiras de especificar um filtro Geralmente elas são feitas por duas formas i especificando a frequência de corte e ordem do filtro ou ii especificando a banda de transição ganho na banda de passagem e atenuação na banda de rejeição A primeira forma deve ser feita usando o código 51 onde se projeta um filtro com frequência corte 05 lembrando que a máxima frequência é 1 que equivale a Nyquest de ordem 3 Note a utilização da função butter do Matlab O termo low indica um filtro passabaixas e o termo s indica que o filtro é feito no domínio do analógico de Laplace Além de low podemos usar high e stop sendo este um rejeitabanda 3 ba butterNOmegaClows O resultado é b01250 e a10 10 050 01250 que matematicamente equivale a 01250 10s3 10s2 050s 01250 A segunda forma de especificação deve ser determinada a banda de transição o ganho na banda de passagem e atenuação na banda de rejeição Este processo de projeto é mostrado no Código 52 Neste caso usamos a biblioteca empregada nesta disciplina e a partir destas especificações determinase a ordem do filtro e sua função matemática Neste exemplo usamos uma banda de transição de 02π a 03π lembrando que 1π é a frequência de Nyquest neste tipo de representação O ripple da banda de passagem é de 7dB e a atenuação na banda de rejeição é 16dB neste exemplo Código 52 Projeto de um filtro Butterworth passabaixas analógico a partir da banda de transição ripple e atenuação 1 wp 02pi 0 a 02pi banda passagem 2 ws 03pi 03pi a 1pi banda rejeição 3 rp 7 ripple banda passagem em dB 4 as 16 atenuação mínima na banda rejeição em dB 5 ripple 10rp20 descobre valor do ripple sem escala log 6 atenuacao 10as20 descobre valor da atenuação sem escala log 7 b a afdbuttwp ws rp as 8 w 0pi1000pi determina um eixo para freq 9 h freqsbaw calcula a resposta em freq do filtro 10 plotwpiabsh O resultado é b 01238 e a 10 09969 04969 01238 que matematicamente equivale a 01238 10s3 09969s2 04969s 01238 Ponto 2 Transformações de frequência de analógico s para digital z Há várias formas de converter um filtro analógico especificado matematicamente no domínio s para um filtro no formato digital que deve ser especificado matematicamente no domínio z Uma destas técnicas de conversão é a i transformada de invariância ao impulso que será ilustrada no código 53 e outra é a ii transformada bilinear A transformação por invariância ao impulso utiliza a função impinvar do Matlab é ilustrada no código 53 Consideramos neste exemplo T1 lembrando que wΩT Código 53 Exemplo de conversão de uma equação no domínio s para o z usando o método de invariância ao impulso 1 PARTE 1 especificações projeto 2 wp 02pi 0 a 02pi banda passagem 3 ws 03pi 03pi a 1pi banda rejeição 4 rp 7 ripple banda passagem em dB 5 as 16 atenuação mínima na banda rejeição em dB 6 T 1 taxa de amostragem da frequencia 7 PARTE 2 desenha filtro analógico 8 b a afdbuttwp ws rp as 9 w 0pi1000pi determina um eixo para freq 10 11 h freqsbaw calcula a resposta em freq do filtro 13 subplot211 plotwpiabsh titleresp freq analogico 14 15 PARTE 3 digitaliza filtro usando inv ao impulso 16 bz az impinvarbaT 17 hz freqzbzazw 18 subplot212 plotwpiabshz titleresp freq digital O resultado é bz 00 00438 00314 e az 10 20233 14675 03690 que matematicamente equivale a transformação 01238 10s3 09969s2 04969s 01238 00438z1 00314z2 10 20233z1 14675z2 03690z3 Já a transformação bilinear emprega a função do matlab bilinear e para ver seu uso use o código 53 substituindo na linha 16 o nome impinvar por bilinear Ambas as funções servem a um mesmo propósito Contudo a transformação bilinear é geralmente mais usada por geralmente apresentar menos distorções na transformação do plano s para o círculo z Ponto 3 Projeto de filtros passaaltas rejeitabanda e passabanda Para o projeto de filtros passaaltas rejeitabandas e passabanda Para projetar estes filtros usaremos a função butter do Matlab do seguinte modo b a butterN wn high projeta um filtro passaaltas com frequência de corte wn wn wc 3dB b a butterN wn projeta um filtro passabandas de ordem 2N se wn é um vetor tal que wnw1 w2 compreende a faixa de frequências que deve passar b a butterN wn stop projeta um filtro rejeitabandas de ordem 2N se wn é um vetor tal que wnw1 w2 e compreenda as faixas de frequências que devem ser rejeitadas Em todos os casos anteriores a frequência é dada em função de π indicando que π referese a máxima frequência do sinal Destacase também que as frequências especificadas anteriormente devem ter amplitudes próximas de 3dB 07 unidades Uma outra função importante é a N wnbuttordwp ws rp as onde especificamos a banda de transição especificadas pelas frequências wp ws o ripple da banda de passagem rp em dB e a atenuação nas banda de rejeição as em dB O resultado da função é a determinação da ordem N do filtro e da frequência de corte wn No caso desta função para um filtro passabandas wp e ws são vetores com dois elementos tal que wpwp1 wp2 e wsws1 ws2 Para entender melhor estes parâmetros veja as figuras abaixo É importante destacar que esta função considera frequências normalizadas na faixa 0 1 onde 1 indica a maior frequência Código 54 Código de um filtro passabaixas de 0 a 3400Hz com transição de 3400 a 4000Hz 1 Fs 40000 2 wp 3400Fs2 3 ws 4000Fs2 4 n wnbuttordwp ws 1 20 5 b a butternwn 6 numfreq 512 7 freqzb a numfreq Fs Código 55 Código de um filtro passaaltas de 4000Hz a 20000Hz sendo esta última a máxima frequência 1 Fs 40000 2 wp 4000Fs2 3 ws 3400Fs2 4 n wnbuttordwp ws 1 20 5 b a butternwnhigh 6 numfreq 512 7 freqzb a numfreq Fs Código 56 Código de um filtro passafaixas de 4000 a 8000Hz com transições de largura 600Hz 1 Fs 40000 2 ws1 3400Fs2 3 wp1 4000Fs2 4 wp2 8000Fs2 5 ws2 8600Fs2 6 n wnbuttordwp1 wp2 ws1 ws2 1 20 7 b a butternwn 8 numfreq 512 9 freqzb a numfreq Fs Código 57 Código de um filtro rejeitafaixas que exclui as faixas de 4500 a 8000Hz e largura de transição de 500Hz 1 Fs 40000 2 wp1 4000Fs2 3 ws1 4500Fs2 4 wp2 8500Fs2 5 ws2 8000Fs2 6 n wnbuttordwp1 wp2 ws1 ws2 1 20 7 b a butternwnstop 8 numfreq 512 9 freqzb a numfreq Fs Ponto 4 Implementação do filtro Para executar uma filtragem basta ter em mãos os coeficientes designados por b e a e usar a função filter conforme ilustra o código 58 na linha 15 Código 58 Execução de um filtragem IIR 1 Fs 40000 2 wp1 4000Fs2 3 ws1 3400Fs2 4 wp2 6000Fs2 5 ws2 6600Fs2 6 n wnbuttordwp1 wp2 ws1 ws2 1 20 7 b a butternwn 8 cria sinal sintético 9 t 01Fs002 10 N maxt1Fs 11 n 0N 12 entradadiscretizada sin2pi900nFs sin2pi14500nFs sin2pi5000nFs sin2pi9000nFs 13 aplica filtro 14 sinalfiltrado filterb a entradadiscretizada 15 calcula espectros e plota sinais 16 fftsinalentrada absfftentradadiscretizadaN 17 fftsinalfiltrado absfftsinalfiltradoN 18 fresol FsN 19 f nfresol 20 subplot121 plot entradadiscretizada hold on plotsinalfiltrador 21 subplot122 plotf1N2 fftsinalentrada1N2 hold on 22 plotf1N2fftsinalfiltrado1N2r 7 freqzb a numfreq Fs Código 55 Código de um filtro passaaltas de 4000Hz a 20000Hz sendo esta última a máxima frequência 1 Fs 40000 2 wp 4000Fs2 3 ws 3400Fs2 4 n wnbuttordwp ws 1 20 5 b a butternwnhigh 6 numfreq 512 7 freqzb a numfreq Fs Código 56 Código de um filtro passafaixas de 4000 a 8000Hz com transições de largura 600Hz 1 Fs 40000 2 ws1 3400Fs2 3 wp1 4000Fs2 4 wp2 8000Fs2 5 ws2 8600Fs2 6 n wnbuttordwp1 wp2 ws1 ws2 1 20 7 b a butternwn 8 numfreq 512 9 freqzb a numfreq Fs Código 57 Código de um filtro rejeitafaixas que exclui as faixas de 4500 a 8000Hz e largura de transição de 500Hz 1 Fs 40000 2 wp1 4000Fs2 3 ws1 4500Fs2 4 wp2 8500Fs2 5 ws2 8000Fs2 6 n wnbuttordwp1 wp2 ws1 ws2 1 20 7 b a butternwnstop 8 numfreq 512 9 freqzb a numfreq Fs Ponto 4 Implementação do filtro Para executar uma filtragem basta ter em mãos os coeficientes designados por b e a e usar a função filter conforme ilustra o código 58 na linha 15 Código 58 Execução de um filtragem IIR 1 Fs 40000 2 wp1 4000Fs2 3 ws1 3400Fs2 4 wp2 6000Fs2 5 ws2 6600Fs2 6 n wnbuttordwp1 wp2 ws1 ws2 1 20 7 b a butternwn 8 cria sinal sintético 9 t 01Fs002 10 N maxt1Fs 11 n 0N 12 entradadiscretizada sin2pi900nFs sin2pi14500nFs sin2pi5000nFs sin2pi9000nFs 13 aplica filtro 14 sinalfiltrado filterb a entradadiscretizada 15 calcula espectros e plota sinais 16 fftsinalentrada absfftentradadiscretizadaN 17 fftsinalfiltrado absfftsinalfiltradoN 18 fresol FsN 19 f nfresol 20 subplot121 plot entradadiscretizada hold on plotsinalfiltrador 21 subplot122 plotf1N2 fftsinalentrada1N2 hold on 22 plotf1N2fftsinalfiltrado1N2r