% 1/4 kitöltési tényezõjû négyszögjel Fourier sorfejtése
% Budai Tamás, Unger Tamás
% 2012

clear;
tic;

%Sorfejtés lépéseinek számát definiáljuk
n = 1000;

%Konstansok definiálása
T = 10;
R = 1;
C = 5;
omega = (2*pi) / (T);

%Idõváltozó definiálása
t  = linspace(0,20,2000);

%A gerjesztõjel definiálása
s = 5 * (heaviside(t-(0.25*T)) - heaviside(t-T) + heaviside(t-(1.25*T)) - heaviside(t-(2*T)));

%A gerjesztõjel kirajzolása 
figure(1);
plot(t,s);
hold all;

%S0 DC-komponens elõre definiálása
S0 = 3.75;

% SkA és SkB elõre lefoglalása (optimalizáció)
SkA = zeros(1,n);
SkB = zeros(1,n);
% A for iteráció következik, ahol k a futóváltozó
for k = 1:n,
    SkA(1,k) = (10 / (k*2*pi)) * ( (sin(k*2*pi)) - (sin(k*(pi/2))));
    SkB(1,k) = (-10 / (k*2*pi)) * ( (cos(k*2*pi)) - (cos(k*(pi/2))));
end

%Sk komplex vektorok hosszainak számítása
Sk_length = zeros(1,k);
for k = 1:n,
    Sk_length(1,k) = sqrt((SkA(1,k)^2) + (SkB(1,k)^2));
end

%Sk komplex vektorok szögeinek számítása
Sk_phase = zeros(1,k);
for k = 1:n,
    Sk_phase(1,k) = -atan2( (SkB(1,k)) , (SkA(1,k)) );
end

%Átviteli függvény
Wk=zeros(1,k);
for k = 1:n,
    Wk(1,k) = ((k*omega*R*C)*1i) / (1+((k*omega*R*C)*1i));
end

Wkreal=zeros(1,k);
Wkimag=zeros(1,k);
phaseWk=zeros(1,k);
Wkabs=zeros(1,k);
for k = 1:n,
    Wkreal(1,k) = real(Wk(1,k));
    Wkimag(1,k) = imag(Wk(1,k));
    phaseWk(1,k) = atan2( Wkimag(1,k) , Wkreal(1,k));
    Wkabs(1,k) = abs(Wk(1,k));
end


% a gerjesztés egyes komponenseinek kiszámítása
sc = zeros(k,size(t,2));
for k =1:n,
    a = Sk_length(1,k)*cos(k*omega*t+Sk_phase(1,k));
    sc(k,:) = a;
end;

%A válaszjel komponenseinek kiszámítása
yc = zeros(k, size(t,2));
for k =1:n,
    b = ((Sk_length(1,k)) * (Wkabs(1,k))) * cos(k*omega*t + ((Sk_phase(1,k)) + (phaseWk(1,k))));
    yc(k,:) = b;
end

% a gerjesztés komponenseinek összegzése
st_f = S0;
for k=1:n       
    st_f=st_f+sc(k,:);
end;

%A válaszjel komponenseinek összegzése
yt_f = 0;
for k=1:n,
    yt_f = yt_f+yc(k,:);
end;

%plot(t,yt_f);
plot(t,st_f);

%for i=1:size(t,2)
%    if s(1,i)==0
%        y2(1,i)=0;
%    else
%        if(yt_f(1,i)<0)
%            y2(1,i)=0;
%        else
%           y2(1,i)=yt_f(1,i);
%        end;
%    end;
%end;
%plot(t,y2);
toc;
plot(t,yt_f);