function fourier2130aa(typeKurve)

% Genererer en funksjon som er en ren sinus eller cosinus, over et intervall 
% på 512 pkt dersom man kaller funksjonen med 1 som innparameter. For 2 som 
% input, får man en sagtann osv. Se nedenfor. 

% Dernest plottes funksjonen vi har generert, det foretas en diskret Fourier
% transformasjon og resultatet vises. Deretter gjennomføres en omvendt
% Fourier transformasjon, og man sammenligner resultatet man da får med det
% opprinnelige. Sist beregnes differansen mellom opprinnelig og 2x
% transformert signal for å kunne se hvor stor forskjellen egentlig er.

% Programmet er skrevet for bruk i FYS2130 våren 2008 av Arnt Inge Vistnes

N = 512;
x = linspace(0,2*pi*(N-1)/N,N);
f_x = zeros(1,N);

switch typeKurve
    case 1  % Sin og cos kurver
        x = linspace(0,2*pi*(N-1)/N,N);
        f_x = cos(10.4*x);
    case 2  % Sagtann, en eneste tann
        for i = 1:N
            f_x(i)=i*2/N;
        end;
    case 3  % En firkantpuls (en eneste puls)
        for i = 1:N/2
            f_x(i) = +1.0;
            f_x(i+N/2) = -1.0;
        end;
    case 4  % Åtte firkant pulser innen det området vi betrakter
        for i = 1:32:N
            for j = 1:16
                f_x((i-1)+j) = +1.0;
                f_x((i-1)+16+j) = -1.0;
            end;
        end;
    case 5  % Firkantpuls med nest høyeste "frekvens" synkron med sampling
    % (to punkt pr opp og to punkt pr ned).
        for i = 1:N
            if (mod((i-1),4) == 0) 
                f_x(i) = +1; 
                f_x(i+1) = +1;
                f_x(i+2) = -1; 
                f_x(i+3) = -1;
            end;
        end;
    case 6  % Firkantpuls med høyeste "frekvens" som kan håndteres ok
    % (bare ett punkt pr opp og ett punkt ned).
        for i = 1:N
            if (mod(i,2) == 0) f_x(i) = -1;
            else f_x(i) = +1;
            end;
        end;
    case 7  % Deltapuls
        f_x = zeros(1,N);
        f_x(200) = 1;
    case 8  % Bølgepakke
        f_x = zeros(1,N);
        delta = 8; % Varier denne f.eks. mellom 4 og 200
        for i = 1:N
            f_x(i) = sin(2*pi*(i-1)/4)*exp(-(i-N/2)*(i-N/2)/(delta*delta));
        end;
    case 9  % Amplitudemodulert bærebølge
        f_x = zeros(1,N);
        modfrekv = 12; 
        for i = 1:N
            f_x(i) = sin(2*pi*(i-1)/4)*(1.0-0.3*sin(2*pi*(i-1)*modfrekv/N));
        end;
end;

% Plotter opprinnelig funksjon
plot(x,f_x,'.-k');
%xlim([-0.2,6.5]);
%ylim([-1.2,1.2]);
title('Opprinnelig funksjon');
figure;

an = zeros(1,N);
bn = zeros(1,N);

for j=1:N
    sum_a = 0;
    sum_b = 0;
    delta_omega = (j-1)*2*pi/N;
    for i=1:N
        sum_a = sum_a + f_x(i)*cos(delta_omega*(i-1));
        sum_b = sum_b + f_x(i)*sin(delta_omega*(i-1));
    end;
    an(j) = sum_a.*(1/N);
    bn(j) = sum_b.*(1/N);
end;
plot(x,an,'.-r');
xlim([-0.2,6.5]);
hold on;
plot(x,bn,'.-b');
title('Fourier-komponenter (cos-rød, sin-blå)');

figure;

% Tegner opp tilbaketransformasjonen:

rekonstruksjon = zeros(1,N);
for j = 1:N
    sum = 0;
    delta_omega = (j-1)*2*pi/N;
    for i = 1:N
        sum = sum + an(i)*cos(delta_omega*(i-1)) + bn(i)*sin(delta_omega*(i-1));
    end;
    rekonstruksjon(j) = sum;
end;

plot(x,rekonstruksjon,'.-g');
%xlim([-0.2,6.5]);
%ylim([-1.2,1.2]);
title('Rekonstruert funksjon');

% Plotter så differansen mellom startfunksjon og den rekonstruerte

figure;
plot(x,f_x-rekonstruksjon,'.-k');
xlim([-0.2,6.5]);
title('Avvik mellom rekonstruksjon og opprinnelig funksjon');

    
