function demosampling4
% Demo yksittäisen kosinisignaalin käytöksestä,
% kun sen taajuutta kasvatetaan. Näkymät aika- ja taajuustasossa.
% Demoa voi kuunnella!!!
%
% IKKUNA No. 1:
% 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.
%
% IKKUNA No. 2:
% Aikatason kuva. Sinisellä alkuperäinen kosinisignaali.
% Punaisilla palloilla kuvattu näytteistetty lukujono.
% Jos taajuus yli fs/2, niin sitten näytetään hitain
% kosini, joka sopii samaan näytteistettyyn lukujonoon.
% Eli juuri se taajuus, joka näkyy spektrissä 0 .. fs/2 välissä.
%
% Paina <enter> päästäksesi eteenpäin. Tällöin kuulet myös kyseisen
% kosinikomponentin ää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=5500 Hz.
% Koska f1 > fs/2, tapahtuu vierastuminen:
%   x(t) = cos(2 pi 5500 t) 
%   x[n] = x(nT_s) = x(n/fs)               %  t <- n/fs
%        = cos(2 pi 5500/10000 n) 
%        = cos(2 pi 5500/10000 n - 2 pi 10000/10000 n)  % väh. 2pi-monikerta
%        = cos(2 pi -4500/10000 n)
%        = cos(2 pi 4500/10000 n)          % kosini on parillinen
%   xr(t) = cos(2 pi 4500 t)               %  n <- t fs
% 
% Eli huomataan, että liian vähän näytteistetystä signaalista x(t)
% saadun diskreettiaikaisen lukujono x[n] sopii MYÖS 
% rekonstruoituun/palautettuun jatkuva-aikaiseen signaaliin xr(t),
% joka onkin ERI kuin alkuperäinen x(t). Tässä tapauksessa,
% kun näytteenottotaajuus on 10000 Hz ja signaalin x(t) kosinin
% taajuus 5500 Hz, niin se vierastuu koska 5500 > 10000/2,
% ja palautetun signaalin xr(t) taajuus on fs-f1=10000-5500=4500.
%
% Jos x(t):llä f1=11000 Hz, niin palautettu taajuus olisi
% f1-fs = 1000 Hz. Jos x(t):llä f1=33000 Hz, niin palautettu
% f1-3*fs = 3000 Hz. Jos x(t):llä f1 = 37000 Hz, niin palautettu
% 4*fs-f1 = 3000 Hz, jne.
%
% Jatkuva-aikaisen signaalin näytteenotto voidaan taajuustasossa
% siis ajatella JOKO
%   (1) KOPIOIMISENA jokaisen 2pi:n välein eli jokainen
%       näytteenottotaajuuden välein, spektrit summataan yhteen
%       jokaisessa kohdassa
%   (2) PEILAAMISENA jokaisen pi:n ympäri eli jokaisen 
%       puolen näytteenottotaajuuden ympäri, spektrit summataan
%       yhteen jokaisessa kohdassa. *)
% 
% Demossa kosinin taajuus muuttuu f0 Hz:stä fi Hz:n askelin 
% 1.2 fs Hz:iin asti. Koodissa f0 = 500, fi = 500 ja fs = 10000.
%
% 
% Jukka Parviainen, T-61.140, 16.4.2003 / 1.4.2004 / 27.4.2004 
% Informaatiotekniikan laboratorio, TKK.
%
% *) VIRHEITÄ:
%  Spektripiikit on aina skaalattu ykköseen, mikä ei aivan pidä
%  paikkaansa. Jos et tiedä miksi, älä huolestu vaan mene eteenpäin.

% ALUSTUS
%
% Alustetaan kaksi ikkunaa, Nro 1 on taajuusikkuna ja Nro 2
% on aikatason ikkuna. Lisäksi siirretään hieman ikkunaa Nro 2
% oikealle; tämä siksi, jos ikkunat ensimmäistä kertaa demottaessa
% ovat toistensa päällä, niin käyttäjä tajuaa, että toinenkin 
% ikkuna avautui.
figure(1); % taajuusikkuna
clf;
figure(2);
clf;
ppp = get(gcf,'Position');
set(gcf,'Position',ppp+[30 30 0 0]);  % siirretään
figure(1); % aktivoidaan taajuustason ikkuna

% näytteenottotaajuudet
Fs  = 10000;   % näytteenottotaajuus
fsc = 200000;  % simuloidaan jatkuva-aikaista

% ensimmäinen taajuus ja kuinka paljon taajuutta kasvatetaan
f0  = 500;
fi  = 500;

% aika-akselit
t   = [0 : 1/Fs : 0.3];   % viimeinen arvo määrää piippauksen pituuden
t   = t(1:1024);          % ...mutta 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
%
% TÄRKEÄT VÄRIT OVAT RIVIT 3 4 ja 5. rivi. YHTEENSÄ 7x3 matriisi
col2 = [[1 0 1];
    [0.2 0.3 0];
    [0 1 0.3];
    [0 0 1];  
    [.5 .5 .5];
    [.7 .2 .3];
    [1 1 0]];
    %[0 1 0.3]    % muita värejä?!
%col = ['b', 'b','b','b','b', 'b','b'];
%col = ['b', 'r','g','b','r', 'g','b'];
%col = ['g', 'r', 'c', 'b', 'm', 'k', 'y'];

% Kolmiospektrin väri tulee myös ylläolevista.
kolmvari = col2(4,:);

% Diskreetin sekvenssin piirtoväri ikkunassa Nro 2
dvari = [0 0 0]; 

% VARSINAINEN  TYÖ:
%
% For-looppi jokaista kosinitaajuutta kohti, lopussa
% on "pause", jolla ohjelman ajo keskeytetään siksi,
% kunnes käyttäjä painaa jotain nappulaa.
for m = [f0:fi:1.2*Fs]    
    f1 = m;
    y   = cos(2*pi*f1*t );   % "todellinen"
    yc  = cos(2*pi*f1*tc);   % "simuloitu" kuvaajia varten

    % lasketaan Fourier-muunnokset ja taajuusakselit
    
    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));

    % kolmiospektriä varten kohta:

    kolm = find(abs(ycF) > 0.8);
    kolm = round([-kolm(1) kolm(1)]+(length(ycF)/2));

    % PIIRROKSET

    % piirretään jatkuva-aikainen ikkunan Nro 1 yläosaan:
    subplot(311); cla;
    p1 = plot(wc, fftshift(ycF));
    set(p1, 'Linewidth', 3);
    grid on;
    axis([-2.2*Fs 2.2*Fs -0.2 1.1]);
    title(['Jatkuvan kosinin f = ' num2str(f1) ' Hz spektri. Skaala -2.2fs .. 2.2fs']);

    % piirretään kolmioaalto
    hold on;
    pk = plot([wc(kolm(1)) 0 wc(kolm(2))],[0 1 0]);
    set(pk, 'LineStyle', '--', 'Color', kolmvari);

    % piirretään diskreettiaikainen ikkunan Nro 1 keskelle:
    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);
    
    % piirretään varsinaiset komponentit
    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);

    % piirretään loput monikerrat, joista tapahtuu vierastumista
    %  kun f_high > fs/2
    for k=[-3:-1 1:3]
        % spektripiikit
        p2 = plot(wwj+k*Fs, fftshift(ycF));
        set(p2, 'Linewidth', 3, 'Color', col2(k+4,:));
        % nuolet alapuolelle
        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);
        % lisätään kolmiospektrit
        pk = plot([wc(kolm(1)) 0 wc(kolm(2))]+k*Fs,[0 1 0]);
        set(pk, 'LineStyle', '--', 'Color', col2(k+4,:));
    end;
    % HUOM! Nuolet kuvastavat näytteenottotaajuuden monikertoja
    % ne eivät ole mitään signaalikomponentteja
    
    % piirretään kolmioaalto
    hold on;
    pk = plot([wc(kolm(1)) 0 wc(kolm(2))],[0 1 0]);
    set(pk, 'LineStyle', '--', 'Color', kolmvari);

    title(['Diskreetin kosinin f = ' num2str(f1) ' Hz spektri, kun fs = ' num2str(Fs) ' Hz']);
    grid on;
    axis([-2.2*Fs 2.2*Fs -0.2 1.1]);

    % Ikkunan Nro 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]);

    % Tässä välissä soitetaan kosinin ääni ja jäädään
    % odottelemaan, että piirrokset ehtivät valmiiksi.
    soundsc(y, Fs);
    

    % AKTIVOIDAAN AIKATASON IKKUNA 
    figure(2); 
    clf;
    
    p1 = plot(tc, yc); 
    set(p1, 'LineWidth', 2); % jatkuva signaali
    hold on;
    
    vari = floor(f1 / (Fs/2));
    switch vari
    case 1,
        % Ollaan kohdassa fs/2 .. fs, eli vierastuu eka kertaa
        tmp = cos(2*pi*(Fs-f1)*tc);
        p3 = plot(tc, tmp);
        set(p3, 'Color', col2(5,:), 'LineWidth', 3);
        p3 = plot(tc, tmp);
        set(p3, 'Color', [0 0 0]);
    legend({'Alkuperainen','Palautettu'});
    case 2,
        % Ollaan kohdassa fs .. 3fs/2, eli vierastuu toka kertaa
        tmp = cos(2*pi*(f1-Fs)*tc);
        p3 = plot(tc, tmp);
        set(p3, 'Color', col2(3,:), 'LineWidth', 3);
        p3 = plot(tc, tmp);
        set(p3, 'Color', [0 0 0]);
    legend({'Alkuperainen','Palautettu'});
    end;
    
    tmpx = [0 : 1/Fs : 0.01];
    tmpy = (-1.5)*ones(1, length(tmpx));
    p4 = plot(tmpx,tmpy,'d');
    set(p4, 'MarkerSize', 16, 'LineWidth', 2, 'Color', dvari);
    p4 = plot(tmpx,tmpy,'-');
    set(p4,  'LineWidth', 2, 'Color', dvari);
    
    t4 = text(0.1/1E3, -1.3, ['fs = ' num2str(Fs) ' Hz, Ts = 1/fs = ' num2str(1/Fs) ' s']);
    set(t4, 'FontSize', 20, 'FontWeight', 'Bold', 'Color', dvari);
    
    p4 = stem(t, y, 'o');
    set(p4, 'MarkerSize', 16, 'LineWidth', 3, 'Color', dvari);

    title(['Aikatason jatkuva kosini f = ' num2str(f1) ' Hz']);
    xlabel('aika (s)');
    
    grid on;
    axis([0 0.0011 -1.8 1.2]);
    shg;
    
    % PALATAAN TAKAISIN IKKUNAAN Nro 1
    figure(1);
    shg;
    
    % Jäädään odottelemaan käyttäjän näppäimenlyöntiä...a
    pause;
end;