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

Registro on-line das aulas

Unidade 1

Unidade 1

Aula 1 (29 jul) Apresentação da disciplina Auto inscrição na Plataforma Moodle de PSD29007 (engtelecom2019-2) Aula 2 e 3 (1 e 5 ago) Transformadas de Fourier Aula 4 (8 ago) 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(211); plot(1000*t(1:200),X(1:200)) title('Signal Corrupted with Zero-Mean Random Noise') xlabel('t (milliseconds)') ylabel('X(t)') %% 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); % 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(212); plot(f,P1) ylim([0 1.05]) title('Single-Sided Amplitude Spectrum of X(t)') xlabel('f (Hz)') ylabel('|P1(f)|') % 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) title('Single-Sided Amplitude Spectrum of S(t)') xlabel('f (Hz)') ylabel('|P1(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. Aula 5 (12 ago) 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 1 EAD Aula 6 (15 ago EAD) Revisão de Sinais e Sistemas no tempo discreto em Matlab: Filtros Digitais ATIVIDADE 1 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). %% 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); ATIVIDADE 2 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. Relate os resultados obtidos no comentário da atividade no Moodle. Aula 7 (19 ago) Exercício - Sinal DTMF com ruído Verifique se o Matlab está reproduzindo corretamente o som. %% Carregando o som clear, close, clc load handel; %% Reproduzindo o som sound(y,Fs) % Reproduzindo o som %soundsc(y,Fs) % Reproduzindo o som %player = audioplayer(y, Fs); %play(player); Usando o Matlab (ou Audacity) para gerar um sinal DTMF correspondente a um número N e adicionar um ruido ao sinal. Opcionalmente utilize um sinal DTMF gravado Utilizar uma frequência de amostragem de 8000Hz de fazer a duração do sinal igual a 2 segundos. Para adicionar o ruído utilize a função y = awgn(x,snr), ou y = x + nivel*randn(n). Observe este sinal no domínio do tempo (DT) e domínio da frequência (DF). %% Carregando o som clear, close, clc [y,Fs] = audioread('DTMF_8kHz.ogg'); %% Reproduzindo o som sound(y,Fs) %% Visualizando o som no DT time = [0:length(y)-1]'/Fs; plot(time',y'); xlabel('segundos'); xlim([0 time(end)]), ylim([-1 1]); %% Visualizando o som no DF Nfreq = length(y); freq = linspace(0,2*pi,Nfreq)'*Fs/pi/2; Y = fft(y,Nfreq)/Nfreq; plot(freq,abs(Y)); xlabel('Hertz'); xlim([0 Fs/2]); Utilizar no Audacity um sinal DTMF ("1234567890") com fa= 8kHz Visualizar no domínio do tempo e frequência. Realizar a filtragem passa-baixas com fc = 1066 Hz, (selecionar a maior atenuação permitida) Realizar a filtragem passa-faixa com f0 = 770 Hz e B = 70Hz (selecionar a maior ordem permitida) Repetir o procedimento anterior para um sinal de ruído branco. Consulte a documentação do Matlab sobre fft, ifft, fftshift, randn Consulte a documentação do Matlab sobre plot, grid, subplot, hold, xlabel, ylabel, title, legend, xlim, ylim, log10, log Consulte a documentação do Matlab sobre text, zp2tf, tf2zp, fftfilt, awgn Ver pag. 141 a 145 e 230 a 235 de

Unidade 2

Unidade 2

Aula 8 (22 ago) 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.$ Análise básica de filtros analógicos com 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}}$ b = [1 1]; a = [1 1 5]; [z1,p1,k]=tf2zp(b,a) z2 = roots(b); p2 = roots(a); zplane(b,a); %% freqs(b,a); %% 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) semilogx(ws,abs(h)); grid on; subplot(212) semilogx(ws,angle(h)/pi*180); grid on; Projeto de filtros analógicos do tipo Butterworth A aproximação de magnitude de filtros analógicos pode ser realizado usando as aproximações de Butterworth, Chebyshev (tipo 1 ou 2) e Cauer. Proposta de exercício Use os polinômios de Butterworth com ordens de 1 a 10 mostrados na tabela abaixo para obter os filtros $H_{n}(s)=1/B_{n}(s)$. 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)$ 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. Aula 9, 10 (16, 26 ago) 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. $A_{p}=3dB$ 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$. Considere que $\epsilon ={\sqrt {10^{0.1A_{p}}-1}}$, teremos $\epsilon =1$ Para projetar o filtro é necessário: 1) determinar a ordem $n$ do filtro: $n\geq {\frac {\log(10^{0.1A_{s}}-1)}{2\log \Omega _{s}}}$ 2) obter os polos do filtro: $p_{k}=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 {1}{D(p)}}$, onde $D(p)=\prod _{k=1}^{n}\left(p-p_{k}\right)$ 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.$ Para qualquer $A_{p}$ Teremos $\epsilon ={\sqrt {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)}$. 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.$ Ver Uso do calculo simbólico na Matlab Ver pag. 186 a 204 de Aula 11 e 12 (29 ago e 2 set) Projeto de filtros analógicos do tipo Chebyshev I. Polinômios de Chebyshev: $C_{n}(\Omega )=cos(ncos^{-1}(\Omega ))\!$ Os polinômios de Chebyshev de primeira ordem são definidos pela relação recursiva: $C_{0}(\Omega )=1\,\!$ $C_{1}(\Omega )=\Omega \,\!$ $C_{n}(\Omega )=2\times \Omega \times C_{n-1}(\Omega )-C_{n-2}(\Omega ).\,\!$ Os primeiros cinco polinômios de Chebyshev de primeira ordem são: $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 \,$ 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, $\epsilon ={\sqrt {10^{0.1A_{p}}-1}}$, $\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)/(10^{0.1A_{p}}-1)}}}{\cosh ^{-1}\Omega _{s}}}$ Em seguida obter 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)$ $\varphi _{2}={\frac {1}{n}}\sinh ^{-1}\left({\frac {1}{\epsilon }}\right)$ $\left|H(j\Omega )\right|^{2}=1{\text{para n impar}},{\frac {1}{1+\epsilon ^{2}}}{\text{para n par}}$ $\left|H(0)\right|^{2}={\frac {H_{0}^{2}}{1+\epsilon ^{2}C_{n}^{2}\left(\Omega _{0}\right)}}$ Para obter a função de transferência: $H(p)={\frac {H0}{D(p)}}$, onde $D(p)=\prod _{k-1}^{n}\left(p-p_{k}\right)$ onde $H_{0}=d_{0}*{\sqrt {\frac {1}{1+\epsilon ^{2}}}}$ se n é par $H_{0}=d_{0}$ se n é impar $d_{0}=\prod _{k-1}^{n}\left(-p_{k}\right)$ é o último termo do denominador D(p) 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'). Ler Comparison of Analog IIR Lowpass Filters em ellip Uso das funções freqs, "zplane", fvtool 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 Aula 13 e 14 (5 e 9 set) Filtros Analógicos: 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}}{Bs}}$ 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 {Bs}{s^{2}+\omega _{0}^{2}}}$ onde $B=\omega _{p2}-\omega _{p1}$ e $\omega _{0}={\sqrt {\omega _{p2}\omega _{p1}}}$ Uso das funções semilogx, semilogy,logspace, linspace. Ver em IIR Filter Design, Special Topics in IIR Filter Design. Funções para projeto do filtro protótipo analógico passa-baixas: buttap, cheb1ap, cheb2ap, ellipap Funções de transformação de frequencia: lp2bp, lp2bs, lp2hp, lp2lp Ver pag. 208 a 218 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_{p1}$ = 811 Hz, $f_{p2}$ = 895,5 Hz $f_{s1}$ = 770 Hz, $f_{s2}$ = 1209 Hz, $A_{p}$ = 1 dB, $A_{r}$ = 30 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_{r}$ = 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. Aula 15 (16 set) 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 entre 0 e 1, onde 1 corresponde a metade da frequência de amostragem $(fa/2)$ 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}}{Bs}}$ 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 {Bs}{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 Aula 16 (19 set) 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

Unidade 3 - Filtros FIR

Aula 16 (23 set) Filtros Digitais: Filtros FIR Filtros de fase linear: simétricos e antisimétricos (Tipo 1, 2, 3 e 4) Filtros de fase linear: propriedades (respostas em frequência possíveis, distribuição dos zeros em simetria quadrantal) Observar as propriedades dos FIR tipo 1, 2, 3 e 4. Observe a resposta de frequência, fase, atraso de grupo, coeficientes e 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)]; ... %% Tipo II - LP, BP % tem um zero em -1 b = [bi flip(bi)]; ... %% Tipo III - BP % tem um zero em 1 e -1 b = [bi 0 -flip(bi)]; ... %% Tipo IV - BP, HP % tem um zero em 1 b = [bi -flip(bi)]; ... Figura 1 - Propriedades do filtro FIR de fase linear (Tipo 1) ATUAL Aula 17 (26 set) Coeficientes da série de Fourier de filtros ideias: LP, HP, BP, BS Passa-baixas (Low-pass) $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.$ Passa-altas (High-pass) $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.$ Passa-faixa (Band-pass) $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.$ Rejeita-banda (Band-stop) $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.$ Janela retangular, fenômeno de Gibbs Estudar no Matlab as funções wintool, wvtool, window Ver pag. 249 a 256 de Ver FIR Filter Design Uso de 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$ Bartlett $w(n)=1-{\frac {\left|n\right|}{M+1}};\qquad -M\leq n\leq M$ Hanning $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)); Tabela 5.1 Janela $A_{sl}$ $A_{s}$ $\Delta \omega$ Retangular 13.3 20.33 0.92$\pi$/M Triangular 26.6 27.41 Bartlett 26.5 27.48 Hann 31.5 44.03 3.11$\pi$/M Bartlett-Hanning 35.9 40.77 Hamming 42.5 54.08 3.32$\pi$/M Bohman 46.0 51.84 7.01$\pi$/M Parzen 53.1 56.89 Blackman 58.1 75.25 5.56$\pi$/M 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$ Projeto de filtro FIR utilizando janelas temporais fixas. Exemplo de projeto Projetar um filtro passa baixas usando uma janela temporal fixa (verificar a janela que atende a especificação) wp = 0.2*pi; Ap = 0.2 dB; Gp = 0 dB ws = 0.3*pi; As = 60 dB; Informar qual o tipo de janela, a ordem obtida, e o valor de wc do projeto final Exemplo de projeto Projetar um filtro LP usando uma janela temporal fixa (hamming, bartlett-hanning, hanning). wp = 0.4*pi; Ap = 1 dB; Gp = 0 dB ws = 0.6*pi; As = 40 dB; Comparar os 3 tipos de janela, a ordem obtida, e o valor de wc em cada projeto. Use como uma estimativa inicial os valores da Tabela 5.1 pag. 268. 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 PASSO 3 - Calcule os coeficientes clp do filtro LP , calcule os valores da janela w e obtenha a resposta ao impulso do filtro h = clp * w. PASSO 4 - Ajuste o ganho do filtro para que a resposta na banda de passagem fique abaixo da especificação maxima. PASSO 5 - Verifique o valor real de Dwr = wAs-wAp, e faça a correção da ordem do filtro em função do desvio constatado. N2 = N*Dwr/Dw. PASSO 6 - Corrija o valor de projeto dos coeficientes Clp do filtro ideal, a janela e a resposta ao impulso. Repita o PASSO 3 até 6, até obter um filtro que atenda as especificações de Dw. PASSO 7 - Desloque a frequência de corte wc de modo a obter o valor correto de wp. wc2 = wp + (wp-wAp). Projeto de filtro FIR. Projete os dois filtros projetados anteriormente como IIR, utilizando 3 janelas diferentes. Compare os filtros obtidos com os filtros IIR. 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) Ver pag. 256 a 265 de 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

Avaliações

• Entrega dos diversas Atividades Extraclasse ao longo do semestre.
• Provas escritas A1 e A2
• Entrega do Projeto Final. O projeto é avaliado nos quesitos:

1) Implementação do Sistema,

2) Documentação,

3) Avaliação Global do aluno no projeto.

• Entrega dos Atividades Extraclasse ao longo do semestre AE1 a AE(N). A entrega, detalhes e prazos de cada AE serão indicados na plataforma Moodle

Referências Bibliográficas

Curso de Engenharia de Telecomunicações