PSD29007-Engtelecom(2016-1) - Prof. Marcos Moecke

De MediaWiki do Campus São José
Ir para: navegação, pesquisa

Registro on-line das aulas

Unidade 1

Aula 1 (22 Mar)
  • Revisão de Sinais e Sistemas no tempo discreto em Matlab:
  • Resposta de sistemas LTI (Experimento 1.1)
  • Relembrar o conceito de equação de diferenças de um sistema LTI discreto e resposta ao impulso.
  • Resposta ao delta de Kronecker do sistema LTI discreto
onde , e logo
%  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.1
alpha = 1.15; N = 256;
x = [1 zeros(1,N)];
y = filter(1,[1 -1/alpha],x);
stem(y);
  • 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 (3 Hz) e um sinal (7 Hz) são idênticas quando amostrado com um sinal de 10 Hz.
%  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')
Aula 2 (24 Mar)
  • Revisão de Sinais e Sistemas no tempo discreto em Matlab:
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')
  1. Verifique a diferença entre os tipos de plots comentados no código.
  2. substitua o denominador de H(z) por dois polos em [-0.8 -0.8].
  3. verifique o que ocorre se forem utilizados polos complexos conjugados [0.3-0.4i 0.3+0.4i 0.1]
  4. verifique o que ocorre se forem utilizados polos complexos não conjugados [0.3-0.4i 0.3+0.8i]
  5. verifique o que ocorre se os polos estiverem fora do circulo unitário [1.2 -0.2]. Interprete este resultado
Aula 3 (29 Mar)
  • Revisão de Sinais e Sistemas no tempo discreto em Matlab:
Aula 4 (31 Mar)
  • Revisão de Sinais e Sistemas no tempo discreto em Matlab:
  • Análise de Sinais (Experimento 3.2) - Análise de um sistema h[n] correspondente a um filtro passa-faixa, utilizando um sinal de entrada x[n] senoidal (ou um sinal r[n] de ruído branco). Análise da entrada x[n] e saída y[n] usando a fft.
Variação do Experimento 3.2
%% Variação do Experimento 3.2 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.
%
% Análise de sinais no domínio da frequência 
% File Exp3_2.m 

fs = 200;   % frequência de amostragem
f_sinal = 10;  A_sinal = 1;   % freqüência e amplitude do sinal 
T = 1;      % Duração do sinal
k_noise = 0;    % Intensidade do ruído
 
time = 0 : 1/fs : (T-1/fs);
L = length(time);
freq = time * fs/T;
 
sinal = A_sinal*sin(2*pi*f_sinal.*time);
noise = k_noise*randn(1,fs*T);
x = sinal + noise;
X = abs(fft(x))/L;
 
figure(1);
subplot(211);plot(time,x);
subplot(212);plot(freq,X);
  1. Acrescente a Figura 1 um plot com a magnitude em dB do sinal no domínio da frequência - 20*log10(X)
  2. Insira nos gráficos títulos para cada subplot, labels para os eixos X e Y, e posicione o texto "F Hz" para indicar o pico nos gráficos 2 e 3, conforme mostrado na figura abaixo.
DTxDF sinal noise.png

Figura 1 - Análise no domínio da frequência do sinal

  1. Varie o valor de k entre 0 e 2 (com passo de 0.1) e analise o sinal no domínio do tempo e no domínio da frequência.
  2. Utilize k = 0.3 e varia a frequência do sinal entre 0 até 200 Hz (com passo de 10 Hz). Interprete os resultados obtidos.
  • Consulte a documentação do Matlab sobre
     grid, subplot, xlabel, ylabel, xlim, ylim, title, log10, log
    
  • Ver pag. 141 a 145 e 230 a 235 de [1]

Unidade 2

Aula 5 (5 Abr)
  • Filtros Analógicos:
  • Aproximação de magnitude de filtros analógicos: do tipo Butterworth.
  • Ver pag. 186 a 204 de [2]
Aula 6 (7 Abr)
  • Filtros Analógicos:
  • Projeto de filtros analógicos passa-baixas: do tipo Butterworth. (continuação)
  • Ver pag. 194 a 204 de [2]
Aula 7 (12 Abr)
  • Filtros Analógicos:
  • Projeto de filtros analógicos passa-baixas: do tipo Butterworth. (continuação)
  • Projeto de filtros analógicos passa-baixas: do tipo Chebyshev I.
  • Ver pag. 204 a 208 de [2]
Aula 8 (14 Abr)
  • Filtros Analógicos:
  • Exemplos de projeto de filtro passa-baixas com frequência de passagem de 16 krad/s com atenuação máxima de 0.3 dB, frequência de rejeição de 20 krad/s com atenuação mínima de 20 dB; e ganho em DC de 3 dB.
%% 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');
Aula 9 (19 Abr)
  • Filtros Digitais: Filtros IIR:
  • Transformação de frequência de filtros analógicos
(passa-baixas -> passa-baixas, passa-baixas -> passa-altas, passa-baixas -> passa-faixa, passa-baixas -> rejeita-faixa)
  • Funções para projeto do filtro protótipo analógico passa-baixas: besselap, buttap, cheb1ap, cheb2ap, ellipap
  • Funções de transformação de frequencia: lp2bp, lp2bs, lp2hp, lp2lp
  • Ver pag. 208 a 218 de [2]
Aula 10 (26 Abr)
  • 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)
  • Ver pag. 219 a 229 de [2]
  • Ver pag. 403 a 415 e 434 a 435 de [1]
Aula 11 (28 Abr)
  • Filtros Digitais: Filtros IIR: Uso do Matlab.
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 e denominador de .

Unidade 3

Aula 12 (3 Mai)
  • Filtros Digitais: Filtros FIR
  • Filtros de fase linear: simétricos e antisimétricos
  • Ver pag. 249 a 256 de [2]
Aula 13 (5 Mai)
  • Filtros Digitais: Filtros FIR
  • Filtros de fase linear: propriedades
  • Coeficientes da série de Fourier de filtros ideias: LP, HP, BP, BS
  • Passa-baixas (Low-pass)
  • Passa-altas (High-pass)
  • Passa-faixa (Band-pass)
  • Rejeita-banda (Band-stop)
  • Janela retangular, fenômeno de Gibbs
  • Estudar no Matlab as funções wintool, wvtool, window
  • Uso de funções de janelamento temporal no projeto de filtros digitais.
  • Tipos de janelas temporais usadas no projeto de filtros digitais.
  • Retangular
  • Bartlett
  • Hanning
  • Hamming
  • Blackman
  • em todas as janelas quando
  • Ver pag. 256 a 265 de [2]
Aula 16 e 17 (10 e 12 Mai)
L = 64; 
wvtool(rectwin(L), triang(L), bartlett(L), hann(L), hamming(L), blackman(L), blackmanharris(L), nuttallwin(L));
Janela
Retangular 13.3 {{{4}}}
Triangular 26.6 {{{4}}}
Barlett 26.5 {{{4}}}
Hann 31.5 {{{4}}}
Barlett-Hanning 35.9 {{{4}}}
Hamming 42.5 {{{4}}}
Bohman 46.0 {{{4}}}
Parzen 53.1 {{{4}}}
Backman 58.1 {{{4}}}
Flat Top 88.0 {{{4}}}
Backman-Harris 92.1 {{{4}}}
Nutfall 93.8 {{{4}}}
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 e da ordem do filtro são obtidos por:

onde é a atenuação do lóbulo lateral e é a largura da banda de transição em rad/amostra.

Utilizando o Matlab é possível estimar esses valores utilizando a função kaiserord. Exemplo da obtenção de um filtro passa baixa com , , atenuação de 40 dB na "stopband"

fsamp = 8000;
fcuts = [1000 1500];
mags = [1 0];
devs = [0.01 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp);


Com os parâmetros é possível projetar o filtro usando a função fir1, que utiliza o método da janela para o projeto do filtro.

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)
  • Ver as funções fir1, kaiserord do Matlab.
  • Início da AE2. (Os alunos devem completar a atividade extra-classe)
  • Ver pag. 266 a 273 de [2]
Aula 18 (17 Mai)
  • Filtros Digitais: Filtros FIR
  • Aula Prática: Uso do [1] Fdatool para projeto de filtro IIR, FIR equiripple e FIR com janela.
%% Exemplo de Filtro 
fp = 3000 Hz;
fr = 4000 Hz;
fs = 20000 Hz;
Ap = 1 dB;
Ar = 40 dB;
Aula 19 (19 Mai)
  • Filtros Digitais: Filtros FIR
  • AE2 - Projeto de Filtro Digitais FIR - MATLAB

Unidade 4

Aula 20 (21 Mai)
  • Realização de Filtros
  • Realização de filtros FIR: Forma Direta.
FIR FD MathWorks.png
Figura 1 - Realização de filtros FIR na Forma Direta
  • Realização de filtros FIR: 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.
FIR FDT MathWorks.png
Figura 2 - Realização de filtros FIR na Forma Transposta
FIR FDT2 MathWorks.png
Figura 3 - Realização de filtros FIR na Forma Transposta
  • Realização de filtros FIR de fase linear: simétrico tipo I e II e antissimétrico tipo III e IV.
FIR Sym2 MathWorks.png
Figura 4 - Realização de filtros FIR de fase linear Simétrico I
FIR Sym1 MathWorks.png
Figura 5 - Realização de filtros FIR de fase linear Simétrico II
FIR AntiSym3 MathWorks.png
Figura 6 - Realização de filtros FIR de fase linear Antisimétrico III
FIR AntiSym4 MathWorks.png
Figura 7 - Realização de filtros FIR de fase linear Antisimétrico IV
  • Realização de Filtros usando o comando realizemdl do MatLab
Fs = 40000;              % 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(Hd,'MapCoeffsToPorts','on');
Aula 22 (24 Mai)
  • Realização de filtros FIR: Cascata, Polifase
  • Vantagens do uso de filtro Polifase:
1) Quando o sinal será subamostrado (downsampling) de "D" amostras após a filtragem, a complexidade da implementação é reduzida de "D" vezes, pois apenas uma das "fases" precisa ser implementada.
2) Para reduzir o harware a ser implementado, é possível implementar apenas uma das "fases" do filtro e trocar "D" vezes os coeficientes.
  • Realização de filtros IIR de 2ª ordem: Forma Direta I e II, e Forma Transposta I e II.
  • Separando H(z) em dois blocos , e obtendo o sinal intermediário W(z) ou Y(z) dependendo da ordem dos blocos.
H1 H2 MathWorks.png
Figura 7 - Separação do filtro IIR H(z) em H1(z) e H2(z)
Com o ordenamento dos blocos e em ordem direta teremos a Forma Direta I:
Podemos obter a realização de na forma direta.
Para obter a realização de , é necessário reescrever a saída em função de e das saídas anteriores e :
IIR FD1 MathWorks.png
Figura 8 - Realização de filtros IIR na Forma Direta I
Com o ordenamento dos blocos e em ordem reversa teremos a Forma Direta II:
IIR FD2a MathWorks.png
Figura 9 - Realização de filtros IIR na Forma Direta II
Considerando que os sinais no centro são idênticos podemos simplificar e obter a Forma Direta II (Canônica):
IIR FD2b MathWorks.png
Figura 10 - Realização de filtros IIR na Forma Direta II Canônica
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.
IIR FT1 MathWorks.png
Figura 11 - Realização de filtros IIR na Forma Transposta I
IIR FT2 MathWorks.png
Figura 12 - Realização de filtros IIR na Forma Transposta II
  • 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.
Aula 23 (31 Mai)
  • Filtros Digitais: Ferramentas do Matlab para projeto
  • Filtros Digitais: Utilização de filtros FIR
  • Um sinal DTMF com duração de 1 segundo com frequência de amostragem de 8 kHz, correspondente aos dígitos 1234567890 ('Dtmf.wav').
  • Um sinal contendo ruído branco com duração de 5 segundo com frequência de amostragem de 8 kHz ('RuidoBranco.wav').
  • Um sinal onda quadra com duração de 2 segundo com frequência de amostragem de 8 kHz e período de 2 ms ('Quadrada.wav').

Utilize o Matlab para gerar o seguinte sinal:

  • Um sinal de varredura de Cosseno entre 0 Hz e 4 kHz com duração de 1 segundo.
Fs = 8000;
t = 0:1/Fs:1-1/Fs;
x = chirp(t,0,1,4000);
spectrogram(x, blackman(128), 120, 200, Fs,'yaxis')
audiowrite('Chirp0-4kHz.wav',x,Fs)
sound(x,Fs)  % Atenção remova o fone de ouvido antes de realizar este procedimento.

Utilizar o Matlab para projetar os seguintes filtros FIR e transmitir os sinais obtidos acima. Em todos os filtros considere a frequência de amostragem como 8 kHz, atenuação máxima na banda de passagem de 0,5 dB, e as bandas de transição como 400 Hz. Use a janela de Kaiser.

  • Filtro passa-baixas com fc = 1,5 kHz e atenuação de 60 dB na rejeição;
  • Filtro passa-altas com fc = 2,5 kHz e atenuação de 30 dB na rejeição;
  • Filtro passa faixa com fc1 = 1,5 kHz e fc2 = 2,5 kHz e atenuação de 80 dB na rejeição;
  • Filtro rejeita faixa com fc1 = 1,9 kHz e fc2 = 2,1 kHz e atenuação de 80 dB na rejeição;

Após obter os filtros, transmita cada um dos sinais gerados no Audacity através do filtro e verifique o resultado obtido analisando os sinal obtidos comparando o espectrograma com a resposta em magnitude do filtro.

[x, Fs] = audioread('Chirp0-4kHz.wav');  % Leitura do sinal
t = (0:length(x)-1)/Fs;  % Vetor de tempo
b = fir1(48,[2000 2100]/Fs); % Filtro Hamming com ordem 48 passa faixa

sound(x,Fs)  % Atenção remova o fone de ouvido antes de realizar este procedimento.
y = filter(b,1,x);
sound(y,Fs)  % Atenção remova o fone de ouvido antes de realizar este procedimento.
subplot(311); spectrogram(x, blackman(128), 100, 200, Fs)
subplot(312); spectrogram(y, blackman(128), 100, 200, Fs)
[Hw, w] = freqz(b,1,2000);
subplot(313); plot(w/pi*Fs/2,20*log10(abs(Hw))); ylim([-100,0]);
Aula 24 (2 Jun)
  • Projeto (coletivo) de um receptor DTMF.
  • cada aluno deverá projetar dois discriminadores de frequências correspondente a uma linha e uma coluna do sistema DTMF.
  • As especificações do discriminador de frequência, mostrado na figura, são:
DiscriminadorDTMF.png
  • A frequência de amostragem f_s do sinal de entrada é de 8 kHz.
  • Os filtros passa banda (BP) deverão ter, inicialmente, uma largura de banda BW' correspondente a 10% da frequência central f_0.
  • Os filtros passa baixa (LP) deverão ter, inicialmente, uma freqüência de passagem f_p de 100 Hz.
  • O circuito retificador deve se implementado pela função abs.

Ver as especificações DTMF em:

Aula 25 (7 Jun)
  • Atraso de grupo em filtros IIR e FIR no Matlab}}
  • O atraso de grupo de um filtro é a medida da atraso médio do filtro em função da frequência do sinal de entrada. Ele é obtido pela primeira derivada da resposta de fase do filtro. Se a resposta em frequencia é , então o atraso de grupo é:
onde é a fase de .
  • Um filtro sem distorção de fase (Não causal) pode ser obtido ao passar uma sequencia x(n) por um filtro H1, tomando a saída do filtro revertida e passando novamente pelo mesmo filtro H1. A saída do último filtro revertida corresponde ao sinal x(n) filtro com fase zero. O filtro obtido desta forma tem as seguintes características:
  • A Distorção de fase nula
  • A função de transferência do filtro é igual a magnitude ao quadrada da função de transferência original do filtro H1.
  • A ordem do filtro é o dobro da ordem do filtro H1.
%% Carregando um sinal de ECG com ruído com duração de 4 segundos.
load noisyecg.mat
x = noisyECG_withTrend;
fa = 500;  %% 2000 amostras em 4 segundos => 500 amostras por segundo.
t = [0:length(x)-1]*1/fa;
plot(t,x);

%% Projetando um filtro passa-baixa tipo IIR  butter com f_passagem = 0.15 rad/s
d = designfilt('lowpassiir', ...
    'PassbandFrequency',0.15,'StopbandFrequency',0.2, ...
    'PassbandRipple',1,'StopbandAttenuation',60, ...
    'DesignMethod','butter');
freqz(d)

%% Filtro de x revertido x e somando com x filtrado. OFF LINE
y = flip(filter(d,flip(filter(d,x))));
y1 = filter(d,x);

figure(2);
subplot(2,1,1)
plot(t, [y y1])
title('Filtered Waveforms')
legend('Zero-phase Filtering','Conventional Filtering')

subplot(2,1,2)
plot(t, [x y])
title('Original Waveform')
legend('noisy ecg ','fitered ecg')
  • Verifique também o resultado da filtragem usando um filtro IIR (ellip, cheby1 ou cheby2) e filtros FIR (equiripple e de janela)
%% Projetando um filtro passa-baixa tipo FIR  equiripple com f_passagem = 0.15 rad/s
d = designfilt('lowpassfir', ...
    'PassbandFrequency',0.15,'StopbandFrequency',0.2, ...
    'PassbandRipple',1,'StopbandAttenuation',60, ...
    'DesignMethod','equiripple');

y = flip(filter(d,flip(filter(d,x))));
  • Note que nos filtros FIR de fase linear o procedimento mais simples é adiantar o sinal de acordo com o atraso de grupo (metade da ordem do filtro), devendo-se tomar cuidado para arredondar a meia amostra nos filtros de ordem impar.
y1 = filter(d,x);
gd = grpdelay(d);
gd1 = ceil(gd(1));
y = [y1(gd1:end); zeros(gd1-1, 1)];

O cálculo do atraso de grupo pode ser realizado utilizando a função grpdelay ou diretamente pela definição da derivada do ângulo em relação a frequência:

%% Calculo do atraso de grupo 
% Método 1 - uso da função grpdelay
[z,p,k] = butter(30,0.2);
sos = zp2sos(z,p,k);
[gd,w]=grpdelay(sos,128);
figure(1)
plot(w/pi,gd),grid on;
% Método 2 - derivada obtida por aproximação discreta
% calculo a cada par de pontos (w2-w1)/delta_w
[h,w] = freqz(sos);
a = unwrap(angle(h));
hold on; plot(w/pi,a,'g');
delta_w = pi/length(a);
plot(w(1:end-1)/pi+delta_w/2,-(a(2:end)-a(1:end-1))/delta_w,'r');
Aula 24 (9 Jun)
  • Projeto de filtros Passa-Tudo para equalizar a fase (atraso de grupo) de filtros IIR.
  • Filtros de fase mínima.
Aula 25 (14 Jun)
  • Ponto Fixo e Ponto Flutuante
  • Aritmética Binária
  • Modulo Arithmetic
  • Casting
  • Fixed-point in MATLAB
  • Quantização do Filtro Digital
filtSpecs = fdesign.lowpass('Fp,Fst,Ap,Ast',0.2,0.25,0.5,40);
%% Design IIR butterworth filter
Hiir = design(filtSpecs,'butter','MatchExactly','passband');
Hiir
Hiir.Arithmetic = 'fixed';
Hiir
%get(Hiir)
%% Design FIR equiripple filter
Hfir = design(filtSpecs,'equiripple')
Hfir
Hfir.Arithmetic = 'fixed';
Hfir

Ver mais em:

Unidade 5

Aula 26 (16 Jun)
  • Uso do Simulink
Introduction: What Is Simulink? 4:42
Constructing and Running a Simple Model 13:45
Simulating a Model 10:10
Aula 27 (21 Jun)
  • Uso do Simulink
  • Bases da modelagem gráfica com Simulink (45 minutos)
Working with MATLAB 9:12 - Pass data between Simulink and MATLAB
Creating Subsystems 6:46 - Simplify your model by grouping blocks into subsystems
  • Usando o Simulink para modelagem de Sistemas Dinâmicos Discretos(60 minutos)
Modeling Discrete Dynamical Systems 19:35 - Learn to use the Integer Delay and Discrete Filter blocks to model difference equations
Use DSP System Toolbox 12:01 - Explore the basics of DSP System Toolbox, and model a simple system
Working with Signals in Simulink 9:59 - Learn frame-based processing and its benefits, and visualize signals in frequency domain
Applying a Filter 9:04 - Learn to model noise and implement a hand-designed filter to remove noise
Designing and Implementing a Filter 10:23 - Design a digital filter using the filter design tool and review signal processing application examples
dspstartup.m command
Aulas XX
  • Projeto do Gerador e Detector de DTMF com Simulink
Aula XX (12 Jul)
  • Uso do HDL Coder
a = fi(-1, true, 8, 0)
a.bin
a = fi(-128, true, 8, 0)
a.bin
a = fi(127, true, 8, 0)
a.bin
  • Ver também a função resize da ieee.numeric_std library.
  function RESIZE (ARG: SIGNED; NEW_SIZE: NATURAL) return SIGNED;
  -- Result: Resizes the SIGNED vector ARG to the specified size.
  --         To create a larger vector, the new [leftmost] bit positions
  --         are filled with the sign bit (ARG'LEFT). When truncating,
  --         the sign bit is retained along with the rightmost part.

  function RESIZE (ARG: UNSIGNED; NEW_SIZE: NATURAL) return UNSIGNED;
  -- Result: Resizes the SIGNED vector ARG to the specified size.
  --         To create a larger vector, the new [leftmost] bit positions
  --         are filled with '0'. When truncating, the leftmost bits
  --         are dropped.
  • Simulação do projeto no ModelSim-ALTERA.
Procedimentos para a Simulação

Abra o ModelSim:

/opt/altera/13.0sp1/modelsim_ase/bin/vsim &

Troque a pasta de trabalho para a pasta onde o Matlab gerou os arquivos .vhd e .do [File > Change Directory] ou digite na janela tcl:

cd  /tmp/mlhdlc_sfir/codegen/mlhdlc_sfir/hdlsrc

Execute o arquivo

do mlhdlc_sfir_fixpt_tb_compile.do

Edite o arquivo mlhdlc_sfir_fixpt_tb_sim.do, comentando as linhas que forçam a saida do Modelsim. Estas linhas foram criadas para a integração direta com o Matlab, mas ela não funciona com a versão do ModelSim que temos disponível.

edit mlhdlc_sfir_fixpt_tb_sim.do
#onerror {quit -f}
#onbreak {quit -f}
...
#quit -f

Agora execute o arquivo mlhdlc_sfir_fixpt_tb_sim.do, o qual irá adicionar os sinais a serem analisados na janela Wave e executará todos os comandos vhdl do arquivo *_tb.vhdl gerado pelo Matlab.

do mlhdlc_sfir_fixpt_tb_sim.do

Ao tinal da simulação a janela Transcript indicará se o teste passou:

# ** Note: **************TEST COMPLETED (PASSED)**************
#    Time: 20030 ns  Iteration: 1  Instance: /mlhdlc_sfir_fixpt_tb

ou se falhou:

# ** Note: **************TEST COMPLETED (FAILED)**************
#    Time: 20030 ns  Iteration: 1  Instance: /mlhdlc_sfir_fixpt_tb

Note ainda que os sinais de entrada (x_in), os sinais de saída (y_out e delayed_xout) e os dois sinais de referência gerados na simulação com o Matlab (y_out_ref e delayed_xout_ref) são mostrados como sequências de bits. Para melhorar a visualização mude o Formato desses sinais para analógico [Clique Direito do Mouse> Format > Analog (automatic)].

ATUAL

AULA XX (14 Jul)
  • Uso do HDL Coder
  • Implementar o filtro indicado abaixo e realizar a simulação com MODELSIM
Exemplo de um Filtro IIR
  • Código da função de filtragem
% FILE mlhdlc_iir_filter.m

function y = mlhdlc_iir_filter(x, sos, g)

% Declare persistent variables and initialize
numSections = numel(sos)/6;
persistent z
if isempty(z)	
	z = zeros(numSections, 2);
end

y = x;
for i=coder.unroll(1:numSections)
    curSOS = sos((i-1)*6+1:i*6);
    [y z(i,:)] = biquad_filter(y, curSOS(1:3), curSOS(4:6), z(i, :));
end
y = y * g;

end

function [y, z] = biquad_filter (x, b, a, z)
% a(1) is assumed to be 1
% Direct-form II implementation

tmp = x - z(1)*a(2) - z(2)*a(3);
y = z(2) * b(3) + z(1) * b(2) + tmp * b(1);
z(2) = z(1);
z(1) = tmp;

end
  • Código do testbench
%   Copyright 2011-2013 The MathWorks, Inc.

% FILE: mlhdlc_iir_filter_tb.m

clear mlhdlc_iir_filter;

% All frequency values are in MHz.

Fs = 100;  % Sampling Frequency

N   = 4;     % Order
Fc1 = 29.5;  % First Cutoff Frequency
Fc2 = 30.5;  % Second Cutoff Frequency
[z,p,k] = butter(N/2,[Fc1 Fc2]/(Fs/2),'bandpass');
[sos1,g1] = zp2sos(z,p,k);	     % Convert to SOS form
[num,den] = zp2tf(z,p,k);

L = 1000; 
Fs = Fs*1e6  %  Passa para MHz
t = (0:L-1)'/Fs;  
x = 0.5*sin(2*pi*30e6*t) + 0.5*cos(2*pi*20e6*t) + + 0.5*cos(2*pi*35e6*t);
rng('default'); % always default to known state  
x = x + .5*randn(size(x));  % noisy signal
y = zeros(size(x));
%%
% Call to the design
sos = sos1.';
g = g1;
for i=1:numel(x)
    y(i) = mlhdlc_iir_filter(x(i), sos(:), g);
end
%%
close all;
figure('Name', [mfilename, '_psd_plot']);
pwelch(x, 128);
hold on;
pwelch(y, 128);
yh = get(gca,'Children');
set(yh(1),'Color','r');
  • Ver também

Avaliações

  • Entrega dos diversos trabalhos ao longo do semestre.
  • Projeto Final. O projeto é avaliado nos quesitos: 1) Implementação do Sistema, 2) Documentação, 3) Avaliação Global do aluno no projeto.

Atividades extra

Neste tópico serão listadas as atividades extras que os alunos da disciplina deverão realizar ao longo do curso. É importante observar o prazo de entrega, pois os conceitos serão reduzidos conforme o atraso na entrega. Para a entrega no prazo os conceitos possíveis são (A, B, C, D). Entrega com até uma semana de atraso (B, C, D). Entrega com até duas semanas de atraso (C ou D). Entrega com mais de duas semanas de atraso (D).

PARA ENTREGAR

AE1 - Projeto de Filtros Digitais IIR (Prazo de entrega 31/05/2016)
Uma das metodologias de projeto de filtros digitais IIR, consiste em: (a) projeto de um filtro protótipo analógico passa-baixas H(p); (b) transformação em frequência do filtro H(p) -> H(s), obtendo o filtro LP, HP, BP, BS, conforme desejado; (c) transformação do filtro analógico em filtro digital H(s) -> H(z) utilizando a transformação Bilinear. Neste exercício avaliativo é solicitado que cada equipe de alunos (2 a 3) realize o projeto do um conjunto de filtros, seguindo os passos acima descritos. Os filtros que cada equipe irá projetar estão listados na tabela abaixo:
Equipe Filtro 1 Filtro 2 Filtro 3
Gustavo, Stephany, Fernando LP - Butter (fp = 5 kHz; fr = 20 kHz, fs = 44 KHz) LP - Cheby I (fp = 15 kHz; fr = 20 kHz, fs = 50 KHz) BP - Cheby II (fr1 = 6 kHz; fp1 = 8.2 kHz, fp2 = 14 kHz; fr2 = 20 kHz, fs = 50 KHz))
Roicenir, Ernani LP - Butter (fp = 8 kHz; fr = 20 kHz, fs = 44 KHz) LP - Cheby I (fp = 15 kHz; fr = 20 kHz, fs = 100 KHz) BP - Butter (fr1 = 6 kHz; fp1 = 8.2 kHz, fp2 = 8.7 kHz; fr2 = 20 kHz, fs = 100 KHz))
Ronaldo, Vinicius LP - Butter (fp = 10 kHz; fr = 20 kHz, fs = 44 KHz) HP - Cheby II (fr = 15 kHz; fp = 20 kHz, Ap = 2 dB, Ar = 40 dB, G_p = 0dB, fs = 80 KHz) BS - Ellip (fr1 = 7 kHz; fp1 = 8.2 kHz, fp2 = 8.7 kHz; fr2 = 10 kHz, fs = 80 KHz))
Giulio, Walter, Tiago LP - Butter (fp = 5 kHz; fr = 20 kHz, fs = 50 KHz) HP - Ellip (fr = 15 kHz; fp = 20 kHz, Ap = 0.8 dB, Ar = 40 dB, G_p = 0dB, fs = 50 KHz) BS - Butter (fr1 = 7 kHz; fp1 = 8.2 kHz, fp2 = 8.7 kHz; fr2 = 10 kHz, fs = 50 KHz))
Matias, Lucas, Markus LP - Butter (fp = 8 kHz; fr = 20 kHz, fs = 50 KHz) LP - Cheby II (fp = 15 kHz; fr = 20 kHz, fs = 200 KHz) BS - Cheby I (fr1 = 7 kHz; fp1 = 8.2 kHz, fp2 = 8.7 kHz; fr2 = 11 kHz, fs = 200 KHz))
onde:
LP - Passa Baixa, HP - Passa Altas, BP - Passa Faixa, BS - Rejeita Faixa
- frequencia de amostragem do sinal de entrada; - frequência de passagem; - frequência de rejeição, - Atenuação máxima na banda de passagem (dB), - Atenuação mínima na banda de rejeição (dB), - Ganho médio na banda de passagem (dB).
Butter - Aproximação tipo Butterworth, Cheby I - Aproximação tipo Chebyshev I, Cheby II - Aproximação tipo Chebyshev Inversa, Ellip - Aproximação tipo Cauer ou Eliptica.
O filtro 1 deve ter Ap = 1 dB, Ar = 60 dB, G_p = 10 dB
Os filtros 2 e 3 devem ter Ap = 0.8 dB, Ar = 40 dB, G_p = 0 dB
  • O projeto do filtro 1 deve apresentar o cálculo da ordem do filtro, dos polos e zeros do filtro, a equação de H(p), H(s), H(z), |H(jw)|^2, a magnitude e fase do filtro nas frequências de passagem e rejeição.
  • Para todos os filtros, apresente de modo gráfico o diagrama de zeros do filtro, o atraso de grupo, e a resposta em frequência do filtro (ganho em dB e fase) dos filtros (a) protótipo H(p), (b) Filtro analógico H(s) e Filtro digital H(z).
  • Obtenha o valor numérico em dB da atenuação em fp e fr. Assegure-se que os filtros obtidos atendem as especificações.
  • A resposta em frequência normalizada ou real devem ser feitas usando o Matlab.
  • Apresentar o gráfico do ganho em dB e da fase em cada caso com escalas corretas e com legendas
  • Deve ser apresentado o diagrama dos pólos e zeros dos filtros H(s) e H(z)
  • Escreva um relatório técnico em PDF mostrando os resultados obtidos e comentando os resultados obtidos.
  • O "Publish" pode ser utilizado, mas o arquivo entregue deve ser PDF e não HTML.
  • Envie o relatório e os arquivos ".m" utilizados para o email "moecke at ifsc.edu.br" com o Assunto: PSD29007 - AE1 - Projeto de Filtros Digitais IIR.
AE2 - Projeto de Filtro Digitais FIR - MATLAB (Prazo de entrega 09/06/2016)
1. Projeto os filtros digitais FIR com fase linear as seguintes características: tipo passa baixas; = 3 kHz; = 4 kHz; = 20 kHz; = 0.01; = 0.01; = 0 dB.
onde:
- frequencia de amostragem do sinal de entrada; - frequência de passagem; - frequência de rejeição, - Atenuação máxima na banda de passagem (linear), - Atenuação mínima na banda de rejeição (linear), - Ganho médio na banda de passagem (dB).
2. Verifique se é possível projetar o filtro usando as janelas do tipo Retangular, Hanning, 'Hamming, Blackmann, Kaiser, Chebyshev.
  • Utilize as funções firpmord e kaiserord para estimar a ordem do filtro.
  • Utilize as funções adequadas do Matlab para obter os coeficientes das janelas e em seguida utilize a função fir1 para obter os filtros.
  • Ajuste a ordem do filtro e frequência de passagem de modo a conseguir que cada filtro atenda as especificações iniciais.
3. Para cada tipo de janela, apresente de modo gráfico a resposta em frequência do filtro (ganho) de menor ordem que atende as especificações. Sobreponha os gráficos inserindo uma legenda adequada (indicando o tipo de janela e ordem). Utilizando um escala em dB (entre 10 dB e -80 dB) e frequência em kHz. Utilize uma mascara com cor diferenciada para indicar claramente a especificação do filtro, e crie um segundo gráfico mostrando claramente a banda de passagem conforme ilustrado nas figuras abaixo:

Resposta em frequência - Ganho em dB Detalhe da banda de passagem da resposta em frequência - Ganho em dB

6. Gere um arquivo "pdf" utilizando o Publish com os resultados e texto explicativo e envie o email "moecke at ifsc.edu.br" com o Assunto: PSD29007 - AE2 - Projeto de Filtro Digitais FIR - MATLAB.

JÁ ENCERRADAS

ESTUDOS SEM ENTREGA DE DOCUMENTAÇÃO

AL1 - Variação do Experimento 1.2

No Experimento 1.2 varie o valor da frequência de amostragem de 6 até 20 Hz e observe:

  1. Em qual frequência deixa de ocorrer recobrimento do sinal 2.
  2. O que ocorre quando a frequência é 6, 7, 14 Hz? Explique
  3. Qual deveria ser a frequência do sinal f_2 para que as amostras tomadas sejam coincidentes como o sinal f_1 para uma frequência de amostragem f_s? Reescreva a equação e verifique no Matlab.

Referências Bibliográficas:

  1. 1,0 1,1 1,2 1,3 1,4 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
  2. 2,0 2,1 2,2 2,3 2,4 2,5 2,6 2,7 2,8 2,9 SHENOI, B. A. Introduction to Digital Signal Processing and Filter Design. 1.ed. New Jersey: John Wiley-Interscience, 2006. 440 p. ISBN 978-0471464822
  3. LATHI, Bhagwandas P. Sinais e Sistemas Lineares. 2. ed. Porto Alegre: Artmed-Bookman, 2007. 856 p. ISBN 978-8560031139