function [a0,A,B] = fourier(l,n,f) //wszystko bezposrednio z wzorow ( http://www.dspguide.com/graphics/E_13_5.gif ), uzycie ponizej
    a0 = 1/l*intg(-l,l,f,1e-2)
    for i=1:n
        function b=f1(x)
            b = f(x)*cos(i*%pi*x/l)
        endfunction
        function c=f2(x)
            c = f(x)*sin(i*%pi*x/l)
        endfunction
        A(i) = 1/l*intg(-l,l,f1,1e-2)
        B(i) = 1/l*intg(-l,l,f2,1e-2)
    end
endfunction

function result = prostokatna(x) //definujemy przykladowa funkcje
    if x-floor(x)<0.5 then
        result = 1
    else
        result = 0
    end
endfunction


function result = trojkatna(x) //definujemy przykladowa funkcje
    if x - floor(x)<0.5 then
        result = x - floor(x)
    else
        result = 1-x + floor(x)
    end
endfunction

function result = piloksztaltna(x) //definujemy przykladowa funkcje
    result = x - floor(x)
endfunction

f = prostokatna //wybor funkcji

l = 0.5 //polowa okresu funkcji
t = 0:0.02:3

N = 2 //do ilu wspolczynnikow rozwijamy, więcej -> dokladniejsze odwzorowanie pierwotnej funkcji

[a0,A,B] = fourier(l,N,f)
subplot(311)
plot(t,B'*sin([1:N]'*t*%pi/l)+A'*cos([1:N]'*t*%pi/l)+a0/2,t,f) //przemnazanie i sumowanie skladowych w jednej linijce

N = 6

[a0,A,B] = fourier(l,N,f)
subplot(312)
plot(t,B'*sin([1:N]'*t*%pi/l)+A'*cos([1:N]'*t*%pi/l)+a0/2,t,f)

N = 10

[a0,A,B] = fourier(l,N,f)
subplot(313)
plot(t,B'*sin([1:N]'*t*%pi/l)+A'*cos([1:N]'*t*%pi/l)+a0/2,t,f)