PSD29007-Engtelecom(2020-2) - Prof. Marcos Moecke

Unidade 1 - Introdução

Unidade 1 - Introdução

Aula 1 (10 nov) APRESENTAÇÃO DA DISCIPLINA Nesta página será atualizado a cada aula o REGISTRO DIÁRIO E AVALIAÇÕES. A PÁGINA DA DISCIPLINA contem os materiais gerais de referência e consulta. Para a realização e entrega das atividades será utilizada a PLATAFORMA MOODLE. Chave para auto inscrição (engtelecom2020-2) Durante a pandemia do COVID19 os encontros síncronos serão realizados através da plataforma Google Meet. Os encontros (havendo a concordância de todos) serão gravados, e ocorrerão nos horários normais das aulas. As aulas terão duração de 1h20 minutos com possibilidade de ainda acrescer 20 minutos previstos para esclarecimento de dúvidas. As gravações ficarão disponíveis por 30 dias no Drive do aluno, mas não poderão ser baixadas por limitação do sistema. Além dos horários de aula síncrona, serão agendados horários de ATENDIMENTO EXTRACLASSE para que os alunos possam tirar dúvidas da disciplina ou das ferramentas de ensino. Para a comunicação entre professor-aluno, além dos avisos do SIGAA, utilizaremos a plataforma SLACK. Sugere-se que os alunos também a utilizem para comunicação entre eles, principalmente nos trabalhos em equipe. Nas aulas de PSD utilizaremos desde o inicio o Matlab e no final o Simulink e HDL Coder. Para facilitar o acesso foi disponibilizada uma VM com essa plataforma para vocês. No inicio vocês podem usar normalmente a nuvem do IFSC, mas quando forem usar o Simulink, a VM tem melhor desempenho, então recomendo que baixem a VM e instalem em seu computador. A página da wiki com as orientações para instalação e uso é Teste de VM para uso pelos alunos no ensino remoto#VM Matlab Aula 2 a 4 (12, 17, 19 nov) Transformadas de Fourier Aula 5 e 6 (24, 26 nov) Revisão de Sinais e Sistemas no tempo discreto em Matlab: Explorar a interface do Matlab. Funções de visualização das variáveis no workspace. Execução de instruções passo a passo. Escrita de script .m Uso da execução das seções de um script. Incremento de valor e execução. EXEMPLOS: Leia com atenção e execute o exemplo (Moving-Avarage Filter) na página de help da função filter. Revisão de Sinais e Sistemas no tempo discreto em Matlab: Leia com atenção o help Using FFT, abra o script clicando no botão [Open this Example]. Execute o script seção após seção. Note o uso da fft para determinar a frequência das manchas solares. Para melhorar o desempenho no Matlab recomendo que leiam a pagina do Help, . Também gostaria de lembra-los que a tecla F9 executa o código destacado no Help. Programação com scripts .m. Leia sobre manchas solares para entender o que são os dados do segundo exemplo. Sinais no dominio do tempo e dominio da frequencia. Uso da função fft Exemplo de uso da FFT %% Signal in Time Domain % Use Fourier transforms to find the frequency components of a signal buried in noise. % Specify the parameters of a signal with a sampling frequency of 1 kHz and a signal duration of 1.5 seconds Fs = 1000; % Sampling frequency T = 1/Fs; % Sampling period L = 1500; % Length of signal t = (0:L-1)*T; % Time vector % Form a signal containing a 50 Hz sinusoid of amplitude 0.7 and a 120 Hz sinusoid of amplitude 1. S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); % Corrupt the signal with zero-mean white noise with a variance of 4. X = S + 2*randn(size(t)); % Plot the noisy signal in the time domain. It is difficult to identify the frequency components by looking at the signal X(t). subplot(311); plot(1000*t(1:200),X(1:200), 'b') title('Signal Corrupted with Zero-Mean Random Noise') xlabel('t (milliseconds)') ylabel('X(t)') hold on plot(1000*t(1:200),S(1:200),'r') hold off % Signal in Frequency Domain % Compute the Fourier transform of the signal. Y = fft(X); % Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L. P2 = abs(Y/L); P1 = P2(1:L/2+1); P1(2:end-1) = 2*P1(2:end-1); A2 = angle(Y); A1 = A2(1:L/2+1); % Define the frequency domain f and plot the single-sided amplitude spectrum P1. % The amplitudes are not exactly at 0.7 and 1, as expected, because of the added noise. % On average, longer signals produce better frequency approximations. f = Fs*(0:(L/2))/L; subplot(3,1,2); plot(f,P1, 'b') title('Single-Sided Amplitude Spectrum of X(t)') xlabel('f (Hz)') ylabel('|P1(f)|') hold on % Now, take the Fourier transform of the original, uncorrupted signal and retrieve the exact amplitudes, 0.7 and 1.0. Y = fft(S); P2 = abs(Y/L); P1 = P2(1:L/2+1); P1(2:end-1) = 2*P1(2:end-1); plot(f,P1, 'r') title('Single-Sided Amplitude Spectrum of S(t)') xlabel('f (Hz)') ylabel('|P1(f)|') hold off % Angulo / fase subplot(3,1,3); plot(f,A1) %ylim([0 1.05]) title('Single-Sided Phase Spectrum of X(t)') xlabel('f (Hz)') ylabel('phase(f)') Amostragem de Sinais (Experimento 1.2) Relembrar teorema da amostragem. Efeito da amostragem abaixo da frequência de Nyquist. Aliasing. Notar que as amostras de um sinal $s_{1}(t)=cos(2\pi \times 3t)$ (3 Hz) e um sinal $s_{2}(t)=cos(2\pi \times 7t)$ (7 Hz) são idênticas quando amostrado com um sinal de 10 Hz. Experimento 1.2 % Exemplos e Experimentos baseados no livro: % DINIZ, P. S. R., DA SILVA, E. A. B., e LIMA NETTO, S. Processamento Digital de Sinais: Projeto e Análise de Sistemas. 2. ed. Porto Alegre: Bookman, 2014. 976 p. ISBN 978-8582601235. %% Experimento 1.2 fs = 10; % frequencia (Hz) de amostragem dos sinais Ts = 1/fs; fase = 0; time = 0:Ts:(1-Ts); f1 = 3; % frequencia (Hz) do sinal s_1 f2 = 7; % frequencia (Hz) do sinal s_2 s_1 = cos(2*pi*f1*time+fase); s_2 = cos(2*pi*f2*time+fase); fsa = 1000; % frequência auxiliar de amostragem usada apenas para representação dos sinais originais Tsa = 1/fsa; time_aux = 0:Tsa:(1-Tsa); figure(1); stem(time,s_1,'ob'); hold on; plot(time_aux, cos(2*pi*f1*time_aux+fase),'--k'); stem(time,s_2,'+r'); plot(time_aux, cos(2*pi*f2*time_aux+fase),'--m'); hold off; legend('s_1 discreto','s_1 contínuo','s_2 discreto','s_2 contínuo') DICAS: No help on-line da Mathworks, usando o botão [Try This Example > Try in your browser], permite executar o código no próprio browser sem ter nenhuma instalação do Matlab. Para verificar que o código realmente é executado mude a amplitude do ruído randômico para 0.1 ou 0.5, insira o comando close all antes da primeira linha, e execute todo o código [Run All] No help do Matlab, usando o botão [Open this Example], é possível executar o código seção a seção. Revisão de Sinais e Sistemas no tempo discreto em Matlab: Filtragem de Sinais (Experimentos 1.3, 2.1 e 2.2) Consulte a documentação do Matlab sobre roots, poly, linspace, logspace, residue, residuez, pretty, latex, freqs, freqz, syms, symfun, zplane. Ver também o Publish para a geração automática de relatórios em html, doc, pdf, latex ou ppt. Ver também Publishing MATLAB Code. Ver pag. 138 a 141 de Variação do Experimento 2.2 % Exemplos e Experimentos baseados no livro: % DINIZ, P. S. R., DA SILVA, E. A. B., e LIMA NETTO, S. Processamento Digital de Sinais: Projeto e Análise de Sistemas. 2. ed. Porto Alegre: Bookman, 2014. 976 p. ISBN 978-8582601235. %% Experimento 2.2 % Resposta em frequencia usando a função freqz N = 1; num = [1 0 0 0]; den = poly([0.8 0.2]) %den = [1 0.6 -0.16]; % modo 1 %[H,w]=freqz(num,den,[0:pi/100:N*pi-pi/100]); %plot(w/pi, abs(H)); % modo 2 %[H,w]=freqz(num,den); %plot(w/pi, abs(H)); % modo 3 %[H,w]=freqz(num, den, 'whole'); %plot(w/pi, abs(H)); % modo 4 freqz(num, den, 'whole'); figure(2); zplane(num,den); %% Resposta em frequencia substituindo z -> e^(jw) syms z Hf(z) = symfun(z^2/(z-0.2)/(z+0.8),z); pretty(Hf) latex(Hf) N = 1; w = [0:pi/100:N*pi-pi/100]; plot(w/pi,abs(Hf(exp(1i*w)))) %title(['$' latex(Hf) '$'],'interpreter','latex') text(0.2,2,['H(z) = ' '$$' latex(Hf) '$$'],'interpreter','latex') xlabel(['w/' '$$' '\pi' '$$'],'interpreter','latex') Verifique a diferença entre os tipos de plots comentados no código. substitua o denominador de H(z) por dois polos em [-0.8 -0.8]. verifique o que ocorre se forem utilizados polos complexos conjugados [0.3-0.4i 0.3+0.4i 0.1] verifique o que ocorre se forem utilizados polos complexos não conjugados [0.3-0.4i 0.3+0.8i] verifique o que ocorre se os polos estiverem fora do circulo unitário [1.2 -0.2]. Interprete este resultado Consulte a documentação do Matlab sobre zeros, ones, plot, stem, subplot, filter. Para usar melhor a interface do Matlab leia também Execução de seções e variação de valores nos scripts, e ainda Uso de gráficos no Matlab. Ver pag. 65 a 71 de Ver também PDF Documentation for MATLAB. Principalmente MATLAB Primer. Aula 7 (1 dez) A filtragem de sinais digitais pode ser realizada de diferentes formas: convolução (y = conv(x,h)), onde x(n) é o sinal de entrada e h(n) é a resposta ao impulso do filtro (sistema linear invariante no tempo), filtragem no domínio do tempo (y = a1.x(n)+ a2.x(n-1)+ .. ak.x(n-k)); no domínio da frequência (y = ifft(fft(x)fft(h)) Variação do Experimento 3.1 %% Variação do Experimento 3.1 do livro: % DINIZ, P. S. R., DA SILVA, E. A. B., e LIMA NETTO, S. Processamento Digital de Sinais: Projeto e Análise de Sistemas. 2. ed. Porto Alegre: Bookman, 2014. 976 p. ISBN 978-8582601235. % FILE: Ex3_1.m % Exemplificando as possiveis formas de realizar a filtragem de um sinal x(n) clc; clear all; close all; %% Definindo valores iniciais Nh = 10; Nx = 20; %Nh = 400; Nx = 10000; x = ones(1,Nx); % A resposta ao impulso de um sistema h(n) % no filtro FIR aos coeficientes b(n) = h(n) h = [1:Nh]; b = h; %% Filtrando o sinal e medindo tempos % OPÇÃO 1 - Filtragem utilizando a convolução % NOTE: length(y) = length(x) + length(h) -1 tic; % iniciar a contagem do tempo y1 = conv(x,h); t(1) = toc; % terminar a contagem e mostrar tempo no console % OPÇÃO 2 - filtragem utilizando a equação recursiva % NOTE: length(y) = length(x) tic; y2 = filter(b,1,x); t(2) = toc; % OPÇÃO 3 - filtragem utilizando a equação recursiva % aumentando o tamanho de x para que length(y3) = length(y1) x3 = [x zeros(1,length(h)-1)]; tic; y3 = filter(h,1,x3); t(3) = toc; length_y = length(x) + length(h) - 1; % OPÇÃO 4 - filtragem utilizando a FFT % a y = IFFT(FFT(x)*FFT(h)) tic; X = fft(x,length_y); H = fft(h,length_y); Y4 = X.*H; y4 = ifft(Y4); t(4) = toc; % OPÇÃO 5 - filtragem utilizando a função fftfilt % a y = IFFT(FFT(x)*FFT(h)) tic y5 = fftfilt(h,x3); t(5) = toc; disp('Comprimento do vetor de saída length(y)') disp([' ' num2str([length(y1) length(y2) length(y3) length(y4) length(y5)])]) disp('Tempo usado na filtragem em micro segundos') disp([' ' num2str(t*1e6) ' us']) %% Plotando o gráfico subplot(411);stem(y1); hold on; stem(y2,'xr'); stem(y3,'+m'); legend('y1', 'y2', 'y3') hold off subplot(412);stem(y1, 'ob');legend('y1') subplot(413);stem(y2, 'xr'); hold on; stem(zeros(size(y1)),'.w');hold off; legend('y2') subplot(414);stem(y3, '+m');legend('y3') Verificar as funções tic e toc Notar a diferença de tempo de processamento entre os processos de filtragem. A situação pode ser muito diferente conforme muda o tamanho do sinal e ordem do filtro (h(n)). Modifique a resposta ao impulso e o sinal de entrada modificando os valores das variáveis de tamanho: Nh = 10, 100, 1000; Nx = 20, 1000, 10000; Em função do sistema operacional e reserva de memória para as variáveis é importante desprezar a primeira medida de tempo. Realize 3 medidas de tempo para cada uma das 5 opções de filtragem, com pelo menos duas combinações de comprimento Nh e Nx. Verifique o tempo de processamento usando a instrução profile profile on Após ativar o profiler, execute o programa e veja o tempo total do script e de cada função chamada. profile viewer Execute no Matlab o código abaixo, e analise os 3 filtros implementados através dos seus zeros e polos. Busque tirar conclusões sobre a influência da posição dos polos e zeros (ver o gráfico do plano z) e correlacione com a resposta de frequência em magnitude (gráfico do freqz). Variação do Experimento 2.3 %% Experimento 2.3 - Filtros Digitais % Exemplos e Experimentos baseados no livro: % DINIZ, P. S. R., DA SILVA, E. A. B., e LIMA NETTO, S. Processamento Digital de Sinais: Projeto e Análise de Sistemas. 2. ed. Porto Alegre: Bookman, 2014. 976 p. ISBN 978-8582601235. % FILE: Exp2_3.m %% 1º filtro p1 = 0.9*exp(1j*pi/4); Z = [1 -1 ]'; P = [p1 p1']'; [num,den] = zp2tf(Z,P,1); [h,w] = freqz(num,den); figure(1); plot(w,abs(h)/max(abs(h))); figure(2); zplane(num,den); %% 2º filtro z1 = exp(1j*pi/8); z2 = exp(1j*3*pi/8); p1 = 0.9*exp(1j*pi/4); Z = [1 -1 z1 z1' z2 z2']'; P = [p1 p1' p1 p1' p1 p1']'; [num,den] = zp2tf(Z,P,1); [h,w] = freqz(num,den); figure(1); plot(w,abs(h)/max(abs(h))); figure(2); zplane(num,den); %% 3º filtro z1 = exp(1j*pi/8); z2 = exp(1j*3*pi/8); p1 = 0.99*exp(1j*pi/4); p2 = 0.9*exp(1j*pi/4 - 1j*pi/30); p3 = 0.9*exp(1j*pi/4 + 1j*pi/30); Z = [1 -1 z1 z1' z2 z2']'; P = [p1 p1' p2 p2' p3 p3']'; [num,den] = zp2tf(Z,P,1); [h,w] = freqz(num,den); figure(1); plot(w,abs(h)/max(abs(h))); figure(2); zplane(num,den);

Unidade 2 - Filtros IIR

Unidade 2 - Filtros IIR

Aula 8 (3 dez) Conceitos Gerais sobre Filtros Analógicos Função de transferência $H(s)={\frac {c_{0}+c_{1}s+c_{2}s^{2}+...+c_{m}s^{m}}{d_{0}+d_{1}s+d_{2}s^{2}+...+d_{n}s^{n}}},m\leq n$ Resposta em frequência: para obter a resposta em frequência é necessário avaliar $H(j\omega )=H(s)\left|{\begin{matrix}\\s=j\omega \end{matrix}}\right.$ $H(j\omega )=\left|H(j\omega )\right|e^{j\phi (\omega )}$ $\left|H(j\omega )\right|^{2}=H(j\omega )H(-j\omega )$ $e^{j2\phi (\omega )}={\frac {H(j\omega )}{H(-j\omega )}}$ O projeto de filtros analógicos é realizado em 2 etapas: projeto de um filtro passa baixas (LP) protótipo normalizado $H(p)$ com frequência de passagem $\Omega _{s}=1$ transformação em frequência para o tipo de filtro (LP, HP, BP ou BS) $H(s)=H(p)\left|{\begin{matrix}\\p=g(s)\end{matrix}}\right.$ Na sequência será mostrado como inicialmente projetar o filtro LP protótipo, e depois as transformações em frequência. No entanto, antes de projetar filtros, vejamos a análise básica de filtros analógicos utilizando o Matlab. Dado um sistema linear invariante no tempo, representado pela função de transferência $H(s)$, obter a resposta de frequência do sistema (Magnitude e Fase). $H(s)={\frac {s+1}{s^{2}+s+5}}$ $H(j\omega )={\frac {s+1}{s^{2}+s+5}}\left|{\begin{matrix}\\s=j\omega \end{matrix}}\right.$ $H(j\omega )={\frac {1+w\,\mathrm {i} }{-w^{2}+w\,\mathrm {i} +5}}$ %%Definição do filtro % Definindo os coeficientes do filtro b = [1 1]; % Numerador a = [1 1 5]; % Denominador % Calculando os zeros (raízes do numerador) e pólos (raízes do denominador) % Método 1 - usando a função tf2zp [z1,p1,k]=tf2zp(b,a) % Método 2 - obtendo as raízes z2 = roots(b); p2 = roots(a); zplane(b,a); %% Obtendo a resposta em frequência % substitituindo a variável complexa s por jw usando a função freqz freqs(b,a); % Usando cálculo simbólico e plotando o gráfico com semilogx syms s w H(s) = (s+1)/(s^2 + s + 5); pretty(H(1j*w)) latex(H(1j*w)) %% ws = logspace(-2, 1, 1000); h = H(1j*ws); subplot(211) plot(ws,abs(h)); grid on; %semilogx(ws,abs(h)); grid on; subplot(212) plot(ws,angle(h)/pi*180); grid on; %semilogx(ws,angle(h)/pi*180); grid on; Para aproximação de magnitude de filtros analógicos o projeto pode usar as aproximações de Butterworth, Chebyshev (tipo 1 ou 2) ou Cauer, mostradas na figura abaixo. Aula 9 e 10 (8 e 10 dez) Projeto de filtros analógicos do tipo Butterworth Os projetos de filtro Butterworth com função de transferência $H_{n}(s)=1/B_{n}(s)$ utilizam os polinômios de Butterworth $B_{n}(s)$ mostrados na tabela a seguir: n Fatores Polinomiais de $B_{n}(s)$ 1 $(s+1)$ 2 $(s^{2}+1.4142s+1)$ 3 $(s+1)(s^{2}+s+1)$ 4 $(s^{2}+0.7654s+1)(s^{2}+1.8478s+1)$ 5 $(s+1)(s^{2}+0.6180s+1)(s^{2}+1.6180s+1)$ 6 $(s^{2}+0.5176s+1)(s^{2}+1.4142s+1)(s^{2}+1.9319s+1)$ 7 $(s+1)(s^{2}+0.4450s+1)(s^{2}+1.2470s+1)(s^{2}+1.8019s+1)$ 8 $(s^{2}+0.3902s+1)(s^{2}+1.1111s+1)(s^{2}+1.6629s+1)(s^{2}+1.9616s+1)$ 9 $(s+1)(s^{2}+0.3473s+1)(s^{2}+s+1)(s^{2}+1.5321s+1)(s^{2}+1.879s+1)$ 10 $(s^{2}+0.3129s+1)(s^{2}+0.9080s+1)(s^{2}+1.4142s+1)(s^{2}+1.7820s+1)(s^{2}+1.9754s+1)$ Proposta de exercício Use os polinômios de Butterworth com ordens de 1 a 10 mostrados na tabela abaixo para obter os filtros . Escolha uma ordem n (entre 5 e 10) Plote a resposta em frequência em escala log da amplitude (em dB) e da fase (em rad/pi). Qual é o ganho do filtro na banda passante? Qual é a frequência de corte (-3dB) do filtro. Qual é o salto de de fase que ocorre em algumas frequências? Qual é o fator de atenuação em dB/decada após a frequência de corte? Faça o diagrama de polos e zeros desse filtro. Procure observar o que ocorre com a posição dos polos do filtro. Calcule o valor do módulo dos pólos. Projeto de filtros analógicos LP protótipo Projeto de filtros analógicos passa baixas (low pass - LP) do tipo Butterworth, considerando: $\omega _{p}$ é a frequência de passagem, $A_{p}=3dB$ é a atenuação em dB na frequência de passagem, $\omega _{s}$ é a frequência de stopband, $A_{s}$ é a atenuação em dB na frequência de stopband. Escalando as frequências em relação a ${\omega _{p}}$, teremos que $\Omega _{s}={\frac {\omega _{s}}{\omega _{p}}}$, e $\Omega _{p}={\frac {\omega _{p}}{\omega _{p}}}=1$ são as frequências de passagem e stopband do filtro protótipo $H(p)$, que tem ganho unitário e frequência de passagem 1. Casos em que o ganho na banda de passagem é $A_{p}=3dB$ Considerando o caso de filtro Butterworth com frequência de passagem ${\Omega _{p}}$ e frequência de stopband (rejeição) de $\Omega _{s}$, com ganho unitário em ${\Omega =0}$ Se considerarmos o caso particular em que na frequência de passagem o ganho (em escala linear) deve ser $G_{p}=1/{\sqrt {2}}=0,707$, que corresponde a um ganho (em escala log) $G_{p}=-3dB$, ou atenuação $A_{p}=3dB$. Obtemos o fator $\epsilon ^{2}={10^{0.1A_{p}}-1}$, ou $\epsilon ={\sqrt {10^{0.1A_{p}}-1}}$, para $A_{p}=3dB$ temos que $\epsilon =1$. Esse fator no caso dos filtros com essa atenuação $A_{p}=3dB$ acaba desaparecendo das equação de projeto. Para atenuações diferentes de 3 dB, ele ajusta a magnitude dos pólos, e afeta a ordem do filtro. Os passos para projetar um filtro analógico Hs(s) são: 0) fazer a normalização da frequência e do ganho. ${\Omega _{p}}=1$, e $\Omega _{s}={\frac {\omega _{s}}{\omega _{p}}}$ para o caso de filtros LP. $A_{p}=-(G_{p}-G_{0})$, $A_{s}=-(G_{s}-G_{0})$. 1) determinar a ordem $n$ do filtro utilizando a equação: $n\geq {\frac {\log(10^{0.1A_{s}}-1)}{2\log \Omega _{s}}}$ 2) e em seguida obter os polos do filtro: $p_{k}=e^{\left[j{\frac {(2k+n-1)}{2n}}\pi \right]},k=1,2,3,...n$ 3) Com os pólos $p_{k}$ botem-se o denominador $D(p)=\prod _{k=1}^{n}\left(p-p_{k}\right)$ da função de transferência do filtro. 4) E assim obtém-se a função de transferência do filtro protótipo $H(p)={\frac {1}{D(p)}}$ 5) Para obter a função de transferência do filtro analógico LP é necessário fazer uma transformação de frequência $H(p)->H(s)$ $H(s)=H(p)\left|{\begin{matrix}\\p={\frac {s}{\omega _{p}}}\end{matrix}}\right.$ 6) Se o ganho na frequência ${\omega =0}$ não for unitário G0, é ainda necessário ajustar o ganho do filtro do fator de ganho. Considerando que o valor do Ganho G0 seja dado em dB, teremos que $20*log10(G_{0(linear)})=G_{0(dB)}$, ou seja $G_{0(linear)}=10^{G_{0(dB)}/20}$ Exemplo Filtro LP Butterworth Projete um filtro Butterworth LP com ganho em $\omega =0$ $G_{topo}=5dB$, frequência de passagem ${\omega _{p}}=1000rad/s$ com ganho no mínimo de $G_{p}=2dB$, frequência de rejeição de $\omega _{s}=5000rad/s$, na qual o ganho deve ser inferior a $G_{s}=-25dB$ dB. Dados de $Hs(s)$, passa-baixas (lowpass-LP) $\omega _{p}=$ $\omega _{s}=$ $G_{topo(dB)}=$ $G_{p(dB)}=$ $G_{s(dB)}=$ Especificações de $Hp(p)$, passa-baixas (lowpass-LP) protótipo $\Omega _{p}=$ $\Omega _{s}=$ $A_{0(dB)}=$ $A_{p(dB)}=$ $A_{s(dB)}=$ Determinação de $Hp(p)$ $\epsilon =$ $n=$ $p_{k}=$ $D(p)=$ $Hp(p)=$ Determinação de $Hs1(s)$ substituindo $p={\frac {s}{\omega _{p}}}$ e corrigindo o ganho em $G_{topo}$ $Hs2(s)=Hs1(s)\cdot 10^{G_{topo}/20}$ Obtida a função de transferência $Hs(s)=N(s)/D(s)$ obtenha a resposta em frequência, substituindo $s=j*\omega$ Obtenha a resposta em frequência, para $\omega =0,\omega _{p},\omega _{s},10\omega _{p}$ Plote o gráfico de $|Hp(p)|$ e $|Hs(s)|$, indicando a máscara de especificação do filtro. Aula 11 (15 dez) Casos em que o ganho na banda de passagem é um valor $A_{p}$ qualquer Teremos $\epsilon ={\sqrt {10^{0.1A_{p}}-1}}$, ou $\epsilon ^{2}=(10^{0.1A_{p}}-1)$ Para projetar o filtro é necessário: 1) determinar a ordem $n$ do filtro: $n\geq {\frac {\log[(10^{0.1A_{s}}-1)/\epsilon ^{2}]}{2\log \Omega _{s}}}$ 2) obter os polos do filtro: $p_{k}=\epsilon ^{(-1/n)}e^{\left[j{\frac {(2k+n-1)}{2n}}\pi \right]},k=1,2,3,...n$ 3) obter a função de transferência: $H(p)={\frac {k}{D(p)}}$, onde $k=\prod _{k=1}^{n}\left(-p_{k}\right)=\epsilon ^{-1}$ e $D(p)=\prod _{k=1}^{n}\left(p-p_{k}\right)$. NOTA: o valor $k$ também pode ser obtido a partir de ${D(p)}$, pois corresponde ao último termo do polinômio ${D(end)}$. 4) No caso de um filtro LP é necessário ainda obter a função de transferência do filtro especificado fazendo a transformação de frequência $H(p)->H(s)$ $H(s)=H(p)\left|{\begin{matrix}\\p={\frac {s}{\omega _{p}}}\end{matrix}}\right.$ Exemplo Filtro LP Butterworth Projete um filtro Butterworth LP com ganho em $\omega =0$ $G_{topo}=0dB$, frequência de passagem ${\omega _{p}}=1000rad/s$ com atenuação máxima de $A_{p}=0.5dB$, frequência de rejeição de $\omega _{s}=5000rad/s$ com atenuação mínima de $A_{s}=30dB$. Dados de $Hs(s)$, passa-baixas (lowpass-LP) $\omega _{p}=$ $\omega _{s}=$ $G_{topo(dB)}=$ $G_{p(dB)}=$ $G_{s(dB)}=$ Especificações de $Hp(p)$, passa-baixas (lowpass-LP) protótipo $\Omega _{p}=$ $\Omega _{s}=$ $G_{0(dB)}=$ $A_{p(dB)}=$ $A_{s(dB)}=$ Determinação de $Hp(p)$ $\epsilon =$ $n=$ $p_{k}=$ $D(p)=$ $Hp(p)=$ Determinação de $Hs1(s)$ substituindo $p={\frac {s}{\omega _{p}}}$ e corrigindo o ganho em $G_{topo}$ $Hs2(s)=Hs1(s)\cdot 10^{G_{topo}/20}$ Obtida a função de transferência $Hs(s)=N(s)/D(s)$ obtenha a resposta em frequência, substituindo $s=j*\omega$ Obtenha a resposta em frequência, para $\omega =0,\omega _{p},\omega _{s},10\omega _{p}$ Plote o gráfico de $|Hp(p)|$ e $|Hs(s)|$, indicando a máscara de especificação do filtro. Ver Fazer mascara de um filtro no Matlab Ver Uso do calculo simbólico na Matlab Ver pag. 186 a 204 de Aula 12 (17 dez) Polinômios Chebyshev de primeira ordem Para o projeto dos filtros do tipo Chebyshev, são utilizados os polinômios de Chebyshev de primeira ordem, os quais são definidos pela equação trigonométrica: $C_{n}(\Omega )=cos(ncos^{-1}(\Omega ))$ Os dois primeiros polinômios são facilmente calculados como: $C_{0}(\Omega )=cos(0\times cos^{-1}(\Omega ))=cos(0)=1$ $C_{1}(\Omega )=cos(1\times cos^{-1}(\Omega ))=\Omega$ O cálculo dos demais termos pode ser feita pela relação recursiva: $C_{n}(\Omega )=2\times \Omega \times C_{n-1}(\Omega )-C_{n-2}(\Omega )$ Portanto o polinômio de grau 2 pode ser obtido por $C_{2}(\Omega )=2\times \Omega \times C_{1}(\Omega )-C_{0}(\Omega )$ $C_{2}(\Omega )=2\times \Omega \times \Omega -1$ $C_{2}(\Omega )=2\Omega ^{2}-1$ Assim os primeiros nove polinômios de Chebyshev de primeira ordem podem ser obtidos: $C_{0}(\Omega )=1\,$ $C_{1}(\Omega )=\Omega \,$ $C_{2}(\Omega )=2\Omega ^{2}-1\,$ $C_{3}(\Omega )=4\Omega ^{3}-3\Omega \,$ $C_{4}(\Omega )=8\Omega ^{4}-8\Omega ^{2}+1\,$ $C_{5}(\Omega )=16\Omega ^{5}-20\Omega ^{3}+5\Omega \,$ $C_{6}(\Omega )=32\Omega ^{6}-48\Omega ^{4}+18\Omega ^{2}-1\,$ $C_{7}(\Omega )=64\Omega ^{7}-112\Omega ^{5}+56\Omega ^{3}-7\Omega \,$ $C_{8}(\Omega )=128\Omega ^{8}-256\Omega ^{6}+160\Omega ^{4}-32\Omega ^{2}+1\,$ Esses polinômios mostram um comportamento oscilatório entre $-1<\Omega <1$. Figura 2.1 - Os primeiros cinco polinômios de Chebyshev de primeira ordem no domínio $-1<\Omega <1$ Fonte: Polinômios de Tchebychev, Wikipedia. Ver mais em Polinômio de Chebyshev Esses polinômios de Chebyshev permitem aproximar a resposta de magnitude ideal de um filtro a medida que a ordem é aumentada. Na banda de passagem o polinômio mantem uma oscilação entre os valores de mínimos e máximos. Depois da frequência de passagem, a resposta de magnitude decai monotonicamente a medida que a frequência aumenta. Por isso os filtros baseados nos polinômios de Chebyshev são também conhecidos como aproximação minimax ou equiriple. Para obter a magnitude quadrada de um protótipo, a função da magnitude possui o polinômio de Chebyshev elevado ao quadrado no denominador: $\left|H(\Omega )\right|^{2}={\frac {H_{0}^{2}}{1+\epsilon ^{2}C_{n}^{2}(\Omega )}}$ O gráfico a seguir mostra a magnitude $10log_{10}\left|H(\Omega )\right|^{2}$ em função de $\Omega$ e $n$ Figura 2.2 - Resposta de magnitude de filtro protótipo Chebyshev Fonte: Elaborado pelo autor. Projeto de filtro protótipo LP do tipo Chebyshev I Determine a ordem mínima necessária considerando: $\omega _{p}$ é a frequência de passagem do filtro LP, $A_{p}$ é a atenuação em dB na frequência de passagem, $\omega _{s}$ é a frequência de stopband do filtro, $A_{s}$ é a atenuação em dB na frequência de stopband, $\Omega _{s}={\frac {\omega _{s}}{\omega _{p}}}$, $\Omega _{p}={\frac {\omega _{p}}{\omega _{p}}}=1$ são as frequências de passagem e stopband do filtro protótipo. $n\geq {\frac {\cosh ^{-1}{\sqrt {(10^{0.1A_{s}}-1)/\epsilon ^{2}}}}{\cosh ^{-1}\Omega _{s}}}$ onde $\epsilon ^{2}=(10^{0.1A_{p}}-1)$ ou $\epsilon ={\sqrt {10^{0.1A_{p}}-1}}$. Obtenha os polos do filtro: $p_{k}=-\sinh(\varphi _{2})\sin(\theta _{k})+j\cosh(\varphi _{2})\cos(\theta _{k})\ \ \ \ \ k=1,2,3,...n$, onde $\theta _{k}=\left({\frac {(2k-1)\pi }{2n}}\right)$ e $\varphi _{2}={\frac {1}{n}}\sinh ^{-1}\left({\frac {1}{\epsilon }}\right)$ Para obter a função de transferência: $H(p)={\frac {H_{0}}{D(p)}}$, onde $D(p)=\prod _{k-1}^{n}\left(p-p_{k}\right)$ onde $H_{0}=H(0)\times d_{0}=\left\{{\begin{matrix}d_{0}&{\text{para}}&n&{\text{impar}}\\{\frac {d_{0}}{\sqrt[{}]{1+\epsilon ^{2}}}}&{\text{para}}&n&{\text{par}}\end{matrix}}\right.$ $d_{0}=\prod _{k-1}^{n}\left(-p_{k}\right)$ é o último termo do denominador $D(p)=d_{n}p^{n}+d_{n-1}p^{n-1}+\cdots +d_{1}p+d_{0}$ $\left|H(0)\right|^{2}=\left\{{\begin{matrix}1&{\text{para}}&n&{\text{impar}}\\{\frac {1}{1+\epsilon ^{2}}}&{\text{para}}&n&{\text{par}}\end{matrix}}\right.$ $\left|H(1)\right|^{2}={\frac {1}{1+\epsilon ^{2}}}$, onde $20log10(\left|H(1)\right|)=A_{p}(dB)$ Aula 13 (22 dez) Projeto de Filtros Analógicos do tipo LP, HP, BP, BS Para o projeto de filtros analógicos é necessário fazer as transformações de frequência indicadas abaixo, as quais devem ser consideradas no momento da determinação dos parâmetros do filtro protótipo LP. Transformações de frequência de filtros analógicos passa-baixas ($\Omega _{p}=1$) -> passa-baixas ($\omega _{p}$) Cálculo do protótipo com $\Omega _{s}={\frac {\omega _{s}}{\omega _{p}}}$ Substituição de variáveis $p={\frac {s}{\omega _{p}}}$ passa-baixas ($\Omega _{p}=1$) -> passa-altas ($\omega _{p}$) Cálculo do protótipo com $\Omega _{s}={\frac {\omega _{p}}{\omega _{s}}}$ Substituição de variáveis $p={\frac {\omega _{p}}{s}}$ passa-baixas ($\Omega _{p}=1$) -> passa-faixa ($\omega _{0}$ e $B$) Cálculo do protótipo com $\Omega _{s}={\Bigg |}{\frac {-\omega _{s}^{2}+\omega _{0}^{2}}{B\ \omega _{s}}}{\Bigg |}$ Substituição de variáveis $p={\frac {s^{2}+\omega _{0}^{2}}{B\ s}}$ onde $B=\omega _{p2}-\omega _{p1}$ e $\omega _{0}={\sqrt {\omega _{p2}\omega _{p1}}}$ passa-baixas ($\Omega _{p}=1$) -> rejeita-faixa ($\omega _{0}$ e $B$) Cálculo do protótipo com $\Omega _{s}={\Bigg |}{\frac {B\ \omega _{s}}{-\omega _{s}^{2}+\omega _{0}^{2}}}{\Bigg |}$ Substituição de variáveis $p={\frac {B\ s}{s^{2}+\omega _{0}^{2}}}$ onde $B=\omega _{p2}-\omega _{p1}$ e $\omega _{0}={\sqrt {\omega _{p2}\omega _{p1}}}$ Exemplo Filtro HP Chebyshev I Projete um filtro Chebyshev I HP com ganho em $\omega =0$ $G_{topo}=5dB$, frequência de rejeição ${\omega _{s}}=2500rad/s$ com atenuação mínima em relação ao topo de $A_{s}=30dB$, frequência de passagem de $\omega _{p}=12500rad/s$ com atenuação máxima em relação ao topo de $A_{p}=0.5dB$. Dados de $Hs(s)$, passa-altas (highpass-HP) $\omega _{p}=$ $\omega _{s}=$ $G_{topo(dB)}=$ $G_{p(dB)}=$ $G_{s(dB)}=$ Especificações de $Hp(p)$, passa-baixas (lowpass-LP) protótipo $\Omega _{p}=$ $\Omega _{s}=$ $A_{0(dB)}=$ $A_{p(dB)}=$ $A_{s(dB)}=$ Determinação de $Hp(p)$ $\epsilon =$ $n=$ $p_{k}=$ $D(p)=$ $d_{0}=$ $H_{0}=$ $Hp(p)=$ Determinação de $Hs1(s)$ HP substituindo $p={\frac {\omega _{p}}{s}}$ e corrigindo o ganho em $G_{topo}$ $Hs2(s)=Hs1(s)\cdot 10^{G_{topo}/20}$ Obtida a função de transferência $Hs(s)=N(s)/D(s)$ obtenha a resposta em frequência, substituindo $s=j*\omega$ Obtenha a resposta em frequência, para $\omega =0,\omega _{p},\omega _{s},10\omega _{p}$ Plote o gráfico de $|Hp(p)|$ e $|Hs(s)|$, indicando a máscara de especificação do filtro. Recesso de Natal, Ano Novo e férias Retorno as aulas em 02/02/2021 Aula 14 (2 fev) Uso das funções buttord, butter, cheb1ord, cheby1, cheb2ord, cheby2, ellipord, ellip para o projeto de filtros analógicos com Matlab (é necessário usar o parâmetro 's'). [n,Wn] = buttord(Wp,Ws,Rp,Rs,'s') Encontra a menor ordem n é a frequência de corte Wn para um filtro Butterworth analógico, com ripple na banda passante (bandpass) menor que Rp dB e atenuação na banda de rejeição (stopband) pelo menos de Rs dB. As frequências angulares Wp e Ws são dadas em rad/s. Se Wp e Ws são escalares, o filtro será um LP ou HP. Se Wp e Ws forem vetores, o filtro será um BP ou BS. [b,a] = butter(n,Wn,ftype,'s') Projeta o filtro Butterworth analógico LP, HP, BP ou BS com frequência angular de corte Wn, dependendo do valor de ftype e do número de elementos de Wn. %% Projeto de filtro passa-alta (HP) usando funções do Matlab Wp = 150; % rad/s Ws = 40; % rad/s Rp = 3; % dB Rs = 60; % dB [n,Wn] = buttord(Wp,Ws,Rp,Rs,'s') [b,a] = butter(n,Wn,'high','s'); [h,w] = freqs(b,a,logspace(1,3,1000)); semilogx(w,20*log10(abs(h)));grid on; hold on; plot([Wp Wn Ws],[-Rp -3 -Rs],'x'); hold off; title(sprintf('Filtro HP Butterworth, n = %d',n)) %% Projeto de filtro passa-faixa (BP) usando funções do Matlab Wp = [100 200]; % rad/s Ws = [50 250]; % rad/s Rp = 3; % dB Rs = 40; % dB [n,Wn] = buttord(Wp,Ws,Rp,Rs,'s'); [b,a] = butter(n,Wn,'s'); freqs(b,a,logspace(1,3,1000)) title(sprintf('Filtro BP Butterworth, n = %d',n)) [n,Wp] = cheb1ord(Wp,Ws,Rp,Rs,'s') designs a lowpass, highpass, bandpass, or bandstop analog Chebyshev Type I filter with cutoff angular frequencies Wp. [b,a] = cheby1(n,Rp,Wp,ftype,'s') designs a lowpass, highpass, bandpass, or bandstop analog Chebyshev Type I filter with passband edge angular frequency Wp and Rp decibels of passband ripple. [n,Ws] = cheb2ord(Wp,Ws,Rp,Rs,'s') designs a lowpass, highpass, bandpass, or bandstop analog Chebyshev Type II filter with cutoff angular frequencies Ws. [b,a] = cheby2(n,Rs,Ws,ftype,'s') designs a lowpass, highpass, bandpass, or bandstop analog Chebyshev Type II filter with stopband edge angular frequency Ws and Rs decibels of stopband attenuation. [n,Wn] = ellipord(Wp,Ws,Rp,Rs,'s') finds the minimum order n and cutoff frequencies Wn for an analog elliptic filter. Specify the frequencies Wp and Ws in radians per second. The passband or the stopband can be infinite. [b,a] = ellip(n,Rp,Rs,Wn,ftype, 's') designs a lowpass, highpass, bandpass, or bandstop analog elliptic filter with passband edge angular frequency Wp, Rp decibels of passband ripple, and Rs decibels of stopband attenuation. Note: The resulting bandpass and bandstop designs are of order 2n. Note: See Limitations [1] for information about numerical issues that affect forming the transfer function. Ver também Comparison of Analog IIR Lowpass Filters. Uso das funções freqs, zplane na análise da resposta em frequência de filtros analógicos. Projeto de filtros analógicos do tipo Butterworth, Chebyshev I e II e Cauer (eliptico) usando funções do Matlab. %% Projeto de filtro passa-baixas usando funções do Matlab %% Especificações do filtro Wp =16000; Ws = 20000; Ap = 0.3; As = 20; G0= 3; % Para analisar o filtro projetado, use fvtool(b,a) para observar plano s, resposta em magnitude, fase e atraso de grupo %% Butterworth [n,Wn] = buttord(Wp, Ws, Ap, As,'s') [b,a] = butter(n,Wn, 's'); %% Chebyshev I n = cheb1ord(Wp, Ws, Ap, As,'s') [b,a] = cheby1(n,Ap, Wp, 's'); %% Chebyshev II n = cheb2ord(Wp, Ws, Ap, As,'s') [b,a] = cheby2(n,As, Ws, 's'); %% Elliptic - Cauer [n, Wn] = ellipord(Wp, Ws, Ap, As,'s') [b,a] = ellip(n,Ap,As, Wn, 's'); Ver pag. 204 a 208 de Exemplos de Filtros Analógicos: Exemplo 1: Filtro passa-baixas ($f_{p}$ = 941Hz, $f_{s}$ = 1209 Hz, $A_{p}$ = 1 dB, $A_{s}$ = 20 dB) Exemplo 2: Filtro passa-altas ($f_{p}$ = 1209 Hz, $f_{s}$ = 941Hz, $A_{p}$ = 1 dB, $A_{s}$ = 20 dB) Exemplo 3: Filtro passa-faixa ($f_{s1}$ = 770 Hz, $f_{p1}$ = 811 Hz, $f_{p2}$ = 895,5 Hz, $f_{s2}$ = 1209 Hz, $A_{p}$ = 1 dB, $A_{s}$ = 30 dB, $G_{topo}$ = 5 dB). Exemplo 4: Filtro rejeita-faixa ($f_{p1}$ = 53 Hz, $f_{s1}$ = 58 Hz, $f_{s2}$ = 62 Hz $f_{p2}$ = 67 Hz, $A_{p}$ = 2 dB, $A_{s}$ = 25 dB) NOTA: No calculo do filtro lembre-se de usar as frequências angulares para $\omega _{p}$, $\omega _{s}$, $B\omega$, $\omega _{0}$. onde $f_{p}$ ($\omega _{p}$) é a frequência de passagem em Hz (rad/s), $f_{s}$ ($\omega _{s}$) é a frequência de rejeição em Hz (rad/s), $f_{0}$ ($\omega _{0}$) é a frequência central em Hz (rad/s), $B$ ($B\omega$) é a largura de banda em Hz (rad/s). Confira os projetos dos filtros plotando as respostas em frequência dos filtros protótipo H(p) e filtro final H(s) de cada um dos exemplos. Comparar o uso das funções plot, stem, fplot, semilogx, semilogy. Estudar as funções auxiliares logspace, linspace, hold, grid, subplot, legend, xlabel, ylabel, title, text. Funções de transformação de frequência: lp2lp, lp2hp,lp2bp, lp2bs. Ver em IIR Filter Design, Special Topics in IIR Filter Design. Ver pag. 208 a 218 de Aula 15 e 16 (4 e 9 fev) Filtros Digitais: Filtros IIR: transformações do tempo contínuo no tempo discreto Transformação Invariante ao Impulso (pode ser usada apenas para filtros com forte atenuação em frequência altas, ex: passa-baixas e passa-faixa) Transformação Bilinear (pode ser usada para todos tipos de filtro) $s=2fa{\frac {z-1}{z+1}}$ Obter a especificação do filtro em angulo $\theta$ entre 0 e 1, onde 1 corresponde a metade da frequência de amostragem $(fa/2)$, ou também o angulo $\pi$ Obter o valor desse angulo predistorcido $\lambda$ para compensar a distorção na frequência causada pela transformação bilinear $\lambda =2tan\left({\frac {\theta \pi }{2}}\right)$, onde $\theta ={\frac {f}{f_{a}/2}}={\frac {\omega }{\omega _{a}/2}}$ passa-baixas ($\Omega _{p}=1$) -> passa-baixas ($\lambda _{p}$) Substituição de variáveis $p={\frac {s}{\lambda _{p}}}$ Cálculo do protótipo com $\Omega _{s}={\frac {\lambda _{s}}{\lambda _{p}}}$ passa-baixas ($\Omega _{p}=1$) -> passa-altas ($\lambda _{p}$) Substituição de variáveis $p={\frac {\lambda _{p}}{s}}$ Cálculo do protótipo com $\Omega _{s}={\frac {\lambda _{p}}{\lambda _{s}}}$ passa-baixas ($\Omega _{p}=1$) -> passa-faixa ($\lambda _{0}$ e $B$) Substituição de variáveis $p={\frac {s^{2}+\lambda _{0}^{2}}{B\ s}}$ Cálculo do protótipo com $\Omega _{s}={\Bigg |}{\frac {-\lambda _{s}^{2}+\lambda _{0}^{2}}{B\ \lambda _{s}}}{\Bigg |}$ onde $B=\lambda _{p2}-\lambda _{p1}$ e $\lambda _{0}={\sqrt {\lambda _{p2}\lambda _{p1}}}$ passa-baixas ($\Omega _{p}=1$) -> rejeita-faixa ($\lambda _{0}$ e $B$) Substituição de variáveis $p={\frac {B\ s}{s^{2}+\lambda _{0}^{2}}}$ Cálculo do protótipo com $\Omega _{s}={\Bigg |}{\frac {B\ \lambda _{s}}{-\lambda _{s}^{2}+\lambda _{0}^{2}}}{\Bigg |}$ onde $B=\lambda _{p2}-\lambda _{p1}$ e $\lambda _{0}={\sqrt {\lambda _{p2}\lambda _{p1}}}$ Ver as funções de discretização usadas no Matlab: bilinear, impinvar Realizar os projetos dos exemplos anteriores, considerando uma frequência de amostragem de 8 kHz. Ver pag. 219 a 229 de Ver pag. 403 a 415 e 434 a 435 de Ver também a ferramenta fvtool, que facilita a análise dos filtros digitais. Exemplo filtro digital BP IIR - transformação bilinear Projete um filtro digital passa-banda IIR com $G_{topo}=5dB$, frequências de passagem $f_{p1}=811Hz$ e $f_{p2}=895,5Hz$ com ganho no mínimo de $G_{p}=2dB$, frequências de rejeição de $f_{s1}=770Hz$ e $f_{s1}=1209Hz$, na qual o ganho deve ser inferior a $G_{s}=-25dB$ dB. Considere dois projetos para uma frequência de amostragem $f_{a1}=4000Hz$ e também para $f_{a2}=8000Hz$. Dados de especificação do filtro passa-banda (band-pass - BP) $f_{a}=$ $f_{s1}=$ $f_{p1}=$ $f_{p2}=$ $f_{s2}=$ $G_{topo(dB)}=$ $G_{p(dB)}=$ $G_{s(dB)}=$ Dados do filtro digital $Hz(z)$, passa-banda (band-pass - BP) $\theta _{s1}=$ $\theta _{p1}=$ $\theta _{p2}=$ $\theta _{s2}=$ $G_{topo(dB)}=$ $G_{p(dB)}=$ $G_{s(dB)}=$ Dados do filtro analógico $Hs(s)$, passa-banda (band-pass - BP) $\lambda _{s1}=$ $\lambda _{p1}=$ $\lambda _{p2}=$ $\lambda _{s2}=$ $\lambda _{0}=$ $B\lambda =$ $G_{topo(dB)}=$ $G_{p(dB)}=$ $G_{s(dB)}=$ Especificações de $Hp(p)$, passa-baixas (lowpass-LP) protótipo $\Omega _{p}=1$ $\Omega _{s}=$ $A_{0(dB)}=0$ $A_{p(dB)}=$ $A_{s(dB)}=$ Determinação de $Hp(p)$ $n=$ $p_{k}=$ $z_{k}=$ $Dp(p)=$ $Np(p)=$ $Hp(p)=$ Determinação de $Hs1(s)$ substituindo $p={\frac {s}{\omega _{p}}}$ e corrigindo o ganho em $G_{topo}$ $Hs2(s)=Hs1(s)\cdot 10^{G_{topo}/20}$ Obtida a função de transferência $Hs(s)=N(s)/D(s)$ obtenha a resposta em frequência, substituindo $s=j*\omega$ Determinação de $Hz(z)$ substituindo $s=2{\frac {z-1}{z+1}}$, onde $fa=1$ Plote o gráfico de $|Hp(p)|$, $|Hs(s)|$ e $|Hz(z)|$, indicando a máscara de especificação do filtro. Aula 17 (11 fev) Ver em IIR Filter Design Uso das funções buttord, butter, cheb1ord, cheby1, cheb2ord, cheby2, ellipord, ellip para o projeto de filtros IIR digitais (sem o parâmetro 's'). O projeto dos filtros digitais IIR baseados na transformada bilinear no Matlab é realizada em dois passos: (1) Determinação da ordem do filtro; (2) Determinação dos coeficientes do numerador $b(n)$ e denominador $a(n)$ de $H(z)$ . Outros tipos de filtros IIR: yulewalk, iirnotch, iirpeak, iircomb, maxflat, invfreqz e outros filtros de modelagem paramétrica. Ver funcionamento em Anti-Causal, Zero-Phase Filter Implementation), filtfilt, fa = 200; fN = fa/2; wo = 60/fN; bw = 10/fN; [b,a] = iirnotch(wo,bw); fvtool(b,a); syms z; N(z) = poly2sym(b,z); D(z) = poly2sym(a,z); H(z) = N(z)/D(z); pretty(vpa(H(z),3)) fa = 8000; fN = fa/2; wo = 941/fN; bw = 100/fN; [b,a] = iirpeak(wo,bw); fvtool(b,a); syms z; N(z) = poly2sym(b,z); D(z) = poly2sym(a,z); H(z) = N(z)/D(z); pretty(vpa(H(z),3)) fa = 8000; fN = fa/2; fo = 1000; bw = 20/fN; [b,a] = iircomb(fa/fo,bw,'peak'); % ou use a flag 'notch' fvtool(b,a); syms z; N(z) = poly2sym(b,z); D(z) = poly2sym(a,z); H(z) = N(z)/D(z); pretty(vpa(H(z),3))

Unidade 3 - Filtros FIR

Unidade 3 - Filtros FIR

Aula 18 (18 fev) FIR - Conceitos gerais A função de transferência de transferência de um filtro digital FIR ${\begin{aligned}H(z)&=b_{0}+b_{1}z+b_{2}z^{2}+...+b_{N}z^{N}\\&=\sum _{i=0}^{N}b_{i}\cdot z^{i}\end{aligned}}$ ${\begin{aligned}H(z^{-1})&=b_{0}+b_{1}z^{-1}+b_{2}z^{-2}+...+b_{N}z^{-N}\\&=\sum _{i=0}^{N}b_{i}\cdot z^{-i}\end{aligned}}$ filtro FIR causal de ordem n mostrado acima pode ser descrito também através da equação de diferenças: ${\begin{aligned}y[n]&=b_{0}x[n]+b_{1}x[n-1]+\cdots +b_{N}x[n-N]\\&=\sum _{i=0}^{N}b_{i}\cdot x[n-i],\end{aligned}}$ Pode-se notar que a saída $y[n]$ do filtro FIR é uma soma ponderada dos N valores mais recentes das entradas $x[n]$ A realização desse filtro pode ser feita através de algoritmos de software ou circuitos digitais usando por exemplo a estrutura: A determinação da resposta ao impulso $h[n]$ do filtro pode ser feita substituindo a entrada $x[n]$ por $\delta [n]$. O resultado é $h[n]=b[n]$, e portanto a resposta ao impulso tem duração igual ao número de coeficientes N+1 (onde N é a ordem do filtro). Esse é o motivo pelo qual o filtro tem o nome de filtro de resposta ao impulso finita (FIR - Finite Impulse Response). O filtro também recebe nomes como filtro transversal, Filtro não recursivo, filtro de média móvel, e tapped delay filter (torneira com atrasos?). A função de transferência também pode ser descrita como: ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}+...+h_{N}z^{-N}\end{aligned}}$ Algumas vantagens que os filtros FIR tem sobre os IIR: 1. É possível projetar FIR com fase linear, ou seja atraso de grupo constante. Esses filtros são desejáveis na transmissão de sinais digitais. O atraso de grupo é definido como as $\tau =-{\frac {d\theta }{dw}}$, onde $\theta$ é a resposta de fase do filtro. 2. As amostras da resposta ao impulso $h[n]$ são os coeficientes do filtros $b[n]$, e portanto não precisam ser calculadas. 3. Os FIR são sempre estáveis pois tem todos os polos na origem. Também é consequência de não ter realimentação. Por isso também não tem ciclo limite que surge nos filtros IIR como resultado da resposta ao impulso de duração infinita associada a representação dos coeficientes e dos sinais com palavras de comprimento finito de bits. 4. O efeito da representação dos coeficientes e dos sinais com palavras de comprimento finito de bits, na resposta de frequência e resposta no domínio do tempo é menor que nos IIR. Os filtro FIR podem ter fase linear ou não ter fase linear. Filtros de fase linear: simétricos e antissimétricos (Tipo 1, 2, 3 e 4) Filtros de fase não linear: são todos que não se enquadram em um dos 4 tipos acima. Os filtros de fase linear possuem algumas propriedades (respostas em frequência possíveis, distribuição dos zeros em simetria quadrantal), conforme é mostrado a seguir. Inicialmente observe em exemplos as propriedades dos FIR tipo 1, 2, 3 e 4. Observe a resposta de frequência, fase, atraso de grupo, coeficientes e a simetria dos zeros em relação ao circulo unitário no diagrama de polos e zeros dos filtros abaixo. N = 10; bi = 2*(rand(1,N)-0.5) %% Tipo I - LP, HP, BS, BP b = [bi (2*rand(1,1)-0.5) flip(bi)]; fvtool(b,1); %% Tipo II - LP, BP % tem um zero em -1 b = [bi flip(bi)]; fvtool(b,1); %% Tipo III - BP % tem um zero em 1 e -1 b = [bi 0 -flip(bi)]; fvtool(b,1); %% Tipo IV - BP, HP % tem um zero em 1 b = [bi -flip(bi)]; fvtool(b,1); Figura 1 - Propriedades do filtro FIR de fase linear (Tipo 1) Ver também fvtool, Controlling FVTool from the MATLAB Command Line e Filter Analysis using FVTool. Aula 19 (23 fev) FIR - Filtros de fase linear FIR tipo I: Considere o exemplo de um filtro simétrico de ordem par (N=6) ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}+h_{3}z^{-3}+h_{4}z^{-4}+h_{5}z^{-5}+h_{6}z^{-6}\end{aligned}}$ Se ele apresenta simetria dos coeficientes $h_{0}=h_{6}$, $h_{1}=h_{5}$ e $h_{2}=h_{4}$, temos que: ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}+h_{3}z^{-3}+h_{2}z^{-4}+h_{1}z^{-5}+h_{0}z^{-6}\\&=h_{0}(1+z^{-6})+h_{1}(z^{-1}+z^{-5})+h_{2}(z^{-2}+z^{-4})+h_{3}z^{-3}\\&=(h_{0}(z^{+3}+z^{-3})+h_{1}(z^{+2}+z^{-2})+h_{2}(z^{+1}+z^{-1})+h_{3})z^{-3}\end{aligned}}$ Para obter a resposta em frequência $H(e^{-jw})$, substituímos $z=e^{jw}$. ${\begin{aligned}H(e^{-jw})&=(h_{0}(e^{+3jw}+e^{-3jw})+h_{1}(e^{+2jw}+e^{-2jw})+h_{2}(e^{+1jw}+e^{-1jw})+h_{3})e^{-3jw}\end{aligned}}$ Aplicando a identidade $2cos(nw)=e^{+jnw}+e^{-jnw}$, obtemos que: ${\begin{aligned}H(e^{-jw})&=(h_{0}2cos(3w)+h_{1}2cos(2w)+h_{2}2cos(1w)+h_{3})e^{-3jw}\\&=H_{R}(w)e^{j\theta (w)}\end{aligned}}$ Portanto a resposta em frequência tem um resposta de magnitude com uma parte real $H_{R}(w)$, e uma fase linear igual a $\theta (w)=-3w$, e portanto o atraso de grupo é $\tau =-{\frac {d\theta }{dw}}=3$, que é a metade da ordem N do filtro. Pode-se generalizar este resultado para qualquer filtro simétrico de ordem par: ${\begin{aligned}H(e^{-jw})&=\left[2\sum _{n=1}^{N/2}h\left({\frac {N}{2}}-n\right)cos(nw)+h\left({\frac {N}{2}}\right)\right]e^{-j(N/2)w}\end{aligned}}$ Na qual se percebe que a fase linear é igual a $\theta (w)=-{\frac {N}{2}}w$, e portanto o atraso de grupo é $\tau =-{\frac {d\theta }{dw}}={\frac {N}{2}}$, metade da ordem N do filtro. FIR tipo II: Considere o exemplo de um filtro simétrico de ordem impar (N=5) ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}+h_{3}z^{-3}+h_{4}z^{-4}+h_{5}z^{-5}\end{aligned}}$ Se ele apresenta simetria dos coeficientes $h_{0}=h_{5}$, $h_{1}=h_{4}$ e $h_{2}=h_{3}$, temos que: ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}+h_{2}z^{-3}+h_{1}z^{-4}+h_{0}z^{-5}\\&=h_{0}(1+z^{-5})+h_{1}(z^{-1}+z^{-4})+h_{2}(z^{-2}+z^{-3})\\&=(h_{0}(z^{+2.5}+z^{-2.5})+h_{1}(z^{+1.5}+z^{-1.5})+h_{2}(z^{+0.5}+z^{-0.5}))z^{-2.5}\end{aligned}}$ Para obter a resposta em frequência: ${\begin{aligned}H(e^{-jw})&=(h_{0}(e^{+2.5jw}+e^{-2.5jw})+h_{1}(e^{+1.5jw}+e^{-1.5jw})+h_{2}(e^{+0.5jw}+e^{-0.5jw}))e^{-2.5jw}\end{aligned}}$ Aplicando a identidade $2cos(nw)=e^{+jnw}+e^{-jnw}$, obtemos que: ${\begin{aligned}H(e^{-jw})&=(2h_{0}cos(2.5w)+2h_{1}cos(1.5w)+2h_{2}cos(0.5w))e^{-2.5jw}\\&=H_{R}(w)e^{j\theta (w)}\end{aligned}}$ Portanto a resposta em frequência tem um resposta de magnitude com uma parte real $H_{R}(w)$, e uma fase linear igual a $\theta (w)=-2.5w$, e portanto o atraso de grupo é $\tau =-{\frac {d\theta }{dw}}=2.5$, que é a metade da ordem N do filtro. Pode-se generalizar este resultado para qualquer filtro simétrico de ordem impar: ${\begin{aligned}H(e^{-jw})&=\left[2\sum _{n=1}^{(N+1)/2}h\left({\frac {N+1}{2}}-n\right)cos\left(\left(n-{\frac {1}{2}}\right)w\right)\right]e^{-j(N/2)w}\end{aligned}}$ Na qual se percebe que a fase linear é igual a $\theta (w)=-{\frac {N}{2}}w$, e portanto o atraso de grupo é $\tau =-{\frac {d\theta }{dw}}={\frac {N}{2}}$, metade da ordem N do filtro. FIR tipo III: Considere o exemplo de um filtro antissimétrico de ordem par (N=6) ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}+h_{3}z^{-3}+h_{4}z^{-4}+h_{5}z^{-5}+h_{6}z^{-6}\end{aligned}}$ Se ele apresenta simetria dos coeficientes $h_{0}=-h_{6}$, $h_{1}=-h_{5}$, $h_{2}=-h_{4}$ e $h_{3}=0$, temos que: ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}-h_{2}z^{-4}-h_{1}z^{-5}-h_{0}z^{-6}\\&=h_{0}(1-z^{-6})+h_{1}(z^{-1}-z^{-5})+h_{2}(z^{-2}-z^{-4})\\&=(h_{0}(z^{+3}-z^{-3})+h_{1}(z^{+2}-z^{-2})+h_{2}(z^{+1}-z^{-1}))z^{-3}\end{aligned}}$ Para obter a resposta em frequência $H(e^{-jw})$, substituímos $z=e^{jw}$. ${\begin{aligned}H(e^{-jw})&=(h_{0}(e^{+3jw}-e^{-3jw})+h_{1}(e^{+2jw}-e^{-2jw})+h_{2}(e^{+1jw}-e^{-1jw}))e^{-3jw}\end{aligned}}$ Aplicando a identidade $2jsin(nw)=e^{+jnw}-e^{-jnw}$, e que $j=e^{j\pi /2}$ obtemos que: ${\begin{aligned}H(e^{-jw})&=2j(h_{0}sin(3w)+h_{1}sin(2w)+h_{2}sin(1w))e^{-3jw}\\&=2(h_{0}sin(3w)+h_{1}sin(2w)+h_{2}sin(1w))e^{-j3w}e^{j\pi /2}\\&=2(h_{0}sin(3w)+h_{1}sin(2w)+h_{2}sin(1w))e^{-j(3w-\pi /2)}\\&=H_{R}(w)e^{j(\theta (w))}\end{aligned}}$ Portanto a resposta em frequência tem um resposta de magnitude com uma parte real $H_{R}(w)$, e uma fase linear igual a $\theta (w)=-(3w-\pi /2)$, e portanto o atraso de grupo é $\tau =-{\frac {d\theta }{dw}}=3$, que é a metade da ordem N do filtro. Pode-se generalizar este resultado para qualquer filtro antissimétrico de ordem par: ${\begin{aligned}H(e^{-jw})=\left[2\sum _{n=1}^{N/2}h\left({\frac {N}{2}}-n\right)sin(nw)\right]e^{-j(Nw-\pi )/2}\end{aligned}}$ Na qual se percebe que a fase linear é igual a $\theta (w)=-{\frac {N}{2}}w+\pi /2$, e portanto o atraso de grupo é $\tau =-{\frac {d\theta }{dw}}={\frac {N}{2}}$, metade da ordem N do filtro. FIR tipo IV: Considere o exemplo de um filtro antissimétrico de ordem impar (N=5) ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}+h_{3}z^{-3}+h_{4}z^{-4}+h_{5}z^{-5}\end{aligned}}$ Se ele apresenta simetria dos coeficientes $h_{0}=-h_{5}$, $h_{1}=-h_{4}$ e $h_{2}=-h_{3}$, temos que: ${\begin{aligned}H(z^{-1})&=h_{0}+h_{1}z^{-1}+h_{2}z^{-2}-h_{2}z^{-3}-h_{1}z^{-4}-h_{0}z^{-5}\\&=h_{0}(1-z^{-5})+h_{1}(z^{-1}-z^{-4})+h_{2}(z^{-2}-z^{-3})\\&=(h_{0}(z^{+2.5}-z^{-2.5})+h_{1}(z^{+1.5}-z^{-1.5})+h_{2}(z^{+0.5}-z^{-0.5}))z^{-2.5}\end{aligned}}$ Para obter a resposta em frequência: ${\begin{aligned}H(e^{-jw})&=(h_{0}(e^{+2.5jw}-e^{-2.5jw})+h_{1}(e^{+1.5jw}-e^{-1.5jw})+h_{2}(e^{+0.5jw}-e^{-0.5jw}))e^{-2.5jw}\end{aligned}}$ E portanto ${\begin{aligned}H(e^{-jw})&=2j(h_{0}sin(2.5w)+h_{1}sin(1.5w)+h_{2}sin(0.5w))e^{-2.5jw}\\&=2(h_{0}sin(2.5w)+h_{1}sin(1.5w)+h_{2}sin(0.5w))e^{-2.5jw}e^{j\pi /2}\\&=2(h_{0}sin(2.5w)+h_{1}sin(1.5w)+h_{2}sin(0.5w))e^{-j(2.5w-\pi /2)}\\&=H_{R}(w)e^{j\theta (w)}\end{aligned}}$ Portanto a resposta em frequência tem um resposta de magnitude com uma parte real $H_{R}(w)$, e uma fase linear igual a $\theta (w)=-(2.5w-\pi /2)$, e portanto o atraso de grupo é $\tau =-{\frac {d\theta }{dw}}=2.5$, que é a metade da ordem N do filtro. Pode-se generalizar este resultado para qualquer filtro antissimétrico de ordem impar: ${\begin{aligned}H(e^{-jw})&=\left[2\sum _{n=1}^{(N+1)/2}h\left({\frac {N+1}{2}}-n\right)sin\left(\left(n-{\frac {1}{2}}\right)w\right)\right]e^{-j(Nw-\pi )/2)}\end{aligned}}$ Na qual se percebe que a fase linear é igual a $\theta (w)=-(Nw-\pi )/2$, e portanto o atraso de grupo é $\tau =-{\frac {d\theta }{dw}}={\frac {N}{2}}$, metade da ordem N do filtro. Propriedades dos filtros FIR de fase linear Como mostrado acima, os filtros que exibem simetria ou antissimetria em seus coeficientes (ou resposta ao impulso), apresentam fase linear (ou atraso de grupo constante). Também foi mostrado que o atraso de grupo é igual a N/2 onde N é a ordem do filtro. Foi demonstrado por Rabiner *** que apenas esses quatro tipos de filtro FIR possuem essa característica, portanto pode-se afirma que "Se e somente se o filtro FIR possui coeficientes simétrico ou antisimétricos ele possui fase linear". Em relação a posição dos zeros, é possível verificar que cada zero sobre o circulo unitário produz uma resposta de magnitude nula na frequencia angular correspondente e um salto de fase de $\pi$. N = 5; bi = 2*(rand(1,N)-0.5) b = [bi (2*rand(1,1)-0.5) flip(bi)]; [h,w] = freqz(b,1,'whole'); figure(1); subplot(421); plot(w/pi,20*log10(abs(h))); grid on; xlabel('\omega / \pi'); ylabel ('magnitude - dB'); subplot(423); plot(w/pi,angle(h)/pi); grid on; xlabel('\omega / \pi'); ylabel ('fase - rad / \pi'); subplot(425); plot(w/pi,unwrap(angle(h))/pi); grid on; xlabel('\omega / \pi'); ylabel ('fase - rad / \pi'); subplot(427); grpdelay(b,1); xlabel('\omega / \pi'); ylabel ('atraso de grupo - amostras'); subplot(4,2,[2,4,6,8]); zplane(b,1); xlabel('real'); ylabel ('imaginario'); Aula 20 (25 fev) Também devido a existência (ou não) de zeros em $z=1\equiv w=0$ e $z=-1\equiv w/\pi =1$, que corresponde a frequência de Nyquist $f_{N}=f_{a}/2$, mostramos que a resposta de magnitude nessas frequencias é nula (ou não). Assim os tipos 1, 2, 3 e 4 de filtros FIR resultam em: Tipo 1: permite projetar filtros LP, HP, BP e BS. Tipo 2: (tem um zero em -1) permite projetar filtros LP e BP. Tipo 3: (tem um zero em 1 e -1) permite projetar filtros BP. Tipo 4: (tem um zero em 1) permite projetar filtros HP e BP. Essa característica é importante conhecer antecipadamente pois implicará no número de coeficientes e na escolha do tipo de (anti)simetria. Por exemplo para filtro BS apenas o Tipo 1 pode ser usado. Outra propriedade a ser destacada é em relação aos zeros do filtro. Em primeiro vamos analisar a consequencia da simetria nos coeficientes $h(n)=h(N-n)$: ${\begin{aligned}H(z)&=\sum _{n=0}^{N}h(n)\cdot z^{-n}\\&=\sum _{n=0}^{N}h(N-n)\cdot z^{-n}\end{aligned}}$ fazendo na segunda equação N-n = m, temos que os limites n = 0 -> m = N, e n = N -> m = 0. ${\begin{aligned}H(z)&=\sum _{m=0}^{N}h(m)\cdot z^{-N+m}\\&=z^{-N}\sum _{m=0}^{N}h(m)\cdot z^{m}\end{aligned}}$ ${\begin{aligned}H(z)=z^{-N}H(z^{-1})\end{aligned}}$ Com a mesma análise para antissimetria nos coeficientes $h(n)=-h(N-n)$: ${\begin{aligned}H(z)=-z^{-N}H(z^{-1})\end{aligned}}$ Nessas duas equações é possível perceber que se $z=z_{-1}$ é um zero então $1/z$ também será um zero de $H(z)$. No caso de zeros reais, se temos um zero $z=a$ então $z=1/a$ também é um zero, exceto se $a=1$ ou $a=-1$. Por outro lado, se todos os coeficientes b(n) do filtro são reais, então os zeros complexos, aparecem em pares complexos conjugados $z_{2}=r_{2}\cdot e^{j\phi }$ e $z_{2}^{*}=r_{2}\cdot e^{-j\phi }$, e seus reciprocos $z_{2}=1/r_{2}\cdot e^{j\phi }$ e $z_{2}^{*}=1/r_{2}\cdot e^{-j\phi }$ também são zeros de $H(z)$. Esse tipo de disposição dos zeros denominamos de simetria quadrantal. Figura 2 - Simetria quadrantal de filtros FIR de fase linear FONTE: pag. 251 a 261 de firtype - Type of linear phase FIR filter - Mathwork Projeto de filtros FIR pelo método da série de Fourier: Usando a representação dos filtros ideais LP, HP, BP, BS, com frequências de corte $w_{ci}$ e ganho unitário na banda de passagem e ganho zero na banda de rejeição, e considerando que a magnitude das respostas em frequência é uma função periódica em $2\pi$, e conhecendo as equações de síntese e análise de um sinal (ou sistema) ${\begin{aligned}H(e^{jw})&=\sum _{n=-\infty }^{\infty }c(n)e^{-jnw}\end{aligned}}$ onde ${\begin{aligned}c(n)&={\frac {1}{2\pi }}\int _{-\pi }^{\pi }H(e^{jw})e^{jnw}\mathrm {d} w\end{aligned}}$ É possível coeficientes da série de Fourier de filtros ideais: LP, HP, BP, BS Figura 3 - Magnitude da resposta em frequência de filtros Passa-baixas (Low-pass) Figura 4 - Resposta em frequência de filtros LP de fase linear ${\begin{aligned}c_{\text{LP}}(n)&={\frac {1}{2\pi }}\int _{-\pi }^{\pi }H_{LP}(e^{jw})e^{jnw}\mathrm {d} w\\&={\frac {1}{2\pi }}\int _{-w_{c}}^{w_{c}}e^{jnw}\mathrm {d} w\\&={\frac {1}{2\pi }}\left({\frac {e^{jnw}}{jn}}\right){\Biggr |}_{-w_{c}}^{w_{c}}\\&={\frac {e^{jnw_{c}}-e^{-jnw_{c}}}{2j\pi n}}\\&={\frac {sin(nw_{c})}{\pi n}}\end{aligned}}$ ${\begin{aligned}c_{\text{LP}}(n)&={\frac {w_{c}}{\pi }}sinc(w_{c}n)\end{aligned}}$ ${\begin{aligned}c_{\text{LP}}(n)&=\left\{{\begin{matrix}{\frac {\omega _{c}}{\pi }};&\qquad n=0\\{\frac {\sin(\omega _{c}n)}{\pi n}};&\qquad \left|n\right|>0\end{matrix}}\right.\end{aligned}}$ De modo semelhante é possível obter os coeficientes dos filtros HP, BP e BS. Passa-altas (High-pass) ${\begin{aligned}c_{\text{HP}}(n)&=sinc(n)-{\frac {w_{c}}{\pi }}sinc(w_{c}n)\end{aligned}}$ ${\begin{aligned}c_{\text{HP}}(n)&=\left\{{\begin{matrix}1-{\frac {\omega _{c}}{\pi }};\qquad n=0\\-{\frac {\sin(\omega _{c}n)}{\pi n}};\qquad \left|n\right|>0\end{matrix}}\right.\end{aligned}}$ Passa-faixa (Band-pass) ${\begin{aligned}c_{\text{BP}}(n)&={\frac {1}{\pi }}(w_{c2}sinc(w_{c2}n)-w_{c1}sinc(w_{c1}n))\end{aligned}}$ ${\begin{aligned}c_{\text{BP}}(n)=\left\{{\begin{matrix}{\frac {\omega _{c2}-\omega _{c1}}{\pi }};\qquad n=0\\{\frac {\sin(\omega _{c2}n)-\sin(\omega _{c1}n)}{\pi n}};\qquad \left|n\right|>0\end{matrix}}\right.\end{aligned}}$ Rejeita-banda (Band-stop) ${\begin{aligned}c_{\text{BS}}(n)&=sinc(n)-{\frac {1}{\pi }}(w_{c2}sinc(w_{c2}n)-w_{c1}sinc(w_{c1}n))\end{aligned}}$ ${\begin{aligned}c_{\text{BS}}(n)=\left\{{\begin{matrix}1-{\frac {\omega _{c2}-\omega _{c1}}{\pi }};\qquad n=0\\-{\frac {\sin(\omega _{c2}n)-\sin(\omega _{c1}n)}{\pi n}};\qquad \left|n\right|>0\end{matrix}}\right.\end{aligned}}$ onde sabe-se que $sinc(n)=\delta (n)$, ou seja $sinc(n)=1$ para $n=0$ e $sinc(n)=0$ para $n\neq 0$. Ver pag. 249 a 256 de Nota importante sobre o uso do Matlab Diferentemente do que se espera, a função sinc(x) é diferente de sin(x)/x, sendo igual a sin(pi*x)/(pi*x). Veja o help do próprio Matlab. sinc Sin(pi*x)/(pi*x) function. sinc(X) returns a matrix whose elements are the sinc of the elements of X, i.e. y = sin(pi*x)/(pi*x) if x ~= 0 = 1 if x == 0 where x is an element of the input matrix and y is the resultant output element. % Example of a sinc function for a linearly spaced vector: t = linspace(-5,5); y = sinc(t); plot(t,y); xlabel('Time (sec)');ylabel('Amplitude'); title('Sinc Function') See also square, sin, cos, chirp, diric, gauspuls, pulstran, rectpuls, and tripuls. Reference page in Help browser doc sinc Aula 21 (2 mar) Filtros FIR com janelas fixas: O uso da janela retangular no "janelamento" dos coeficientes da série de Fourier, resulta no fenômeno de Gibbs na magnitude da resposta em frequência, conforme mostrado a seguir. ${\begin{aligned}w(n)&=1;\qquad -M\leq n\leq M\\&=0;\qquad \left|n\right|\geq M\end{aligned}}$ Aplicando a equação de síntese da série obtemos: ${\begin{aligned}\Psi (e^{jw})&=\sum _{n=-M}^{M}e^{-jnw}\\&={\frac {e^{-jMw}(e^{j(2M+1)w}-1)}{e^{jw}-1}}\\&={\frac {e^{-jMw}(e^{j(2M+1)w/2}-e^{-j(2M+1)w/2})e^{j(2M+1)w/2}}{(e^{jw/2}-e^{-jw/2})e^{jw/2}}}\\&={\frac {sin[(2M+1)w/2]}{sin(w/2)}}\end{aligned}}$ Note que esta função $\Psi (e^{jw})={\frac {sin[(2M+1)w/2]}{sin(w/2)}}$ tem um máximo $(2M+1)$ em $w=0$, e cruza o zero em $(2M+1)w/2]=\pm \pi$, portanto a lagura do lóbulo central é de $4\pi /(2M+1)$. Além disso percebe-se que se aumentamos o tamanho da janela retangular (2M+1), a largura do lóbulo central é reduzida proporcionalmente. Figura 5 - função $\Psi (e^{jw})$ Ao fazer o "janelamento" dos coeficientes da série de Fourier da resposta em frequência do filtro ideal, estamos multiplicando a série de coeficientes pelo janela retangular $h(n)=w(n)\times C_{LP}$, conforme mostra a figura a seguir. Figura 6 - Janelamento temporal (rectwin) dos coeficientes Essa multiplicação no domínio do tempo corresponde a uma convolução no domínio da frequência. ${\begin{aligned}H_{M}(e^{j\omega })&={\frac {1}{2\pi }}\int _{-\pi }^{\pi }H_{LP}(e^{j\theta })\Psi (e^{j(\omega -\theta })\mathrm {d} \theta \\&={\frac {1}{2\pi }}\int _{-\pi }^{\pi }H_{LP}(\theta )\Psi (\omega -\theta )\mathrm {d} \theta \end{aligned}}$ A qual é mostrada graficamente na figura a seguir. Figura 7 - Convolução da resposta do filtro ideal H_{LP} com a função $\Psi (e^{jw})$ Figura 8 - Aproximação da resposta de magnitude com janela retangular Para reduzir o ripple devido ao corte dos coeficientes, são usadas as funções de janelamento temporal no projeto de filtros digitais. Tipos de janelas temporais usadas no projeto de filtros digitais. Retangular $w(n)=1;\qquad -M\leq n\leq M$ Triangular $w(n)=1-{\frac {\left|n\right|}{M+1}};\qquad -M\leq n\leq M$ Bartlett $w(n)=1-{\frac {\left|n\right|}{M}};\qquad -M\leq n\leq M$ Hann $w(n)=0.5+0.5\cos \left({\frac {2\pi n}{2M+1}}\right),-M\leq n\leq M$ Hamming $w(n)=0.54+0.46\cos \left({\frac {2\pi n}{2M+1}}\right);\qquad -M\leq n\leq M$ Blackman $w(n)=0.42+0.5\cos \left({\frac {2\pi n}{2M+1}}\right)+0.08\cos \left({\frac {4\pi n}{2M+1}}\right);\qquad -M\leq n\leq M$ em todas as janelas $w\left(n\right)=0$ quando $\left|n\right|\geq M$ onde $M$ é $N/2$ para $N$ par e $(N+1)/2$ para $N$ impar Uso de janelas fixas no Matlab : rect, triang, bartlett, hann, hamming, blackman, blackmanharris, nuttall. ver também apodization function L = 64; wvtool(rectwin(L), triang(L), bartlett(L), hann(L), hamming(L), blackman(L), blackmanharris(L), nuttallwin(L)); Exercício: Efeito do uso das janelas no filtro FIR Escreva um código que implemente o cálculo dos coeficientes da série de Fourier dos filtros ideais LP, HP, BP, BS, e aplique as janelas temporal Retangular, Bartlett, Hann, Hamming e Blackman, conforme equações acima. Verifique o efeito das diferentes janelas temporais sobre a magnitude e fase da resposta em frequência, sobre a resposta ao impulso e a posição dos zeros no plano z. Use wc como a frequência de corte do filtro LP e HP, e wc1 e wc2 como frequências de corte dos filtros BS e BP Escolha uma ordem N para o filtro (valores entre 10 e 100 são normais) e depois verifique as mudanças que ocorrem com a variação de N. A resposta ao impulso h(n) dos filtros é dada pela multiplicação dos coeficientes da série de Fourier "C" pela janela temporal "w". Os coeficientes do filtro são obtidos fazendo b(n) = h(n) = C(n).*w(n). Dicas: Use multiplicação elemento a elemento do Matlab (.*), ver Array vs. Matrix Operations. O número de coeficientes sempre será igual a (N+1)=(2M+1). O ripple na banda de rejeição é sempre proporcional ao ripple na banda passante (visualizar com a escala linear de magnitude). Use fvtool para visualizar os gráficos de um filtro digital H(z). Essa ferramenta não serve para filtros analógicos pois utiliza freqz e não freqs. Note que a banda de transição aumenta a medida que o ripple diminui. Note que o aumento da ordem do filtro reduz a banda de transição, mas "quase" não afeta a amplitude do ripple. Varie o wc e veja o efeito no filtro. Aula 22 (4 mar) Use o código abaixo e verifique o efeito das diferentes janelas temporais sobre a magnitude da resposta em frequência, sobre a resposta ao impulso, posição dos zeros no plano z, etc. No código o wc é a frequência de corte do filtro LP, N é a ordem do filtro, CLP é são os coeficientes da série de Fourier do filtro LP ideal multiplicados pela janela retangular, bRET, bHAM e bBLACK são os coeficientes dos filtros usando respectivamente as janelas retangular, Hamming e Blackman. Note que: 1) O número de coeficientes sempre será igual a (N+1)=(2M+1). 2) A função ylim([-0.1 0.1]) foi usada para destacar o ripple na banda passante. 3) O ripple na banda de rejeição é sempre proporcional ao ripple na banda passante (visualizar com a escala linear de magnitude). 4) A banda de transição aumenta a medida que o ripple diminui. 5) O aumento da ordem do filtro reduz a banda de transição, mas "quase" não afeta a amplitude do ripple. N = 32; wc = 0.5; M = N/2; CLP = wc*sinc(wc*(-M:M)); bRET = CLP; bHAM = CLP.*hamming(2*M+1)'; bBLACK = CLP.*blackman(2*M+1)'; fvtool(bRET,1,bHAM,1,bBLACK,1); legend('rectwin', 'Hamming', 'Blackman'); ylim([-0.1 0.1]) Estudar no Matlab as funções wvtool, window Ver Window Designer, Design and analyze spectral windows. Ver FIR Filter Design Ver pag. 263 a 268 de Ver este artigo de 1978 que compara diversas janelas temporais, aplicadas não a filtragem mas a análise de sinais. F. J. Harris, "On the use of windows for harmonic analysis with the discrete Fourier transform," in Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978, doi: 10.1109/PROC.1978.10837. Ver este artigo de 1978 que aplica janelas aos filtros FIR L. Rabiner, B. Gold and C. McGonegal, "An approach to the approximation problem for nonrecursive digital filters," in IEEE Transactions on Audio and Electroacoustics, vol. 18, no. 2, pp. 83-106, June 1970, doi: 10.1109/TAU.1970.1162092. Projeto de filtro FIR utilizando janelas temporais fixas. Exemplos de projeto Projetar um filtro LP usando uma janela temporal fixa (verificar a janela que atende a especificação) wp = 0.2*pi; Ap = 1 dB; Gp = 0 dB ws = 0.4*pi; As = 40 dB; Para o projeto do filtro, o primeiro passo é escolher uma janela que atenda a atenuação na banda de passagem e na banda de rejeição. Em seguida é necessário determinar a ordem do filtro que atende a especificação de largura de banda de transição. Por último será necessário ajustar o valor de wc para que o filtro esteja dentro das especificações. Ao final do projeto, deverá ser informado o tipo de janela escolhida, a ordem do filtro, se é do tipo 1, 2, 3 ou 4, e o valor de wc. Aula 23 (9 mar) Projetar um filtro HP usando uma janela temporal fixa (hamming, bartlett-hanning, hanning). ws = 0.4*pi; Ap = 0.2 dB; Gp = 0 dB wp = 0.6*pi; As = 50 dB; Comparar os 3 tipos de janela, a ordem obtida, e o valor de wc em cada projeto. Para o projeto dos filtros LP FIR de janela fixa, uma possível solução para reduzir o número de passos é: PASSO 1 - Escolher o tipo de janela de acordo com a atenuação do lóbulo lateral Asl e As. PASSO 2 - Estimar a ordem N1 do filtro considerando os parâmetros Dw = |ws -wp| PASSO 3 - Calcule os coeficientes clp do filtro LP considerando N1 e wc1 = |ws + wp|/2. Calcule os valores da janela win e obtenha a resposta ao impulso do filtro h = clp * win. PASSO 4 - Verifique o valor medido de Dwm = wsm-wpm, e faça a correção da ordem do filtro em função do desvio constatado. N2 = N1*Dwm/Dw. PASSO 5 - Refaça os cálculos dos coeficientes Clp do filtro ideal, da janela e da resposta ao impulso para a nova ordem N2. Repita o PASSO 3 até 5, até obter um filtro com a menor ordem que atenda as especificações de Dw. PASSO 6 - Desloque a frequência de corte wc de modo a ter a banda de transições posicionada corretamente entre $wp$ e $ws$. $wc2=wp+(wp-wAp)$. PASSO 7 - Ajustar o ganho de topo do filtro $Glin=10^{(G_{TopoMedido}-G_{TopoEspecificado}/20)}$ Tabela - Estimativa da atenuação do lóbulo lateral $A_{sl}$ da janela, atenuação do primeiro lóbulo lateral do filtro $A_{s}$, e largura da banda de transição $\Delta \omega$, para um filtro LP FIR de janela fixa. Janela $A_{sl}$ $A_{s}$ $\Delta \omega$ Retangular 13.3 20.33 0.92$\pi$/N Triangular 26.6 27.41 Bartlett 26.5 27.48 Hann 31.5 44.03 3.11$\pi$/N Bartlett-Hanning 35.9 40.77 Hamming 42.5 54.08 3.32$\pi$/N Bohman 46.0 51.84 7.01$\pi$/N Parzen 53.1 56.89 Blackman 58.1 75.25 5.56$\pi$/N Flat Top 88.0 106.3 Blackman-Harris 92.1 108.8 Nutfall 93.8 109.7 Dados acima obtidos para um filtro passa baixas de ordem N = 64 com $\omega _{c}=0.5\pi$, mas podem ser utilizados para estimativa de atenuação e ordem em outros filtros. Ver artigos: A new window and comparison to standard windows Yeong Ho Ha ; Pearce, J.A. IEEE Transactions on Acoustics, Speech, and Signal Processing, Feb. 1989, Vol.37(2), pp.298-301. Some windows with very good sidelobe behavior Nuttall, A. IEEE Transactions on Acoustics, Speech, and Signal Processing, February 1981, Vol.29(1), pp.84-91 Projeto de filtro FIR No Matlab, além das funções que calculam as janelas, já mencionadas anteriormente, também tem a função fir1 que implementa os passos descritos anteriormente, calculando os coeficientes da série de Fourier dos filtros ideais LP, HP, BP ou BS, e aplica a janela indicada. O uso dessa função simplifica o projeto. N = <ordem> h_fir = fir1(N,Wn,hamming(N+1)); [Hw,w] =freqz(h_fir); plot(w/pi,20*log10(abs(Hw))) title(['hamming N = ' num2str(N)]) %fvtool(h_fir,1) Projetar os filtros LP, HP e BP de acordo com as especificações dadas para o projeto de filtro IIR na AE1 - Projeto de Filtro Analógico e Filtro Digital IIR, considerando uma frequência de amostragem fa > que 2 * fmax especificada. Faça seu projeto usando 3 janelas fixas diferentes. Compare os filtros FIR obtidos com os filtros IIR, considerando, resposta em frequência, fase, atraso de grupo, estabilidade, procedimento de projeto, exatidão dos cálculos x ajuste de parâmetros, duração da resposta ao impulso. Ver pag. 256 a 265 de Aula 24 (11 mar) Filtros Digitais: Filtros FIR Projeto de filtro FIR utilizando janelas temporais ajustáveis Uso de janelas ajustáveis no Matlab: kaiser, chebyshev, gauss, tukey, taylor. L = 64; r = 60; % Chebyshev e Tukey alpha = 3; % Gauss betha = 8; % Kaiser nbar = 10; % Taylor wvtool(kaiser(L,betha), chebwin(L,r), gausswin(L,alpha),tukeywin(L,r), taylorwin(L,nbar,-r)); Para a janela de Kaiser, a estimação do fator $\beta$ e da ordem do filtro $N$ são obtidos por: $\beta =\left\{{\begin{matrix}0.1102(\alpha -8.7),&\alpha >50,\\0.5842(\alpha -21)^{0.4}+0.07886(\alpha -21),&50\geq \alpha \geq 21,\\0,&\alpha <21.\end{matrix}}\right.$ $N={\frac {\alpha -8}{2.285\Delta \omega }}+1.$ onde $\alpha$ é a atenuação do lóbulo lateral e $\Delta \omega$ é a largura da banda de transição em rad/amostra. A janela de Kaiser é definida por: $w(n)={\frac {I_{0}\left(\beta {\sqrt {1-({\frac {n-\alpha }{\alpha }})^{2}}}\right)}{I_{0}(\beta )}}$ onde :$I_{0}(x)=1+\sum _{k=1}^{\infty }{\left({\frac {({\frac {x}{2}})^{k}}{k!}}\right)}^{2}$ é a função de Bessel de ordem zero [2] Utilizando o Matlab é possível estimar esses valores utilizando a função kaiserord. Exemplo da obtenção de um filtro passa baixa com $f_{pass}=1000Hz$, $f_{stop}=1500Hz$, $f_{amostragem}=8000Hz$ atenuação de 40 dB na "stopband" %% Calculo do filtro de kaiser, sem ajustes % Especificaçao fsamp = 8000; fcuts = [1000 1500]; Ap = 1; As = 40; ftype = 'low'; fN = fsamp/2; wp = fcuts(1)/fN; ws = fcuts(2)/fN; Dw = abs(ws-wp); % Calculo da janela de Kaiser beta = 0.5842*(As-21)^0.4+0.07886*(As-21); n = ceil((As-8)/(2.285*Dw*pi)+1); Wn = (wp+ws)/2; wkaiser = kaiser(n+1,beta); Forma alternativa de projeto usando a função kaiserord fsamp = 8000; fcuts = [1000 1500]; Ap = 1; As = 40; mags = [1 0]; devs = [1-10^(-Ap/20) 10^(-As/20)]; [n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp); wkaiser = kaiser(n+1,beta); h_fir = fir1(n,Wn,ftype,wkaiser,'noscale'); A partir das especificações do filtro é possível obter um projeto usando a função fir1. Essa função basicamente aplica o método da janela ao filtro ideal especificado pela(s) frequência(s) de corte $W_{n}$. h_fir = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); [Hw,w] =freqz(h_fir); plot(w*fsamp/2/pi,20*log10(abs(Hw))) title(['Kaiser filter N = ' num2str(n)]) %fvtool(h_fir,1) Como resultado do projeto a partir das equações de Kaiser é obtido o filtro abaixo: Figura 9 - Filtro LP com janela de Kaiser, sem ajustes. No entanto realizando ajustes tanto do ganho no topo $G_{topo}$, na largura da banda de transição $\Delta \omega$, e na ordem do filtro $n$, é possível reduzir essa ordem obtendo o filtro abaixo: Figura 10 - Filtro LP com janela de Kaiser, com ajustes. Ver as funções fir1, kaiserord do Matlab. Ver pag. 266 a 273 de Uso das funções window e fir1 do Matlab para projeto de filtro FIR Uso de janelas ajustáveis no Matlab: kaiser, chebyshev, gauss, tukey, taylor. L = 64; r = 60; % Chebyshev e Tukey alpha = 3; % Gauss betha = 8; % Kaiser nbar = 10; % Taylor wvtool(kaiser(L,betha), chebwin(L,r), gausswin(L,alpha),tukeywin(L,r), taylorwin(L,nbar,-r)); Figura 11 - Filtro LP com janela de Chebyshev, Taylor e Gaussiana. Aula 25 (16 mar) Filtro de Parks-McClellan. Funções firpm e firpmord, [https://ieeexplore-ieee-org.ez130.periodicos.capes.gov.br/stamp/stamp.jsp?tp=&arnumber=1083419 Chebyshev Approximation for Nonrecursive Digital Filters with Linear Phase]. Remez exchange algorithm - o básico em Wolfram Alpha Remez exchange algorithm - o detalhe com implementação em Matlab. Como resolver as anomalias na banda de transição. [3] Exemplo do projeto de um filtro passa-baixas, com minima ordem (Filtro de Parks-McClellan) com frequência de passagem de 1000 Hz e frequência de rejeição de 1500 Hz, dada uma frequência de amostragem de 8000 Hz. Considere que a atenuação na banda de rejeição é de no mínimo 40 dB e o ripple máximo na banda passante é de 1 dB. fa = 8000; Ap = 1; Ar = 40; fp = 1000; fr = 1500; f = [fp fr]; a = [1 0]; dev = [(10^(Ap/20)-1)/(10^(Ap/20)+1) 10^(-Ar/20)]; [n,fo,ao,w] = firpmord(f,a,dev,fa); b = firpm(n,fo,ao,w); [h,w] = freqz(b,1,1024,fa); plot(w, 20*log10(abs(h))); hold on; plot([0 fr fr fa/2], [Ap/2 Ap/2 -Ar -Ar],':m') plot([0 fp fp], [-Ap/2 -Ap/2 -(Ar+30)],':m'); ylim([-(Ar+30) Ap/2+10]) Exemplo do projeto de um filtro passa-faixa, com minima ordem (Filtro de Parks-McClellan) com frequências de passagem de 1500 e 1700 Hz e e frequências de rejeição de 1000 e 3000 Hz, dada uma frequência de amostragem de 8000 Hz. Considere que a atenuação na primeira banda de rejeição é de no mínimo 40 dB e 60 dB na segunda banda de rejeição. O ripple máximo na banda passante é de 1 dB e fa = 8000; Uso do [4] para projeto de filtro IIR, FIR equiripple e FIR com janela. Uso do [5] no Simulink. Exemplo do projeto de um filtro passa-baixas, com minima ordem (Filtro de Parks-McClellan) com frequência de passagem de 1000 Hz e frequência de rejeição de 1500 Hz, dada uma frequência de amostragem de 8000 Hz. Considere que a atenuação na banda de rejeição é de no mínimo 40 dB e o ripple máximo na banda passante é de 0.4 dB. Aula 26 (18 mar) Projeto de filtro FIR projetar os filtros usando: 1) Janela fixa 2) Janela ajustável 3) Parks-McClellan. garantir que o filtro seja de menor ordem em cada caso, mas que esteja dentro das especificações. se necessário ajustar os valores de fs, fp, Ap, As, e a ordem do filtro, indicando o critério utilizado para o ajuste. %% Projetar o filtro passa baixas fp = 1200 Hz; fs = 1380 Hz; fa = 8000 Hz; Ap = 1 dB; Ar = 50 dB; %% Projetar o filtro passa altas fs = 1200 Hz; fp = 1380 Hz; fa = 8000 Hz; Ap = 1 dB; Ar = 50 dB; %% Projetar o filtro passa faixa fs1 = 800 Hz; fp1 = 900 Hz; fp2 = 1000 Hz; fs2 = 1300 Hz; fa = 8000 Hz; Ap = 1 dB; Ar = 50 dB; %% Projetar o filtro rejeita faixa fp1 = 800 Hz; fs1 = 900 Hz; fs2 = 1000 Hz; fp2 = 1300 Hz; fa = 8000 Hz; Ap = 1 dB; Ar = 50 dB;

Unidade 4 - Realização de Filtros

Unidade 4 - Realização de Filtros

Aula 27 (23 mar) Realização de Filtros FIR A função de transferência de transferência de um filtro digital FIR ${\begin{aligned}H(z)&=b_{0}+b_{1}z+b_{2}z^{2}+...+b_{N}z^{N}\\&=\sum _{i=0}^{N}b_{i}\cdot z^{i}\end{aligned}}$ ${\begin{aligned}H(z)&=b_{0}+b_{1}z^{-1}+b_{2}z^{-2}+...+b_{N}z^{-N}\end{aligned}}$ Como $H(z)=Y(z)/X(z)$ ${\begin{aligned}Y(z)&=b_{0}X(z)+b_{1}z^{-1}X(z)+b_{2}z^{-2}X(z)+...+b_{N}z^{-N}X(z)\end{aligned}}$ Conhecendo a transformada Z inversa de $X(z)=Z\{x[n]\}$, e a propriedade do atraso $X(z)z^{-k}=Z\{x[n-k]\}$, o filtro FIR causal de ordem $N$ mostrado acima pode ser descrito também através da equação de diferenças: ${\begin{aligned}y[n]&=b_{0}x[n]+b_{1}x[n-1]+\cdots +b_{N}x[n-N]\\&=\sum _{i=0}^{N}b_{i}\cdot x[n-i],\end{aligned}}$ Pode-se notar que a saída $y[n]$ do filtro FIR é uma soma ponderada dos N valores mais recentes das entradas $x[n]$ A realização desse filtro pode ser feita através de algoritmos de software ou circuitos digitais usando por exemplo alguma das estruturas mostradas a seguir. A título de exemplo vamos considerar um filtro FIR de ordem 4, e para permitir uma notação de vetores com os índices do Matlab (maiores que 0), a função de de transferência e sua equação de diferenças são mostradas a seguir: ${\begin{aligned}H(z)&={b(1)+b(2)z^{-1}+b(3)z^{-2}+b(4)z^{-3}+b(5)z^{-4}}\\Y(z)&={b(1)+b(2)z^{-1}+b(3)z^{-2}+b(4)z^{-3}+b(5)z^{-4}}X(z)\\y[n]&={b(1)x[n]+b(2)x[n-1]+b(3)x[n-2]+b(4)x[n-3]+b(5)x[n-4]}\end{aligned}}$ Realização de filtros FIR na Forma Direta Para exemplificar as diferentes realizações utilizaremos com base um filtro de ordem 4 representado pela função de transferência ${\begin{aligned}H(z)&=h(1)+h(2)z^{-1}+h(3)z^{-2}+h(4)z^{-3}+h(5)z^{-4}\\\end{aligned}}$ A implementação direta desse sistema discreto no tempo pode ser representada pelo modelo esquemático a seguir. Figura 4.1 - Realização de filtros FIR na Forma Direta Fonte: Elaborado pelo autor. Realização de filtros FIR na Forma Transposta A transposição consiste na inversão do fluxo de todos os sinais, substituição de nós de soma por derivações e as derivações por soma. A entrada e saída também devem ser invertidas. A realização da transposição não altera o sistema implementado. Figura 4.2 - Realização de filtros FIR na Forma Transposta Fonte: Elaborado pelo autor. Figura 4.3 - Realização de filtros FIR na Forma Transposta Fonte: Elaborado pelo autor. Realização de filtros FIR de fase linear Os filtros FIR de fase linear podem ser com coeficientes simétricos (tipo I e II) ou antissimétricos (tipo III e IV). Tipo I - Simétrico de ordem par ${\begin{aligned}H(z)&=h(1)+h(2)z^{-1}+h(3)z^{-2}+h(4)z^{-3}+h(5)z^{-4}\\&=h(1)(1+z^{-4})+h(2)(z^{-1}+z^{-3})+h(3)z^{-2}\end{aligned}}$ Figura 4.4 - Realização de filtros FIR de fase linear Simétrico I Fonte: Elaborado pelo autor. Tipo II - Simétrico de ordem impar ${\begin{aligned}H(z)&=h(1)+h(2)z^{-1}+h(3)z^{-2}+h(4)z^{-3}+h(5)z^{-4}+h(6)z^{-5}\\&=h(1)(1+z^{-5})+h(2)(z^{-1}+z^{-4})+h(3)(z^{-2}+z^{-3})\end{aligned}}$ Figura 4.5 - Realização de filtros FIR de fase linear Simétrico II Fonte: Elaborado pelo autor. Tipo III - Antissimétrico de ordem par ${\begin{aligned}H(z)&=h(1)+h(2)z^{-1}+h(3)z^{-2}+h(4)z^{-3}+h(5)z^{-4}\\&=h(1)(1-z^{-4})+h(2)(z^{-1}-z^{-3})+0z^{-2}\\&=h(1)(1-z^{-4})+h(2)(z^{-1}-z^{-3})\end{aligned}}$ Figura 4.6 - Realização de filtros FIR de fase linear Antisimétrico III Fonte: Elaborado pelo autor. Tipo IV - Antissimétrico de ordem impar ${\begin{aligned}H(z)&=h(1)+h(2)z^{-1}+h(3)z^{-2}+h(4)z^{-3}+h(5)z^{-4}+h(6)z^{-5}\\&=h(1)(1-z^{-5})+h(2)(z^{-1}-z^{-4})+h(3)(z^{-2}-z^{-3})\end{aligned}}$ Figura 4.7 - Realização de filtros FIR de fase linear Antisimétrico IV Fonte: Elaborado pelo autor. Uso do Matlab e Simulink Realização de Filtros FIR usando o FDATool Estudar estrutura de filtros disrcetos FIR no Matlab, Filter Realization Wizard - Reference, Filter Realization Wizard - User Guide. Ver pag. 303 a 312 de . Aula 28 (25 mar) Realização de filtros IIR de 2ª ordem: Forma Direta I e II, e Forma Transposta I e II. $H(z)={\frac {Y(z)}{X(z)}},H(z)={\frac {b_{0}z^{2}+b_{1}z^{1}+b_{2}}{z^{2}+a_{1}z^{1}+a_{2}}},H(z)={\frac {b_{0}+b_{1}z^{-1}+b_{2}z^{-2}}{1+a_{1}z^{-1}+a_{2}z^{-2}}}$ Separando H(z) em dois blocos $\ H(z)=H_{1}(z)H_{2}(z)$, e obtendo o sinal intermediário W(z) ou Y(z) dependendo da ordem dos blocos. Figura 4.8 - Separação do filtro IIR H(z) em H1(z) e H2(z) Fonte: Elaborado pelo autor. Com o ordenamento dos blocos $\ H_{1}(z)$ e $\ H_{2}(z)$ em ordem direta teremos a Forma Direta I: $H_{1}(z)={\frac {W(z)}{X(z)}}=b_{0}+b_{1}z^{-1}+b_{2}z^{-2}$ $H_{2}(z)={\frac {Y(z)}{W(z)}}={\frac {1}{1+a_{1}z^{-1}+a_{2}z^{-2}}}$ Podemos obter a realização de $\ H_{1}(z)$ na forma direta. $\ W(z)=(b_{0}+b_{1}z^{-1}+b_{2}z^{-2})X(z)$ Para obter a realização de $\ H_{2}(z)$ , é necessário reescrever a saída $\ Y(z)$ em função de $\ W(z)$ e das saídas anteriores $\ Y(z)z^{-1}$ e $\ Y(z)z^{-2}$: $\ Y(z)={\frac {W(z)}{1+a_{1}z^{-1}+a_{2}z^{-2}}}$ $\ Y(z)({1+a_{1}z^{-1}+a_{2}z^{-2}})=W(z)$ $\ Y(z)=W(z)-a_{1}Y(z)z^{-1}-a_{2}Y(z)z^{-2}$ Figura 4.9 - Realização de filtros IIR na Forma Direta I Fonte: Elaborado pelo autor. Com o ordenamento dos blocos $\ H_{2}(z)$ e $\ H_{1}(z)$ em ordem reversa teremos a Forma Direta II: $H_{2}(z)={\frac {V(z)}{X(z)}}={\frac {1}{1+a_{1}z^{-1}+a_{2}z^{-2}}}$ $\ V(z)=X(z)-a_{1}V(z)z^{-1}-a_{2}V(z)z^{-2}$ $H_{1}(z)={\frac {Y(z)}{V(z)}}=b_{0}+b_{1}z^{-1}+b_{2}z^{-2}$ $\ Y(z)=(b_{0}+b_{1}z^{-1}+b_{2}z^{-2})V(z)$ Figura 4.10 - Realização de filtros IIR na Forma Direta II Fonte: Elaborado pelo autor. Considerando que os sinais no centro são idênticos podemos simplificar e obter a Forma Direta II (Canônica): Figura 4.11 - Realização de filtros IIR na Forma Direta II Canônica Fonte: Elaborado pelo autor. Considerando as regras de transposição podemos obter a forma transposta I e II. A transposição consiste na inversão do fluxo de todos os sinais, substituição de nós de soma por derivações e as derivações por soma. A entrada e saída também devem ser invertidas. A realização da transposição não altera o sistema implementado. Figura 4.12 - Realização de filtros IIR na Forma Transposta I Fonte: Elaborado pelo autor. Figura 4.13 - Realização de filtros IIR na Forma Transposta II Fonte: Elaborado pelo autor. Realização de filtros IIR de ordem maior que 2: Forma Direta I e II, Transposta I e II, Cascata, Paralela Os filtros IIR de ordem superior a 2 podem ser implementados nas FD I ou II e na FT I ou II. No entanto nessa configuração tendem a ficar instáveis ao terem os coeficientes quantizados, e também terem uma significativa alteração da resposta em frequência. Para reduzir esses problemas uma possível solução é a decomposição em filtros de 2ª ordem para serem associados na forma em Cascata ou Paralela. Realização de Filtros usando o comando realizemdl,[6] do MatLab O que são os objetos de sistema usados pela função realizemdl? Ver What Are System Objects? e [ Fs = 30000; % Sampling Frequency Fpass = 12000; % Passband Frequency Fstop = 13000; % Stopband Frequency Dpass = 0.01; % Passband Ripple Dstop = 0.01; % Stopband Attenuation flag = 'scale'; % Sampling Flag % Calculate the order from the parameters using KAISERORD. [N,Wn,BETA,TYPE] = kaiserord([Fpass Fstop]/(Fs/2), [1 0], [Dstop Dpass]); % Calculate the coefficients using the FIR1 function. b = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag); hFIR = dsp.FIRFilter; hFIR.Numerator = b; % Para definir diretamente os coeficientes realizemdl(hFIR) % Para definir os coeficientes através de uma matriz de entrada realizemdl(hFIR,'MapCoeffsToPorts','on');

Unidade 5 - Projeto Final

Unidade 5 - Projeto Final

Uso do Simulink para processamento de sinais Aula 29 (30 mar) Uso dos blocos de simulação - Sources sinewave/DSP, chirp/DSP, random noise. Uso dos blocos de simulação - Sink time scope e spectrum analyzer. Outros blocos mux, demux, sum, product, matrix concatenate. Veja também esse exemplo com a criação de um modelo com filtro passa-baixas e passa-altas com a análise no domínio do tempo (scope) Digital Filter Design Block, e outro exemplo com análise no domínio da frequência (Spectrum Analyser) Filter Frames of a Noisy Sine Wave Signal in Simulink. Uma opção para a visualização de mais de um sinal no DF em um único [Spectrum Analyser] é utilizar o bloco [Matrix Concatenate] conforme mostrado no modelo em Display Frequency-Domain Data in Spectrum Analyzer. Aula 30 (1 abr) Para entender a diferença entre processamento por amostra e processamento por quadro, leia o texto em (Sample- and Frame-Based Concepts), ou Conceitos de Sistemas baseados em amostras e quadros. Para converter entre fluxos de dados de diferentes taxas de frames, ou entre sample-based e frame-based, os blocos buffer e unbuffer devem ser utilizados. O Simulink possui muitos blocos que permitem a implementação d filtros. Para os filtros IIR e FIR recomendo o uso dos blocos Digital Filter Design, Filter Realization Wizard. Para aprender o básico do Simulink recomendo assistir a estes vídeos da Mathworks Getting Started with Simulink -(duração 3m26s), ensina como usar a biblioteca de blocos [library browser] e o editor de modelos [model editor]. É mostrado um exemplo inserindo um gerador senoidal, um bloco de ganho e um osciloscópio. Detalhes de como conectar os blocos e uso de vetores nos parâmetros da simulação são também mostrados. Simulink for New Users -(duração 54m07s), esse vídeo ensina a usar o ambiente de simulação de forma mais efetiva, indo desde a criação de um modelo até o compartilhamento do modelo com uma equipe e criação de templates para uso como base de projetos. O uso de dashboard é mostrado de forma resumida. Como definir parâmetros de uma simulação Para definir parâmetros para o simulador, tais como frequências de amostragem, amostras por quadro, outras frequências, amplitudes, etc, o recomendado é usar o Model Explorer, [View > Model Explorer] [CTRL+H], então selecione [Model Workspace] e adicione os parâmetros e seu valor. Use o nome do parâmetro nos blocos do Simulink. Figura 5.1 - Definição de Parâmetros no Simulink Fonte: Elaborado pelo autor. Qual solver utilizar em DSP? Para escolher e configurar o tipo de solver que o simulador irá utilizar, selecione [Simulation > Model Configuration Parâmeters], e escolha o solver. Normalmente o melhor opção de solver para sistemas discretos com uma única taxa de amostragem é [Type: Fixed-setup] e [Solver: discrete] Como ativar o Profiler? Para comparar o desempenho de um sistema para processamento por amostras, e para diferentes tamanhos de quadros, é possível utilizar o Profiler. Para a medição dos tempos foi utilizado sempre o mesmo computador, e ativado o Profiler [Analysis > Performance Tools > Show Profiler Report]. Faça uma simulação de alguns segundos, e se quiser use a média de N simulações. Normalmente a primeira deve ser descartada. Como organizar o modelo? Agrupar todas as entradas em um subsistema (input_source) Agrupar as saídas (medições) em um subsistema (measuarements) Agrupar todos os blocos do sistema desenvolvido sob um subsistema (DUT). ver: Create a Subsystem, Create Subsystem from Selection, Navigate Model Hierarchies pode interessar a vocês também os vídeos: [7], [8] AE3 - Uso do Simulink para processamento de sinais Atividades Use o Simulink para construir o modelo de um sistema com filtros conforme descrito na página Digital Filter Design Block. Acrescente um bloco Matrix Concatenate para visualizar a saída no Spectrum Analyser (ver dica em Display Frequency-Domain Data in Spectrum Analyzer. Use o Model Explorer e Model Workspace, para definir no mínimo os parâmetros fsine, fa e spf, que serão utilizados para Sample time = 1/fa e Samples per frame = spf, e no bloco Sine Wave, defina Frequency (Hz) = fsine Escolha um valor de fa entre 8 kHz e 44,1 kHz, e fsine = 1 kHz, 2 kHz, 3 kHz. Defina as atenuações dos filtros HP e LP como Apass = 1 dB e Astop = 40 dB. Ative o Profiler para medir o tempo de processamento do modelo. Escolha como Solver o [Type: Fixed-setup] e [Solver: discrete] Faça a simulação por pelo menos 40 segundos Figura 5.2 - AE3 - Modelo do Sistema Fonte: Elaborado pelo autor. Figura 5.3 - AE3 - Domínio do Tempo (DT) Fonte: Elaborado pelo autor. Figura 5.4 - AE3 - Domínio da Frequência (DF) Fonte: Elaborado pelo autor. Figura 5.5 - AE3 - Domínio da Frequência (DF) com max-hold ativo Fonte: Elaborado pelo autor. Entregas O arquivo AE3.slx com o modelo completo. A impressão da imagem do modelo completo com sinais destacados por cores, e tamanhos e tipos de sinais mostrados no diagrama. A impressão da imagem da tela do Time Scope entre 20 e 20,1 segundos. A impressão da imagem da tela do Spectrum Analyser ao final da simulação. Tempo de total de processamento mostrado pelo Profiler para spf = 1, 16 e 256. Use o campo de comentário para registrar esse dado. Informe os parâmetros utilizados para spf, fa, fsine. Use o campo de comentário para registrar esse dado. Informe a ordem do filtro HP e LP que você utilizou. Use o campo de comentário para registrar esse dado. Analise os gráficos no DT, DF e explique o formato no tempo e o espectro do sinal y(n). Aula 31 (6 abr) Passos adicionais para converter de Matlab/Simulink para HDL. Ver Basic Guidelines for Modeling HDL Algorithm in Simulink Ver Guidelines for Model Setup and Checking Model Compatibility Ver Modeling with Simulink, Stateflow, and MATLAB Function Blocks Executar o hdlsetup('nome_modelo') para configurar o modelo antes da conversão. Se quiser customizar o comando edite o arquivo hdlsetup.m e salve com myhdlsetup.m edit hdlsetup.m Experimentar com esse circuito contador até M. Create Simulink Model for HDL Code Generation Figura 5.6 - Divisão HDL_DUT e testbench Fonte: Elaborado pelo autor. Figura 5.7 - Subsistema HDL_DUT Fonte: Elaborado pelo autor. NOTA: Para funcionar a simulação é importante que o solver esteja setado para discreto, pois o contador é um sistema discreto pois utiliza os blocos Delay. Antes de simular aplique o hdlsetup ao modelo. hdlsetup('nome_modelo') A invés de fixar o [Output data type] dos blocos [Constant] em uint32 conforme mostrado na página da mathworks, escolha um tipo de dado da entrada count_threshold, e deixe os blocos internos como [Inherit: Inherit via back propagation], para que eles se ajustem conforme a necessidade. Verifique os códigos VHDL gerados quando é utilizado double, int8 e uint32. Para gerar o código VHDL pode se digitado makehdl('HDL_DUT'), ou também usar a interface gráfica, com clique-direito sobre o bloco HDL_DUT e selecionando [HDL Code > Generate HDL for Subsystem], ou ainda pelo menu [Code > HDL Code > Generate HDL]. Para que o testbench seja separado em arquivos do tipo .dat, em menu [Code > HDL Workflow Advisor], no item 3.1.3 Set Testbench Option, marque [x] Multi-file test bench e [x] Use file I/O for read/write bench data. Para gerar o código VHDL para o testbench usar a interface gráfica pelo menu Code > HDL Code > Generate testbench]. Para fazer a simulação do sistema em VHDL, abra o Modelsim e vá para a pasta do projeto. cd hdl_prj/hdlsrc/hdlcoder_simple_up_counter/ Compile os arquivos do projeto do HDL_DUT_compile.do Compile os arquivos do testbench do HDL_DUT_tb_compile.do Realize a simulação automática do contador. do HDL_DUT_tb_sim.do Figura 5.8 - Subsistema HDL_DUT simulado no Simulink Fonte: Elaborado pelo autor. Figura 5.9 - Subsistema HDL_DUT simulado no Modelsim Fonte: Elaborado pelo autor. Observe que o sinal simulado em VHDL é identico ao sinal simulado pelo Simulink, além disso observe que a frase abaixo aparece no console confirmando que o sistema HDL_DUT está correto. O próximo passo seria usar o VDHL para implementar o sistema. ** Note: **************TEST COMPLETED (PASSED)************** Orientação sobre o Projeto Final Aula 32 (8 abr) Orientação sobre o projeto final. Dicas sobre uso de alguns blocos do Simulink Como atividade de estudo realizar a montagem do modelo indicado na figura abaixo e fazer a simulação, usar fa = 8000Hz. PASSO 1 - Entrar no Simulink Abra o Matlab e crie uma pasta PSD29007 Entre no Simulink criando um novo modelo digitando simulink Em seguida salve o modelo vazio com o nome aula34.slx PASSO 1 - Definir os parâmetros que serão utilizados nos blocos usando o Model Explorer fa = 8000 % para Sampling frequency, e usar 1/fa para o Sample time. spf = 1 (ou 16 ou 256 ou ...) % Samples per frame; v1 = 0.2 % Sine amplitude; f1 = 500 % Sine frequency; G = 0.05 % Gain; PASSO 2 - Construir o circuito com as fontes (source) de sinal Insira os blocos que deseja utilizar para construir a seleção da fonte (ver o diagrama de blocos) Figura D.1 - Modelagem da seleção das fontes de sinal Fonte: Elaborado pelo autor. Configure os parâmetros de cada bloco. Configurar o bloco sinewave Amplitude — Amplitude of sine waves v1 Frequency (Hz) — Frequency of each sine wave f1 Phase offset (rad) — Phase offset 0 (default) Sample mode — Continuous or discrete sampling mode Discrete (default) Output complexity — Real or complex waveform Real (default) Computation method — Method for computing discrete-time sinusoids Trigonometric fcn (default) Sample time — Sample period 1/fa Samples per frame — Samples per frame spf Resetting states when re-enabled — State behavior inside enabled subsystems Restart at time zero (default) Configurar o bloco chirp Frequency sweep — Type of frequency sweep Linear (default) Sweep mode — Sweep mode Unidirectional (default) Initial frequency (Hz) — Initial frequency 1 Target frequency (Hz) — Target frequency fa/2 Target time (s) — Target time of sweep 2 Sweep time (s) — Sweep time 2 Initial phase (rad) — Initial phase of cosine output 0 (default) Sample time — Output sample period 1/fa Samples per frame — Samples per frame spf Output data type — Output data type Double (default) Configurar o bloco random source Source type — Uniform or Gaussian Uniform (default) Minimum — Minimum value of uniform distribution -1 Maximum — Maximum value of uniform distribution 1 (default) Repeatability — Repeatability of block output Not repeatable Inherit output port attributes — Inherit output port parameters from downstream block off (default) Sample mode — Discrete or continuous Discrete (default) Sample time — Output sample period 1/fa Samples per frame — Samples per output frame spf Output data type — Output data type Double (default) Complexity — Complexity of output Real (default) PASSO 3 - Construir o circuito com os receptores (sink) de sinal Insira os blocos que deseja utilizar para receptores (ver o diagrama de blocos) Figura D.2 - Modelagem das medições de sinal Fonte: Elaborado pelo autor. Configure os parâmetros de cada bloco. Configurar o bloco time scope Abra a configuração [View > Configuration properties] Na aba [Main] Open at simulation start — Specify when scope window opens on (default for Time Scope) Display the full path — Display block path on scope title bar off (default) Number of input ports — Number of input ports on scope block 2 Layout — Number and arrangement of displays 1 coluna 2 linhas Sample time — Simulation interval between scope updates -1 (for inherited) (default) Input processing — Channel or element signal processing Columns as channels (frame based) (default for Time Scope) Maximize axes — Maximize size of plots Auto (default for Time Scope) Na aba [Time] Time span — Length of x-axis to display 2 Time span overrun action — Display data beyond visible x-axis Wrap (default) Time units — x-axis units Metric (default for Time Scope) Time display offset — x-axis offset 0 (default) | scalar | vector Time-axis labels — Display of x-axis labels All (default for Time Scope) Show time-axis label — Display or hide x-axis labels on (default for Time Scope) Na aba [Display] Active display — Selected display 1 (default) Title — Display name %<SignalLabel> (default) Show legend — Display signal legend on Show grid — Show internal grid lines on (default) Plot signals as magnitude and phase — Split display into magnitude and phase plots off (default) Y-limits (Minimum) — Minimum y-axis value -10 (default) Y-limits (Maximum) — Maximum y-axis value 10 (default) Y-label — Y-axis label Amplitude (default for Time Scope) Na aba [Logging] Limit data points to last — Limit buffered data values off and 5000 (default) Decimation — Reduce amount of scope data to display and save off, 2 (default) Log data to workspace — Save data to MATLAB workspace off (default) Variable name — Name of saved data variable ScopeData (default) Save format — MATLAB variable format Dataset (default) Configurar o bloco spectrum analyser Abra a configuração [View > Spectrum Settings] Main option Type — Type of spectrum to display Power (default) View — Spectrum view Spectrum (default) Sample rate — Sample rate of the input signal in hertz Inherited (default) Full frequency span — Use entire Nyquist frequency interval on (default) RBW (Hz) — Resolution bandwidth Auto (default) Window Options Overlap (%) — Segment overlap percentage 0 (default) Window — Windowing method Hann (default) NENBW — Normalized effective noise bandwidth Trace Options Units — Spectrum units dBm (default) Averages — Number of spectral averages 1 (default) Reference load — Reference load 1 (default) Scale — Scale of frequency axis Linear (default) Offset — Constant frequency offset 0 (default) Normal trace — Normal trace view off Max hold trace — Maximum hold trace view on Min hold trace — Minimum hold trace view off (default) Two-sided spectrum — Enable two-sided spectrum view off (default) Ver mais em: Configure Spectrum Analyzer PASSO 4 - Desenvolver o sistema entre estes blocos O ideal é antes de inserir qualquer sistema verificar se a conexão direta das fontes com os receptores está funcionando. Criar um Filtro usando o [Digital Filter Design] Use o bloco Digital Filter Design para projetar e implementar um filtro digital FIR ou IIR. De um duplo click sobre o bloco, e escolha as configurações que desejar. Após configurar o filtro ou importar os coeficientes projetados no Matlab, clique em [Design Filter]. No menu [Edit > Convert Structure], escolha a implementação desejada. Para melhorar a visualização dos sinais no Simulink, faça os seguintes procedimentos: Ative a visualização das dimensões das portas [Display > Signal & Ports > Signal Dimensions] ou [ALT]+D+S+D+[ENTER] Ative a visualização do tipo de dados das portas [Display > Signal & Ports > Port Data Types] ou [ALT]+D+S+D+D+[ENTER] Ative a visualização da cor para indicar os tempos de amostragem sinais [Display > Sample Time > Colors] ou [ALT]+D+T+C+[ENTER] Velocidade de simulação Note que se for usado um bloco To Audio Device, a velocidade de simulação será determinada pelo tempo real de amostragem do sistema. Para acelerar a simulação é necessário remover ou comentar (comment out) os dispositivos de tempo real como o To Audio Device ou From Audio Device. Conforme já mencionado a simulação usando frames no lugar de samples também acelera a simulação. Experimente utilizar spf = 1, 2, 256, 1024, e perceba a diferença de velocidade de simulação. Tipo de entrada para analisar o filtro. Experimente as entradas senoidais, chirp, random noise e mistura dessas, assim como sinais constantes, e observe os sinais de entrada e saída no DT e DF. Para medir o tempo de simulação ative o Profiler. Exemplo de modelo de um sistema e resultados da simulação Figura D.3 - Modelagem de um sistema em Simulink Fonte: Elaborado pelo autor. DICAS: Fazer a análise no DT e DF para um ruido uniforme (-1 a 1) Figura D.4 - Análise no DT com fonte de ruído randomico Fonte: Elaborado pelo autor. Figura D.5 - Análise no DF com fonte de ruído randomico Fonte: Elaborado pelo autor. Fazer a análise no DT e DF para um sinal de chirp (0 a 4000Hz) Figura D.6 - Análise no DT com fonte de chirp Fonte: Elaborado pelo autor. Figura D.7 - Análise no DF com fonte de chirp Fonte: Elaborado pelo autor. Dicas sobre o uso do HDL Coder com Filtros digitais PASSO 1 - Criar um subsistema DUT da parte que se deseja converter em HDL. PASSO 2 - Configurar todos os blocos dentro desse DUT como processamento por amostra. Configurar também as entradas de sinal como PASSO 3 - Use o Filter Realization Wizard para realizar os filtros usando a quantização em ponto fixo. PASSO 4 - Faça a conversão para HDL. Veja detalhes em [Limitações do HDL coder]. Atividades do projeto final Implementar um modelo de entrada e saída que permita testar o sistema de processamento de sinais (DUT) a ser inserido entre eles. O DUT deverá ser composto de no mínimo 4 filtros, projetados conforme já realizado nas atividades AE1 e AE2. 1 filtro IIR HP ou LP (usando aproximação de Chebychev tipo 1 ou Butterworth) 1 filtro IIR BS ou BP (usando aproximação de Chebychev tipo 2 ou Eliptico) 1 filtro FIR HP ou LP (algoritmo de Park-McCleallan ou janela fixa ou janela ajustável) 1 filtro FIR BS ou BP (algoritmo de Park-McCleallan ou janela fixa ou janela ajustável) A frequência de amostragem e frequências de passagem e rejeição, assim como os ganhos, e as conexões entre o sinal de entrada e a(s) saída(s) fica a critério do aluno. O sistema deverá ser simulado no Simulink, mostrando claramente que o sistema realiza a função desejada, através da analise no domínio do tempo (DT) e no domínio da frequencia (DF). Após realizadas as simulações os filtros usados deverão ser agrupados em um subsistema para gerar o código VHDL do sistema DUT e também o testbench para este DUT. Para a geração do VHDL a partir do Simulink, siga o procedimento indicado em aula, mas também confira o link [Limitações do HDL coder]. Entregas Um arquivo comprimido (.zip .rar .tgz .gz) contendo o modelo no Simulink, o DUT em VHDL com o testbench Um relatório contendo a descrição do sistema e os resultados obtidos, com análise de cada figura inserida.