function [erosum] = triangpulse(M)
% TRIANGPULSE approksimoi kolmioaaltoa 2M+1:llä 
% Fourier-kertoimien a_(-M) .. a_(+M) muodostamilla termeillä. 
% Aloita M=0 ja kasvata tämän jälkeen vähitellen suuremmaksi.
%
% [erosum] = triangpulse(M)
%    M        (1x1)  Fourier-kertoimien määrä 2M+1
%    erosum   (1x1)  neliöllinen virhesumma (x_1(t)-x_2(t))^2
%
% Visualisointi: ylimmässä kuvassa ovat kolmioaalto ja siitä
% (2M+1):llä kertoimella approksoitu käyrä. Keskimmäisessä
% kuvaajassa on niiden välinen erotus. Alimmassa on Fourier-
% kertoimien itseisarvot. Sininen viittaa oikeaan kolmioaaltoon,
% punainen väri (2M+1):llä kertoimella approksimoituun kolmioaaltoon,
% ja musta väri niiden erotukseen.

% Muuttuja Ak viittaa nyt kolmioaallon kertoimiin, erityisesti
% nyt niihin (2M+1) kertoimeen. Koodin lopussa on muuttuja Ck,
% joka on käytännössä sama, mutta (2(M+10)+1) kertoimeen laajemmin.
%
% x1 ja t1 viittavat kolmioaaltoon "symbolisesti" ja
% x3 ja t3 samaan kolmioaaltoon "laskennallisesti", kun lasketaan
% neliöllistä virhesummaa.
%
% x2 ja t2 viittaavat (2M+1):llä kertoimella luotuun Fourier-
% sarjaesitykseen.
%
% Jukka Parviainen, 2004, 17.2.2005

if (M>20)
  error('Muuttujan M maksimiarvo on rajoitettu 20:een.');
end;

disp(['Muuttuja M = ' num2str(M) ' eli ' num2str(2*M+1) ' kerrointa']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fourier-kertoimet Ak

for k= -M : M                    % Fourier-kertoimet a_k
  if (mod(k,2)==1)
    Ak(k+M+1) = -2/((k*pi).^2);  % Ak(1) == a_(-M), jne...
  elseif k==0
    Ak(M+1) = 0.5;               % Ak(M+1) == a_0
  else
    Ak(k+M+1) = 0;
  end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fourier-kertoimet Ck

C = M+10;
for k= -C : C
  if (mod(k,2)==1)
    Ck(k+C+1) = -2/((k*pi).^2);  
  elseif k==0
    Ck(C+1) = 0.5;
  else
    Ck(k+C+1) = 0;
  end;
end;
minCk = 1.1*min(abs(Ck(:)));
maxCk = 1.1*max(abs(Ck(:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Aaltomuodon x- ja t-vektorit

t1 = [-2 -1 0 1 2 3 4 5];        % kolmioaallon x-aks. arvot
x1 = [ 0  1 0 1 0 1 0 1];        % kolmioaalton y-aks. arvot

t2 = linspace(-2, 5, 14000);     % aika-akseli [-2 .. 5] approksivoivalle funkt.
x2 = 0;
T  = 2;                          % T_0, Omega = 2pi f = 2pi/T_0
for k= -M : M                    % luodaan x2(t)
  x2 = x2 + Ak(k+M+1)*exp(j*k*2*pi/T*t2);
end;
x2 = real(x2);                   % ignore complex artifacts

x3up    = linspace(0, 1-(1/2000), 2000);  % eron laskemista varten kolmiosign.
x3down  = linspace(1, (1/2000), 2000);
x3      = [x3up x3down x3up x3down x3up x3down x3up];
ero     = x3 - x2;
erosumT = sum(abs(ero).^2)/(t1(end)-t1(1));
erosum  = (1/T)*erosumT;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualisointi

figure(1); clf;

subplot(311);
plot(t1, x1, 'b', t2, x2, 'r');
title(['Signaali x_1(t) ja Fourier-sarjalla ' ...
  'approksimoitu x_2(t), ' num2str(2*M+1) ' F-kerrointa']);
legend({'alkuper. x_1(t)','F-sarjan appr. x_2(t)'});
grid on; 
axis([-2 5 -0.2 1.2]);     % zoomataan

subplot(312);
plot(t2, ero, 'k');
title(['Erotus e(t) = x_1(t) - x_2(t); Neliosumma: (1/T) \int_T |e(t)|^2 dt = ' num2str(erosum)]);
legend({'x_1(t) - x_2(t)'});
grid on; 
axis([-2 5 1.2*min(ero(:)) 1.2*max(ero(:))]);     % zoomataan

subplot(313);
dummy1 = plot([0.001 0.0012],[0.1 0.1002]);
set(dummy1,'Color','b','LineWidth',3);
hold on;
dummy2 = plot([0.001 0.0012],[0.1 0.1002]);
set(dummy2,'Color','r','LineWidth',3);
legend({'Kolmioaalto |a_k_1|','Appr. |a_k_2|'});

for k = 1:2*C+1
   l1 = line([k-1-C k-1-C],[0 abs(Ck(k))]);
   set(l1,'LineWidth',3);
end;
for k = 1:2*M+1
   l1 = line([k-1-M k-1-M]+0.15,[0 abs(Ak(k))]);
   set(l1,'Color',[1 0 0],'LineWidth',3);
end;
title(['Fourier-kertoimet |a_k_1| ja |a_k_2|']);
set(gca,'XTick',[-M:M]);
axis([-C-1 C+1 minCk maxCk]);
grid on;

shg;  % bring front
