function demosampling2(fA, fB, fC, fI, ap, ar, Fs)
% DEMOSAMPLING2 - demo kolmen kosinin muodostamasta signaalista 
% näytteistyksessä.
%
% Joko ilman muuttujia tai muuttujat seuraavasti:
%    Kosini #1   1  . cos(2pi fA+(k . fI) t)    % default fA = 400
%    Kosini #2   ap . cos(2pi fB+(k . fI) t)    % default ap = .7, fB = 900
%    Kosini #3   ar . cos(2pi fC+(k . fI) t)    % default ar = .4, fC = 1750
%    Intervalli fI = 600Hz.
%    Näytteenottotaajuus fs = 10 kHz.
%
% Demo useamman kosinisignaalin käytöksestä,
% kun sen taajuutta kasvatetaan. Demoa voi kuunnella!!!
% Yksittäistä kosinia varten DEMOSAMPLING1.
%
% Ylin ikkuna: jatkuva-aikaisen signaalin x(t) spektri. Yksi kosini
%   tuottaa kaksi piikkiä spektriin symmetrisesti y-akselin 
%   molemmin puolin. Spektri ei ole jaksollinen.
% Keski-ikkuna: diskreettiaikaisen signaalin/sekvenssin x[n] 
%   spektri -2.2 fs .. 2.2 fs. Spektri on jaksollinen. 
%   Näytteenotto voidaan ajatella  jatkuva-aikaisen signaalin
%   spektrin kopiointina näytteenottotaajuuden välein (2pi) ja
%   summattuna komponentit yhteen. Väli 0 .. fs/2 on tummennettu,
%   koska sille kaistalle osuvat kosinit voidaan  kuulla rekonstruktiossa.
% Alin ikkuna: diskreettiaikaisen signaalin spektri 0 .. 0.5 fs
%   Siis zoomaus keski-ikkunasta, se osa signaalia, joka voidaan
%   havaita. Kun diskreettiaikainen sekvenssi x[n] palautetaan
%   jatkuva-aikaiseksi xr(t), vain tässä alueessa olevat komponentit
%   (ja niiden vastinkappaleet vastaavalla negatiivisella taajuudella)
%   voidaan havaita.
%
% Paina <enter> päästäksesi eteenpäin. Tällöin kuulet myös kyseisten
% kosinien muodostaman äänen. 
%
% Ohjelmalla demonstroidaan vierastumisilmiötä (aliasing).
% Jos kosinin taajuus on yli fs/2:n, se vierastuu matalemmalle taajuudelle.
% Matemaattisesti tämän voi ajatella esimerkillä, jossa fs=10000 Hz,
% ja alkuperäinen jatkuva-aikaisen kosinin taajuus f1=5200 Hz.
% Koska f1 > fs/2, tapahtuu vierastuminen:
%   x(t) = cos(2 pi 5200 t) 
%   x[n] = x(nT_s) = x(n/fs)               %  t <- n/fs
%        = cos(2 pi 5200/10000 n) 
%        = cos(2 pi 5200/10000 n - 2 pi 10000/10000 n)  % väh. 2pi-monikerta
%        = cos(2 pi -4800/10000 n)
%        = cos(2 pi 4800/10000 n)          % kosini on parillinen
%   xr(t) = cos(2 pi 4800 t)               %  n <- t fs
% 
%
% 
% Jukka Parviainen, T-61.140, 16.4.2003 / 1.4.2004
% Informaatiotekniikan laboratorio, TKK.

% ensimmäinen taajuus ja kuinka paljon taajuutta kasvatetaan
if nargin == 0
    f0  = 400;
    fi  = 600;
    fp  = 500; 
    ap  = 0.7;
    fr  = 1350; 
    ar  = 0.4;
    Fs  = 10000;   % näytteenottotaajuus
elseif (nargin == 7)
    f0 = fA;
    fi = fI;
    fp = fB-fA;
    fr = fC-fA;
else
    error('Vaara maara muuttujia, joko 0 tai 7, "help demosampling2"');
end;

clf;

% ylimääräinen KORKEA näytteenottotaajuus kuvaajia varten
fsc = 200000;
% aika-akselit
t   = [0 : 1023]/Fs;          % Otetaankin vain 1024 näytettä = 1024/Fs sek.
tc  = [0 : 1/fsc : 0.1];
tc  = tc(1:2^14);         % oikeita vähän reilummin...

% spektrin värjäys - ongelmia dataprojektin kanssa :(
% käytetään joko col tai col2
col2 = [[1 0 1];
    [0.2 0.3 0];
    [1 1 0];
    [0 0 1];
    [.5 .5 .5];
    [.7 .2 .3];
    [0 1 0.3]];
%col = ['b', 'b','b','b','b', 'b','b'];
%col = ['b', 'r','g','b','r', 'g','b'];
%col = ['g', 'r', 'c', 'b', 'm', 'k', 'y'];


% Varsinainen työ:
for m = [f0:fi:1.2*Fs]    
    f1 = m;
    y   = cos(2*pi*f1*t ) + ap*cos(2*pi*(f1+fp)*t ) + ar*cos(2*pi*(f1+fr)*t );
    yc  = cos(2*pi*f1*tc) + ap*cos(2*pi*(f1+fp)*tc) + ar*cos(2*pi*(f1+fr)*tc);

    % spektrit
    
    ycF = abs(fft(hamming(length(yc))'.*yc));
    ycF = ycF/max(ycF);
    wc  = linspace(-fsc/2, fsc/2, length(ycF));

    yF = abs(fft(hamming(length(y))'.*y));
    yF = yF/max(yF);
    w = linspace(0, Fs, length(yF));

    % piirretään jatkuva-aikainen
    subplot(311);
    p1 = plot(wc, fftshift(ycF));
    set(p1, 'Linewidth', 3);
    grid on;
    axis([-2.2*Fs 2.2*Fs -0.2 1.1]);
    if ((abs(ar)>0) | (abs(ap)>0))
       title(['JATKUVA-AIKAISEN signaalin cos(2\pi \cdot ' num2str(f1) ' \cdot t) + ' num2str(ap) ' cos(2\pi \cdot ' num2str(f1+fp) ' \cdot t) + ' num2str(ar) ' cos(2\pi \cdot ' num2str(f1+fr) ' \cdot t) SPEKTRI']);
    else
       title(['JATKUVA-AIKAISEN signaalin cos(2\pi \cdot ' num2str(f1) ' \cdot t) SPEKTRI']);
    end;

    % piirretään diskreettiaikainen
    subplot(312); cla;
    % kiinnostava kaista 0 .. pi tummenetaan
    ff = fill([0 Fs/2 Fs/2 0],[0 0 1.1 1.1],[0.8 0.8 0.8]);
    hold on;
    % li = line([0 Fs/2], [-0.18 -0.18]);  % tämä jana taitaa sekoittaa
    % set(li, 'LineWidth', 3);
    
    k=0;
    wwj = linspace(-fsc/2 , fsc/2 , length(ycF));
    p2 = plot(wwj, fftshift(ycF));
    set(p2, 'Linewidth', 3, 'Color', col2(k+4,:));

    li = line([0 fi -fi; 0 0 0], [-0.1 -0.15 -0.15; -0.2 -0.1 -0.1]);
    set(li, 'Color', col2(4,:),'LineWidth',5);

    for k=[-3:-1 1:3]
        p2 = plot(wwj+k*Fs, fftshift(ycF));
        set(p2, 'Linewidth', 3, 'Color', col2(k+4,:));
        li = line([k*Fs k*Fs+fi k*Fs-fi; k*Fs k*Fs k*Fs],[-0.1 -0.15 -0.15; -0.2 -0.1 -0.1]);
        set(li, 'Color', col2(k+4,:),'LineWidth',5);
    end;
    % HUOM! Nuolet kuvastavat näytteenottotaajuuden monikertoja
    % ne eivät ole mitään signaalikomponentteja
    
    title(['Diskreetin sekvenssin spektri, kun fs = ' num2str(Fs) ' Hz']);
    grid on;
    axis([-2.2*Fs 2.2*Fs -0.2 1.1]);

    % Alimpaan kuvaajaan tulee kiinnotava kaista 0 .. pi eli 0 .. fs/2
    subplot(313); 
    p3 = plot(w, yF);
    set(p3, 'Linewidth', 3);
    vari = floor(f1 / (Fs/2));
    switch vari
    case 0,
        set(p3, 'Color', col2(4,:));
    case 1,
        set(p3, 'Color', col2(5,:));
    case 2,
        set(p3, 'Color', col2(3,:));
    otherwise,
        set(p3, 'Color', col2(4,:));
    end;
    axis([0 Fs/2 -0.2 1.1]);
    grid on;
    title(['Diskreettiaikainen spektri 0 .. ' num2str(Fs/2) ' Hz - palautetun signaalin spektri']);
    set(gca, 'Color', [.8 .8 .8]);

    shg

    soundsc(y, Fs);
    pause;
end;

