%% prog 044 filtros

clearvars;
close all;
clc;

%% ruido blanco 
%ruido blanco es una secuencia aletorea de sonidos que posee todas las frecuencias

fs = 44100;   %frecuencia de mestreo
N = 2*fs+1;   %numero de muestras en la secuencia entre 0 y 2 segundos

%usaremos el comando rand() que da numeros aleatoreos entre 0 y 1
x = 2*rand(N,1) - 1;

%tiempo
t = (0:N-1)/fs;

%ejes
Tmin = t(1);
Tmax = t(N);
Amin = -1.1;
Amax =  1.1; 

figure(1)
plot(t,x);
title('ruido blanco');
xlabel('tiempo t (s)');
ylabel('presion sonora p (Pa)');
axis([Tmin, Tmax, Amin, Amax]);
grid on;
box on;

sound(x,fs);

%% filtro pasa bajo

fcpb = 100;           %frecuencia de corte
Q = 1/sqrt(2);        %factor de calidad o agudeza

K = tan(pi*fcpb/fs);  %factor de normalizacion

%coeficientes del filtro pasa bajo
b0 = (Q*K^2)/(Q*K^2 + Q + K);
b1 = (2*Q*K^2)/(Q*K^2 + Q + K);
b2 = (Q*K^2)/(Q*K^2 + Q + K);
a1 = (2*Q*(K^2-1))/(Q*K^2 + Q + K);
a2 = (Q*K^2 + Q - K)/(Q*K^2 + Q + K);

%prelocalizamos la seņal de salida
ypb = zeros(N,1);

%procesamiento principal
for n = 3:N
    ypb(n) = b0*x(n) + b1*x(n-1) + b2*x(n-2) - a1*ypb(n-1) - a2*ypb(n-2);
end;


figure(2)
plot(t,ypb);
title('ruido blanco pasa bajo');
xlabel('tiempo t (s)');
ylabel('presion sonora p (Pa)');
axis([Tmin, Tmax, Amin, Amax]);
grid on;
box on;

%normalizamos para que se escuche bien
ypb = ypb/max(abs(ypb));

sound(ypb,fs);


%% filtro pasa alto

fcpa = 10000;         %frecuencia de corte
Q = 1/sqrt(2);        %factor de calidad o agudeza

K = tan(pi*fcpa/fs);  %factor de normalizacion

%coeficientes del filtro pasa alto
b0 = ( Q )/( Q*K^2 + K + Q );
b1 = ( -2*Q )/( Q*K^2 + K + Q );
b2 = ( Q )/( Q*K^2 + K + Q );
a1 =  ( 2*Q*(K^2-1) )/( Q*K^2 + K + Q );
a2 = ( Q*K^2 - K + Q)/( Q*K^2 + K + Q );

%prelocalizamos la seņal de salida
ypa = zeros(N,1);

%procesamiento principal
for n = 3:N
    ypa(n) = b0*x(n) + b1*x(n-1) + b2*x(n-2) - a1*ypa(n-1) - a2*ypa(n-2);
end;


figure(3)
plot(t,ypa);
title('ruido blanco pasa alto');
xlabel('tiempo t (s)');
ylabel('presion sonora p (Pa)');
axis([Tmin, Tmax, Amin, Amax]);
grid on;
box on;

%normalizamos para que se escuche bien
ypa = ypa/max(abs(ypa));

sound(ypa,fs);

%tarea usando filtro pasa bajo y pasa alto hacer filtro pasa banda
%y un filtro rechaza banda



%% graficos en la frecuencia fft no va a ser evaluado en trabajos o pruebas

%el comando fft pasa del tiempo a la frecuencia
%letras mayusculas y minusculas son variables distintas
%minusculas para el dominio del tiempo y mayusculas para el dominio de la frecuencia

X   = abs(fft(x)/N);
YPB = abs(fft(ypb)/N);
YPA = abs(fft(ypa)/N);

df = fs/N;          %intervalode frecuencia
f = 0:df:(fs - df); %frecuencia

%ejes
minf = f(1);
maxf = f(round(fs/2));


figure(4)
subplot(2,1,1)
semilogy(f(1:round(fs/2)),X(1:round(fs/2)),f(1:round(fs/2)),YPB(1:round(fs/2)));
title('espectro violin y violin pasa bajo');
xlabel('frecuencia f(Hz)');
ylabel('presion sonora P(f) (Pa)');
legend('X(f)','YPB(t)');
axis([minf maxf 0 Amax]);
grid on;
box on;


subplot(2,1,2)
semilogy(f(1:round(fs/2)),X(1:round(fs/2)),f(1:round(fs/2)),YPA(1:round(fs/2)));
title('espectro violin y violin pasa alto');
xlabel('frecuencia f(Hz)');
ylabel('presion sonora P(f) (Pa)');
legend('X(t)','YPA(t)');
axis([minf maxf 0 Amax]);
grid on;
box on;