%% Series de Fourier
% Dercargar FourierSeriesDemo.zip desde http://users.ece.gatech.edu/mcclella/matlabGUIs/
% Decomprimir y ejecutar
% No añadir al path, cambiar al directorio
%fseriesdemo
%
% Descargar la función apxfourier
f = evalin(symengine,'heaviside(x+0.5)-heaviside(x-0.5)');
%syms x real;f = heaviside(x+0.5)-heaviside(x-0.5);
N=20;T=4;
C=apxfourier(f,-T,N)
subplot(1,2,1);stem(sqrt(sum(C.^2)));title('Dominio de la Frecuencia');
xlabel('Termino k-esimo');ylabel('\itAmplitud c_k');axis([0,N,0,1]);
subplot(1,2,2);stem(atan2(C(1,:),C(2,:)));title('Dominio de la Frecuencia');
xlabel('Termino k-esimo');ylabel('\itDesfase \theta');axis([0,N,-3.5,3.5]);
disp('pulsa una tecla para el siguiente apartado');pause
close all
% usando la función numérica
fn1=@(x) (x>=-0.5 & x<=0.5);X=linspace(-T/2,T/2);Y=fn1(X);plot(X,Y);
C=apxfourier(fn1,-T,N)
% Representación
Y2=ones(size(X))*norm(C(:,1))/2;w=2*pi/T;
for k=2:N;
    Y2=Y2+C(1,k)*cos((k-1)*w*X)+C(2,k)*sin((k-1)*w*X);
    plot(X,[Y;Y2]);grid on ;title(['Armónico ',int2str(k-1)]);
    legend('Señal','Fourier');
    disp('Pulse una tecla para añadir el siguiente armónico');pause
end
disp('pulsa una tecla para el siguiente apartado');pause
close all
%
% Ejercicio 1, aumentamos el periodo y el numero de terminos
% usando la función simbolica anterior
N=50;T=20;
C=apxfourier(f,-T,N)
subplot(1,2,1);stem(sqrt(sum(C.^2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('Termino k-esimo');ylabel('\itAmplitud c_k');axis([0,N,0,1]);
subplot(1,2,2);stem(atan2(C(1,:),C(2,:)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('Termino k-esimo');ylabel('\itDesfase \theta');axis([0,N,-3.5,3.5]);
disp('pulsa una tecla para el siguiente apartado');pause

figure(3)
N=100;T=40;
C=apxfourier(f,-T,N)
subplot(1,2,1);stem(sqrt(sum(C.^2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('Termino k-esimo');ylabel('\itAmplitud c_k');axis([0,N,0,1]);
subplot(1,2,2);stem(atan2(C(1,:),C(2,:)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('Termino k-esimo');ylabel('\itDesfase \theta');axis([0,N,-3.5,3.5]);
disp('pulsa una tecla para el siguiente apartado');pause
%
% Ejercicio 2: funcion diente de sierra
fn2=@(x) (x>=-0.5 & x<=0.5)*2.*(0.5+x);X=linspace(-T/2,T/2,101);plot(X,fn2(X));
C=apxfourier(fn2,-T,N)

%
% Ejercicio 3: funcion triangular
fn3=@(x) (x>=-0.5 & x<=0)*2.*(0.5+x)+(x>0 & x<=0.5)*2.*(0.5-x);X=linspace(-T/2,T/2,101);plot(X,fn3(X));
C=apxfourier(fn3,-T,N)
disp('pulsa una tecla para el siguiente apartado');pause
close all
% usando la función simbolica anterior
%
%% Transformada de Fourier
clear all
close all
syms clear
% Definimos una función simbolica con una onda cuadrada en [-0.5, 0.5]
f = evalin(symengine,'heaviside(x+0.5)-heaviside(x-0.5)');ezplot(f,-5,5)
%syms x real;f = heaviside(x+0.5)-heaviside(x-0.5);ezplot(f,-5,5)
F= fourier(f)
a=40;W=linspace(0,a,201);FW=subs(F,W);
subplot(1,2,1);plot(W,abs(FW));title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itAmplitud');
subplot(1,2,2);plot(W,atan2(imag(FW),real(FW)));title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itDesfase \theta');
%
% Función diente de sierra
f = evalin(symengine,'(heaviside(x+0.5)-heaviside(x-0.5))*2*x');ezplot(f,-5,5)
%syms x real;f = (heaviside(x+0.5)-heaviside(x-0.5))*2*x;ezplot(f,-5,5)
F= fourier(f)
a=40;W=linspace(0,a,200);
FW=subs(F,W);
%
% Función triangulo
f = evalin(symengine,'(heaviside(x+0.5)-heaviside(x))*(1+2*x)+(heaviside(x)-heaviside(x-0.5))*(1-2*x)');ezplot(f,-5,5)
F= fourier(f)
a=40;W=linspace(0,a,200);
FW=subs(F,W);
%
f = evalin (symengine, 'piecewise ([x<=-1,0],[x>-1 and x<=0,1+x],[x>0 and x<=1,1-x],[x>1,0])'); % No admite ezplot
disp('pulsa una tecla para el siguiente apartado');pause
close all
%% Transformada rápida de Fourier
% Caso periódico
clear all
close all
syms clear
syms x real
F=@(x) (x>=-0.5 & x<=0.5);
N=32;T=4;t=linspace(-T/2,T/2,N);w0=2*pi/T;
X=F(t);
Y=fft(X);
subplot(1,2,1);stem(abs(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itAmplitud');
subplot(1,2,2);stem(angle(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itDesfase \theta');
disp('pulsa una tecla para el siguiente apartado');pause
close all
%
% Ejercicio 2
N=128;T=4;t=linspace(-T/2,T/2,N);w0=2*pi/T;
X=F(t);
Y=fft(X);
subplot(1,2,1);stem(abs(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itAmplitud');
subplot(1,2,2);stem(angle(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itDesfase \theta');
disp('pulsa una tecla para el siguiente apartado');pause
close all
%
N=128;T=20;t=linspace(-T/2,T/2,N);w0=2*pi/T;
X=F(t);
Y=fft(X);
subplot(1,2,1);stem(abs(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itAmplitud');
subplot(1,2,2);stem(angle(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itDesfase \theta');
disp('pulsa una tecla para el siguiente apartado');pause
close all
%
N=128;T=40;t=linspace(-T/2,T/2,N);w0=2*pi/T;
X=F(t);
Y=fft(X);
subplot(1,2,1);stem(abs(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itAmplitud');
subplot(1,2,2);stem(angle(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itDesfase \theta');
disp('pulsa una tecla para el siguiente apartado');pause
close all
%
%% Ejemplo
clear all
close all
T=10;w=2*pi/T;N=8;
F=@(t) 3*sin(w*t)+2*cos(4*w*t)-0.5*cos(5*w*t)
%F=inline('3*sin(w*t)+2*cos(4*w*t)-0.5*cos(5*w*t)') MAL, crea F(t,w)
%F=inline(eval('3*sin(w*t)+2*cos(4*w*t)-0.5*cos(5*w*t)'))
C=apxfourier(F,T,N);
subplot(1,2,1);stem(sqrt(sum(C.^2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('Termino k-esimo');ylabel('\itAmplitud c_k');axis([0,N,0,3.5]);
subplot(1,2,2);stem(atan2(C(2,:),C(1,:)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('Termino k-esimo');ylabel('\itDesfase \theta');axis([0,N,-3.5,3.5]);
disp('pulsa una tecla para el siguiente apartado');pause
close all
%
N=16;Ts=linspace(0,10,N);
fs=N/10;
X=F(Ts);
Y=fft(X);
subplot(2,1,1);plot(Ts,X,'Marker','.');title('Señal muestreada');
xlabel('t');ylabel('X(t)');
subplot(2,2,3);stem(abs(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itAmplitud');
subplot(2,2,4);stem(angle(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
xlabel('w');ylabel('\itDesfase \theta');
disp('pulsa una tecla para el siguiente apartado');pause
close all
%
% Ejercicio 3
%
N=[16,32,64,128];
for i=1:length(N);
    Ts=linspace(0,T,N(i));fs=N(i)/T;
    X=F(Ts);
    Y=fft(X);
    figure(i);
    subplot(2,1,1);plot(Ts,X,'Marker','.');title(['Señal muestreada a ',num2str(fs,'%g'),' Hz']);
    xlabel('t');ylabel('X(t)');
    subplot(2,2,3);stem(abs(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
    xlabel('w');ylabel('\itAmplitud');
    subplot(2,2,4);stem(angle(Y(1:N/2)),'Marker','.');title('Dominio de la Frecuencia');
    xlabel('w');ylabel('\itDesfase \theta');
end
disp('pulsa una tecla para el siguiente apartado');pause
close all
%
%% Filtros
% 
%
YY=Y;
k=5;YY(k)=0;YY(end-k+1)=0;
XX=ifft(YY);
plot(Ts,[X;XX]);
legend('original','filtrada');