%% Estudo sobre filtros FIR
%
% Material usado na aula de PSD do CST em Sistemas de Telecomunicações do IFSC Campus São José
%
% Autor: Professor Diego da Silva de Medeiros
%

clear all; close all; clc

%% Resposta em frequência do filtro
%
% Para falar sobre filtros FIR, iniciando projetando um filtro passa baixas com frequência de corte wc
%

% Frequência de corte do filtro
wc = 0.1*pi;

% Resolução do eixo das frequências (número de amostras)
resolucao = 10000;

% Eixo das frequências
w = linspace(-pi,pi,resolucao);

% Resposta em frequência do filtro
H = u(w+wc) - u(w-wc);

% Eixo da frequência a ser mostrado
eixoF = [-0.5 0.5 -0.5 1.5];

% Plot dos resultados
figure; plot(w/pi,H); axis(eixoF); grid on
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude do filtro')

%% Resposta ao impulso do filtro
%
% Fazendo a transformada inversa de Fourier, encontra-se a resposta ao impulso h(n) do filtro:
%
%      h(n) = (wc/pi) Sinc(wc n)
%
% É necessário fazer um ajuste na função sinc(x), pois ela para o Matlab é Sen(pi x)/(pi x)
%

% Número de amostras do sinal no tempo
N = resolucao;

% Criação do eixo do tempo
n = -N/2:N/2 - 1;

% Definição da resposta ao impulso do filtro (ver anotações no caderno)
h = (wc/pi) * sinc(wc*n/pi);

% Tamanho do ciclo da função sinc(x) calculada
c = pi/wc;
% Número de ciclos da função sinc(x) a ser mostrados (número par)
Nc = 8;
% Eixo do tempo a ser mostrado
eixoN = [-Nc*c, Nc*c, -(wc/pi)*0.3, wc/pi];

% Plot dos resultados
figure; stem(n,h,'r','Marker','none'); hold on; plot(n,h); grid on; axis(eixoN)
xlabel('n'); ylabel('h[n]'); title('Resposta ao impulso do filtro')

%% Resposta ao impulso do filtro truncado
%
% Para projetar o filtro passa baixas através de um filtro FIR, é necessário limitar o número de amostras da resposta ao impulso.
% Uma forma de fazer isso é simplesmente truncar a resposta para valores maiores que um dado M.
%

% Número de ciclos mantidos na resposta truncada
Mc = 4;
% Número de amostras da resposta ao impulso truncada (2*M+1)
M = Mc*c;

% Definição da resposta ao impulso do filtro truncada, usando janelas
ht = (wc/pi) * sinc(wc*n/pi) .* (u(n+M) - u(n-M-1));

% Plot dos resultados
figure; stem(n,ht,'r','Marker','none'); hold on; plot(n,ht); grid on; axis(eixoN)
xlabel('n'); ylabel('h[n]'); title('Resposta ao impulso do filtro truncada')

%% Resposta em frequência do filtro truncado
%
% Verificando o efeito do truncamento na resposta em frequência do filtro
%

% Cálculo da resposta em frequência do filtro
Ht = fftshift(fft(ht));

% Plot dos resultados
figure;  semilogy(w/pi,abs(Ht))
axis(eixoF); grid on
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude do filtro')

%% Resposta em frequência da janela retangular
%
% Verificando a resposta em frequência da janela retangular pra explicar o efeito do truncamento da resposta ao impulso na
% resposta em frequência do filtro
%

% Cálculo da resposta em frequência do filtro (ver anotações)
Wr = 2*M*sinc(w*M/pi);

% Plot dos resultados
figure; plot(w/pi,Wr); grid on;
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude da janela retangular de truncamento')

%% Plotando os sinais da convolução
%
% A resposta em frequência do filtro truncado é a convolução da resposta em frequência do filtro ideal com a resposta em
% frequência da janela. Plota os dois sinais e o resultado para fins de comparação.
%

figure

% Filtro ideal
subplot(3,1,1); plot(w/pi,H); grid on
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude do filtro')

subplot(3,1,2); plot(w/pi,Wr); grid on;
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude da janela retangular de truncamento')

subplot(3,1,3); plot(w/pi,abs(Ht)); grid on
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude do filtro')

%% Outras janelas
%
% A fim de reduzir o fenômeno de Gibbs causado pela janela retangular, outras formas de truncamento foram desenvolvidas. Estas se
% dão através de janelas. Aqui várias dessas janelas serão apresentadas
%

% Para realizar a limitação das janelas no intervalo de -M a M, cada janela é multiplicada por uma retangular neste intervalo
janela = u(n+M) - u(n-M-1);

% Eixo do tempo a ser mostrado
eixoN = [-Nc*c, Nc*c, -0.3, 1.1];

%% Janela de Bartlett (triangular)

% Janela
w_bar = 1 - abs(n)/(M+1);
% Limitação de -M a M
w_bar = w_bar .* janela;

% Calculando a resposta em frequência da janela
W_bar = fftshift(fft(w_bar));

% Eixo da frequência a ser mostrado
eixoF = [-0.4 0.4 -1 max(abs(W_bar))];

% Plotando as janelas
figure;
subplot(2,1,1); stem(n,w_bar,'r','Marker','none'); hold on; plot(n,w_bar); grid on; axis(eixoN)
xlabel('n'); ylabel('w[n]'); title('Janela de Bartlett')
subplot(2,1,2); semilogy(w/pi,abs(W_bar)); grid on; axis(eixoF)
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude da janela de Bartlett')

%% Janela de Hann

% Janela de Hann
w_han = (1/2)*(1 + cos(2*pi*n/(2*M+1)));
% Limitação de -M a M
w_han = w_han .* janela;

% Calculando a resposta em frequência da janela
W_han = fftshift(fft(w_han));


% Eixo da frequência a ser mostrado
eixoF = [-0.4 0.4 -1 max(abs(W_han))];

% Plotando as janelas
figure;
subplot(2,1,1); stem(n,w_han,'r','Marker','none'); hold on; plot(n,w_han); grid on; axis(eixoN)
xlabel('n'); ylabel('w[n]'); title('Janela de Hann')
subplot(2,1,2); semilogy(w/pi,abs(W_han)); grid on; axis(eixoF)
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude da janela de Hann')

%% Janela de Hamming

% Janela de Hamming
w_ham = 0.54 + 0.46*cos(2*pi*n/(2*M+1));
% Limitação de -M a M
w_ham = w_ham .* janela;

% Calculando a resposta em frequência da janela
W_ham = fftshift(fft(w_ham));


% Eixo da frequência a ser mostrado
eixoF = [-0.4 0.4 -1 max(abs(W_ham))];

% Plotando as janelas
figure;
subplot(2,1,1); stem(n,w_ham,'r','Marker','none'); hold on; plot(n,w_ham); grid on; axis(eixoN)
xlabel('n'); ylabel('w[n]'); title('Janela de Hamming')
subplot(2,1,2); semilogy(w/pi,abs(W_ham)); grid on; axis(eixoF)
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude da janela de Hamming')

%% Janela de Blackman

% Janela de Blackman
w_bla = 0.42 + 0.5*cos(2*pi*n/(2*M+1)) + 0.08*cos(4*pi*n/(2*M+1));
% Limitação de -M a M
w_bla = w_bla .* janela;

% Calculando a resposta em frequência da janela
W_bla = fftshift(fft(w_bla));


% Eixo da frequência a ser mostrado
eixoF = [-0.4 0.4 -1 max(abs(W_bla))];

% Plotando as janelas
figure;
subplot(2,1,1); stem(n,w_bla,'r','Marker','none'); hold on; plot(n,w_bla); grid on; axis(eixoN)
xlabel('n'); ylabel('w[n]'); title('Janela de Blackman')
subplot(2,1,2); semilogy(w/pi,abs(W_bla)); grid on; axis(eixoF)
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude da janela de Blackman')


%% Resposta em frequência do filtro janelado
%
% Verificando o efeito da multiplicação da resposta ao impulso do filtro pelas janelas

% Multiplicando a janela pelo filtro ideal
h_bar = h .* w_bar;
h_han = h .* w_han;
h_ham = h .* w_ham;
h_bla = h .* w_bla;

% Calculando a resposta em frequência dos filtros janelados
H_bar = fftshift(fft(h_bar));
H_han = fftshift(fft(h_han));
H_ham = fftshift(fft(h_ham));
H_bla = fftshift(fft(h_bla));

% Plotando as respostas em frequência
eixoF = [-1 1 10^(-10) 1];
figure
subplot(2,2,1); semilogy(w/pi,abs(H_bar)); grid on; axis(eixoF)
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude do filtro Bartlett')
subplot(2,2,2); semilogy(w/pi,abs(H_han)); grid on; axis(eixoF)
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude do filtro Hann')
subplot(2,2,3); semilogy(w/pi,abs(H_ham)); grid on; axis(eixoF)
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude do filtro Hamming')
subplot(2,2,4); semilogy(w/pi,abs(H_bla)); grid on; axis(eixoF)
xlabel('Frequência normalizada [\omega/\pi]'); ylabel('|H(\omega)|'); title('Resposta em magnitude do filtro Blackman')

















