%% prog006 mi primer fft

clearvars;
close all;
clc;

%% inicializacion

%datos inciales
f = 5000;       %frecuencia señal de audio
fs = 44100;     %frecuencia de muestreo
dt = 1/fs;      %periodo de muestreo en el tiempo
A = 1.0;        %amplitud
t = 0:dt:100;    %tiempo
N = length(t);  %longitud del vector tiempo senhal

%ejes graficos tiempo
xmint = min(t);
xmaxt = max(t);
ymint = -1.5*A;
ymaxt =  1.5*A;

%% senhal en el tiempo

%onda senoidal en el tiempo
x = A*sin(2*pi*f*t);

%onda cosenoidal
%x = A*cos(2*pi*f*t);

%pulso cuadrado
%x(1:N) = 0;
%x(1: round(N/1e5)) = 5e4;

%grafico
figure(1)
plot(t,x);
title('onda senoidal');
xlabel('t (s)');
ylabel('x (Pa)');
axis( [xmint xmaxt ymint ymaxt]);
grid on;
box on;

%% senhal en el dominio de la frecuencia

%conversion al domio de la frecuencia usando FFT
X = fft(x)/N;

df = fs/N;          %intervalo de muestreo en frecuencia
f = 0:df:(fs - df) ;  %frecuencia
xminf = min(f);
xmaxf = max(f);
yminf = -A;
ymaxf =  A;

figure(2)

subplot(2,1,1);
plot(f,real(X));
title('parte real espectro onda senoidal bilateral');
xlabel('f (Hz)');
ylabel('real(X)(Pa)');
axis( [xminf xmaxf yminf ymaxf]);
grid on;
box on;

subplot(2,1,2);
plot(f,imag(X));
title('parte imaginaria espectro onda senoidal bilateral');
xlabel('f (Hz)');
ylabel('imag(X)(Pa)');
axis( [xminf xmaxf yminf ymaxf]);
grid on;
box on;

figure(3)
plot(f,abs(X));
title('módulo espectro onda senoidal bilateral');
xlabel('f (Hz)');
ylabel('|X|(Pa)');
axis( [xminf xmaxf yminf ymaxf]);
grid on;
box on;

figure(4)
plot(f,2*abs(X));
title('módulo espectro onda senoidal unilateral escala lineal');
xlabel('f (Hz)');
ylabel('|X|(Pa)');
axis( [xminf round(xmaxf/2) yminf ymaxf]);
grid on;
box on;

figure(5)
semilogy(f,2*abs(X));
title('módulo espectro onda senoidal unilateral escala logarírmica');
xlabel('f (Hz)');
ylabel('|X|(Pa)');
axis( [xminf round(xmaxf/2) 1e-6 ymaxf]);
grid on;
box on;


%% senhal contaminada por ruido

%ruido blanco aleatoreo usamos el comando rand(N)
%que genera una secuencia de numesros aleatoreos de longitud N

y = 2*(rand(1,N) - 0.5);

%grafico
figure(6)
plot(t,y);
title('ruido');
xlabel('t (s)');
ylabel('y (Pa)');
axis( [xmint xmaxt ymint ymaxt]);
grid on;
box on;

%transformada de fourier
Y = fft(y);
Y = Y/N;

figure(7)
semilogy(f,2*abs(Y));
title('módulo espectro ruido unilateral escala logarírmica');
xlabel('f (Hz)');
ylabel('|Y|(Pa)');
axis( [xminf round(xmaxf/2) 1e-6 ymaxf]);
grid on;
box on;

%la fft no se puede aplicar a una señal aleatorea a fin de determinar el
%espectro

%densidad espectral de potencia
Nit= 4*4410;
Nov = 0.0;

[Sxx,fxx] = pwelch(x,Nit,round(Nov*Nit),Nit,fs);%tarea leer el help
[Syy,fyy] = pwelch(y,Nit,round(Nov*Nit),Nit,fs);%tarea leer el help
                                                %pwelch y tfestimate

figure(8)
semilogy(fxx,Sxx)
title('ondasenoidal PSD Input');
xlabel('f (Hz)');
ylabel('Sxx(f) (Pa)');
legend('Sxx(f) input');
axis([0 fs/2 1e-6 1]);
grid on;
box on;

figure(9)
semilogy(fyy,Syy)
title('ruido PSD Input');
xlabel('f (Hz)');
ylabel('Syy(f) (Pa)');
legend('Syy(f) input');
axis([0 fs/2 1e-6 1]);
grid on;
box on;