%TEKNILLINEN KORKEAKOULU
%Tietoliikennelaboratorio
%S-72.060 Signaalit ja järjestelmät
%Seppo Saastamoinen 2.3.2005 ver. 1.0
%MAKRO: pcm2005.m
clear; close all; %Nollataan muistit ja kuvat

Teekkari1 = 'Teemu Teekkari 12345A';
Teekkari2 = 'Liisa Teekkari 54321B';

%PCM

%Tässä työssä tutustutaan AD/DA-muunnoksiin.
%Muunnettava signaali voi olla itse tuotettu tai omaa ääntä.
%Muunnettavan signaaliin voidaan summata häiriösignaali.
%Ennen AD-muunninta voidaan kytkeä päälle 10 kHz alipäästösuodatin,
%joka mukailee ideaalista mallia.
%Mahdollisen suodattimen jälkeen voidaan signaalia haluttaessa
%kompressoida A-lain mukaan.
%AD-muuntimen muunnosalueen rajat voidaan valita mieleisiksi.
%Myös muuntimen bittisyys voidaan valita halutuksi (1-N).
%Tietoverkkoa ei tässä makrossa ole, vaan AD-muuntimen lähtö kytketään
%sellaisenaan DA-muuntimen tuloon.
%Muunnettu signaali DA-muunnetaan ja muunnos on mahdollista
%alipäästösuodattaa ja dekompressoida.
%Laskennallisessa osiossa voidaan määrittää mm. kvantisointivirhe
%(-kohina) niin aika-alueessa kuin tehollisarvona.
%Makron lopuksi tulostetaan kaikkien mittauspisteiden
%arvot aika- ja taajuusalueessa.

%Listauksen selityksiä:
%-------------------- = osion raja
%Signaalit aika-alueessa = pienet kirjaimet, esim. in, xt jne.
%Signaalit taajuusalueessa = SUURET KIRJAIMET, esim. IN, Xf jne.

%Tässä makrossa on for-silmukoita ja ehtoja (if, elseif, else),
%joiden avulla saat kytkettyä jonkun ominaisuuden päälle/pois.

%-------------------- 
%Tulosignaaliosio

%Tulosignaali voi olla generoitu tai ääntä.
%Jos generoidaan signaali itse => aani = 0,
%jos käytetään omaa ääntä => aani = 1.
aani = 0; 

if aani == 0;
   
%Generoitu tulosignaali

Fs = 22050; %Näytevauhti, sopii äänikortille.
Ts = 1/Fs; %Näyteväli.

t = 0:Ts:1-Ts; %Näytteistettävän signaalin aikavektori, esimerkiksi 1 s.

f = 100; %Esimerkkisignaalin taajuus.

%Tähän alle tulee tulosignaalin kaava:
in = 1*cos(2*pi*f*t); %+0.2*cos(2*pi*15000*t);
%Esimerkkisignaali voi sisältää häiriön (15000 Hz).
in =in'; %Transponoidaan vektori

else
   
%Haetaan äänitetty signaali tulosignaaliksi;
[in,Fs] = wavread('c:\Temp\Sinkut\nimi.wav'); 

Ts = 1/Fs;

end

%Lasketaan arvoja jatkoa varten:
N = max(size(in)); %Tulosignaalin näytemäärä.
T = N*Ts; %Näytteistettävän signaalin kesto.
t = 0:Ts:T-Ts; %Äänisignaalin aikavektori
fo = 1/T; %Taajuuden askellusväli spektrilaskennassa.

INf = 2*abs(fft(in))/N; %Tulosignaalin 1-puoleinen spektri.

%--------------------
%Tulosignaalin alipäästösuodatus
%1 = päällä, 0 = ei suodatusta
suodatus = 1; 

%Alipäästösuodatin on tarpeen (pakollinen), jotteivät mahdolliset
%häiriösignaalit aiheuttaisi laskostumistuloksia taajuusalueessa.

A = 0.5; %Suodattimen vahvistus. 
B = 10000; %Suodattimen ylärajataajuus.
t1 = -T/20:Ts:T/20; %Aikavektori sinciä varten.
%Lähes ideaalisen alipäästösuodattimen impulssivaste on
h1t = A*(sin(pi*2*B*t1)./(pi*2*B*t1)); %Eli h1(t) = A*sinc(2Bt).
N1 = max(size(h1t)); %Lähes ideaalisen alipäästösuodattimen impulssivasteen näytemäärä.
fo1 = 1/(N1*Ts); %Taajuuden laskentaväli suodatinta koskien.
H1f = 2*abs(fft(h1t))/N1; %Lähes ideaalisen alipäästösuodattimen siirtofunktion itseisarvo.

if suodatus == 1

x11 = conv(h1t,in); %Suodattimen lähtö AD-muuntimelle.
N2 = max(size(x11)); %Suodatetun signaalin näytemäärä.
%Korjataan suodattimen lähtö samanaikaiseksi in-vektorin kanssa eli
x1 = x11(N1:N2); 
%Korjataan suodattimen lähtösignaalin amplitudi samaksi kuin tulosignaalin amplitudi eli
x1 = x1.*(max(in(2*N1:5*N1))/max(x1(2*N1:5*N1))); %Outo lauseke tarkoittaa sitä, että
%katsotaan maksimiarvot signaaleiden keskialueilta, jolloin konvoluutiosta johtuvat
%alku- ja loppupään heilahtelut eivät vaikuta.
%Tämä kompensoi suodattimen vahvistuksen/vaimennuksen pois.

else; %Ei suodatusta.
   x1 = in; N2 = N;
   
end
%Suodatettu tai suodattamaton tulosignaali taajuusalueessa, 1-puoleinen spektri.
X1f = 2*abs(fft(x1))/N; 
   
%--------------------
%Kompressointi

%HUOM! Jos kompressointi on päällä, tulee myös dekompression olla päällä!

%1 = päällä, 0 = pois päältä.
kompressointi = 0; 

if kompressointi == 1;
   
A = 87.6; %A-lain mukainen kerroin.

%Lasketaan x1:lle uudet arvot vektoriin comp.

comp = x1; %Kopioidaan aluksi x1:n arvot comp-vektoriin.

for i = 1:length(x1); %Kelataan koko signaali arvoittain alusta loppuun.
   if abs(x1(i,1)) < 1/A; %Pienten arvojen laskenta.
      comp(i,1) = x1(i,1)*A/(1+log(A));
   elseif x1(i,1) > 0 %Positiivisten suurten arvojen laskenta.
      comp(i,1) = (1+log(x1(i,1)*A))/(1+log(A));
   else %Negatiivisten suurten arvojen laskenta.
      comp(i,1) = -1*(1+log(abs(x1(i,1))*A))/(1+log(A));
   end
end

%Jos ei kompressointia:
else
comp = x1;

end

%Kompression lähtösignaali taajuusalueessa, 1-puoleinen spektri.
COMPf = 2*abs(fft(comp))/N; 

%--------------------
%AD-muunnos
%Näytealueen jänniteen yläraja UP,
%joka saa suurimman binaari- tai numeroarvon, esim. 1111B eli 15.
UP = 1; 
%Näytealueen alaraja,
%joka saa pienimmän binaari- tai numeroarvon, esim. 0000B eli 0.
LOW = -1; 

%AD-muuntimen bittisyys. Tätä voit muuttaa!
bitit = 4; 

kvantti = (UP-LOW)/(2^bitit-1); %Kvantin koko eli muunnoksen jänniteporras.

x2 = round((comp-LOW)./kvantti); %AD-muuntimen tulos desimaalisena.

%--------------------
%Tietoverkko

%Tässä osiossa voisi olla televerkkoa mallintava 300 - 3400 Hz
%kaistanpäästösuodatin. Nyt kuitenkin on verkko korvattu johdolla eli
%AD-muuntimen lähtö on suoraan kytketty DA-muuntimen tuloon

%--------------------
%DA-muunnos
out1 = LOW+x2.*kvantti; %DA-muuntimen lähtö

%DA-muunnettu signaali taajuusalueessa, 1-puoleinen spektri.
OUT1f = 2*abs(fft(out1))/N; 

%--------------------
%Dekompressio

%1 = dekompressio tehdään, 0 = ei tehdä
dekompressio = 0; 

if dekompressio == 1

%Lasketaan out1:lle uudet arvot vektoriin decomp.

decomp = out1; %Kopioidaan aluksi out1:n arvot comp-vektroriin.

for i = 1:max(size(out1)); %Kelataan koko signaali arvoittain alusta loppuun.
   if abs(out1(i,1)) < 1/(1+log(A)); %Pienten arvojen laskenta.
      decomp(i,1) = out1(i,1)*(1+log(A))/A;
   elseif out1(i,1) > 0 %Positiivisten suurten arvojen laskenta.
      decomp(i,1) = (1/A)*exp((1+log(A))*abs(out1(i,1))-1);
   else %Negatiivisten suurten arvojen laskenta.
      decomp(i,1) = -1*(1/A)*exp((1+log(A))*abs(out1(i,1))-1);
   end
end 

else; %Jos ei dekompressointia:

decomp = out1;

end

%Dekompressoitu signaali taajuusalueessa, 1-puoleinen spektri.
DECOMPf = 2*abs(fft(decomp))/N; 

%--------------------
%Lähtösignaalin suodatus

%PIDÄ TÄMÄ SUODATIN KOKO AJAN PÄÄLLÄ!

lahtosuodatus = 1; %1 = suodatus päällä, 0 = ei suodatusta

%Sama suodatin kuin tulopuolella eli
h1t = A*(sin(pi*2*B*t1)./(pi*2*B*t1)); %Eli h1(t) = A*sinc(2Bt).

if lahtosuodatus == 1

out = conv(h1t,decomp); %Alipäästösuodattimen lähtösignaali

N2 = max(size(out)); %Suodatetun signaalin näytemäärä.

%Korjataan suodattimen lähtö samanaikaiseksi suodattimen tulosignaalin kanssa eli
out2 = out(N1:N2); 

%Korjataan suodattimen lähtösignaalin amplitudi samaksi kuin tulosignaalin amplitudi eli
out2 = out2.*(max(decomp(2*N1:5*N1))/max(out2(2*N1:5*N1)));
%Outo lauseke tarkoittaa sitä, että katsotaan maksimiarvot signaaleiden keskialueilta,
%jolloin konvoluutiosta johtuvat alku- ja loppupään heilahtelut eivät vaikuta.
%Tämä kompensoi suodattimen vahvistuksen/vaimennuksen pois.

else %Jos ei suodatusta
out2 = decomp;

end

%Korjataan vielä signaalia, koska konvoluutiosta johtuen signaali loppupää on väärin aika-alueessa.
%Eli otetaan vektoriin ns. tervyt osuus signaalista.

out3 = out2(1:N2-2*N1);

N3 = max(size(out3)); %Korjatun signaalin näytemäärä
t3 = 0:Ts:N3*Ts-Ts; %Korjatun signaalin aikavektori.
T3 = N3*Ts; %Korjatun vektorin kesto.
f3 = 0:1/T3:Fs-1/T3; %Korjatun signaalin taajuusvektori

OUT3f = 2*abs(fft(out3))/N3; %Suodatettu tulosignaali taajuusalueessa, 1-puoleinen spektri.

%Korjataan vielä tulovektori vastaamaan pituudeltaan lähtövektoria eli

in3 = in(1:N2-2*N1);

%--------------------
%Laskutoimituksia
d = in3 - out3; %Kvantisointivirhe eli tulo - lähtö

virhe = std(d); %Kvantisointivirheen tehollisarvo eli keskihajonta.

%--------------------
%Tulostusosio
%Tulostetaan 3 arkkia (Figure) käyräkuvia

%Haetaan päiväys ja kellonaika tulostusta varten eli
h = fix(clock);

figure %1. arkki

subplot(3,2,1); plot(t,in);
axis([0 max(t) 1.2*min(in) 1.2*max(in)])
xlabel('Aika/s'); ylabel('Amplitudi/V')
title('Tulosignaali(t)'); grid on; zoom on;

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%VALITSE JOKU ALLAOLEVISTA SIGNAALEISTA TODEKSI, NIIN KAIKKIIN
%OTSIKOIHIN (3 KPL) TULOSTUVAT TEKSTIT AUTOMAATTISESTI OIKEIN!

%signaali = '1*cos(2*pi*100*t), kompanderi on käytössä, 4-bittinen muunnos'
%signaali = '1*cos(2*pi*100*t), kompanderi ei ole käytössä, 4-bittinen muunnos'
%signaali = '10e-3*cos(2*pi*100*t), kompanderi on käytössä, 4-bittinen muunnos'
signaali = '10e-3*cos(2*pi*100*t), kompanderi ei ole käytössä, 4-bittinen muunnos'
%signaali = '10e-3*cos(2*pi*100*t), kompanderi on käytössä, 8-bittinen muunnos'
%signaali = '10e-3*cos(2*pi*100*t), kompanderi ei ole käytössä, 8-bittinen muunnos'
%signaali = '10e-3*cos(2*pi*100*t)+0.1*cos(2*pi*15000*t), kompanderi on käytössä, 8-bittinen muunnos'
%signaali = '10e-3*cos(2*pi*100*t)+0.1*cos(2*pi*15000*t), kompanderi ei ole käytössä, 8-bittinen muunnos'
%signaali = 'Ääninäyte'

title({['Kurssi: S-72.060 K-2005 Laboratoriotyö 4, osa A: PCM, sivu 1'];...
      ['Signaali: ' signaali];...    
      ['Tekijät: ' Teekkari1 ' ja ' Teekkari2];...
      [int2str(h(3)) '.' int2str(h(2)) '.' int2str(h(1)) ' klo '...
          int2str(h(4)) ':' int2str(h(5)) ':' int2str(h(6))]});

subplot(3,2,2); f = 0:fo:Fs-fo; %Taajuusvektori.
plot(f,INf); axis([0 Fs/2 0 1.2*max(INf)]);
xlabel('Taajuus/Hz'); ylabel('Amplitudi/V');
title('Tulosignaalin X(f)'); grid on; zoom on;

subplot(3,2,3); plot(t1,h1t);
axis([min(t1) max(t1) 1.2*min(h1t) 1.2*max(h1t)]);
xlabel('Aika/s'); ylabel('Amplitudi/V');
title('Impulssivaste h(t)'); grid on; zoom on;

subplot(3,2,4); f1 = 0:fo1:Fs-fo1;
plot(f1,H1f); axis([0 Fs/2 0 1.2*max(H1f)]);
xlabel('Taajuus/Hz'); ylabel('Amplitudi/V');
title('H1(f)'); grid on; zoom on;

subplot(3,2,5); plot(t,x1); axis([0 max(t) 1.2*min(x1) 1.2*max(x1)])
xlabel('Aika/s'); ylabel('Amplitudi/V');
title('Suodatettu signaali x1(t)'); grid on; zoom on;

subplot(3,2,6); plot(f,X1f);
axis([0 Fs/2 0 1.2*max(X1f)]);
xlabel('Taajuus/Hz'); ylabel('Amplitudi/V')
title('Suodatetun signaalin X1(f)'); grid on; zoom on;


figure %2. arkki

subplot(2,3,1); plot(t,x1);
axis([0 max(t) 1.2*min(x1) 1.2*max(x1)]);
xlabel('Aika/s'); ylabel('Amplitudi/V');
title('Suodatettu tulosignaali x1(t)'); grid on; zoom on;

title({['Kurssi: S-72.060 K-2005 Laboratoriotyö 4, osa A: PCM, sivu 2'];...
      ['Signaali: ' signaali];...    
      ['Tekijät: ' Teekkari1 ' ja ' Teekkari2];...
      [int2str(h(3)) '.' int2str(h(2)) '.' int2str(h(1)) ' klo '...
          int2str(h(4)) ':' int2str(h(5)) ':' int2str(h(6))]});

xa = linspace(-1,1,length(x1))'; %suora välillä -1...+1

if kompressointi == 1;

alaki = zeros(size(x1)); %alustetaan a-lain vektori

for i = 1:length(xa); %Kelataan koko signaali arvoittain alusta loppuun.
   if abs(xa(i)) < 1/A; %Pienten arvojen laskenta.
      alaki(i,1) = xa(i)*A/(1+log(A));
   elseif xa(i) > 0 %Positiivisten suurten arvojen laskenta.
      alaki(i) = (1+log(xa(i,1)*A))/(1+log(A));
   else %Negatiivisten suurten arvojen laskenta.
      alaki(i) = -1*(1+log(abs(xa(i,1))*A))/(1+log(A));
   end
end
else
    alaki = xa;
end;

subplot(2,3,2); plot(xa,alaki);
axis([1.2*min(xa) 1.2*max(xa) 1.2*min(alaki) 1.2*max(alaki)]);
xlabel('Kompression tulo/V'); ylabel('Kompression lähtö/V')
title('A-lain kompression ominaiskäyrä'); grid on; zoom on;

subplot(2,3,3); plot(t,comp);
axis([0 max(t) 1.2*min(comp) 1.2*max(comp)]);
xlabel('Aika/s'); ylabel('Amplitudi/V')
title('Kompressoitu signaali comp(t)'); grid on; zoom on;

subplot(2,3,4); plot(t,comp);
axis([0 max(t) 1.2*min(comp) 1.2*max(comp)]);
xlabel('Aika/s'); ylabel('Amplitudi/V');
title('Kompressoitu signaali comp(t)'); grid on; zoom on;

subplot(2,3,5); plot(t,out1);
axis([0 max(t) 1.2*min(out1) 1.2*max(out1)]);
xlabel('Aika/s'); ylabel('Amplitudi/V')
title('AD-DA-muunnettu signaali out1(t)'); grid on; zoom on;

subplot(2,3,6); plot(t,comp-out1);
axis([0 max(t) 1.2*min(comp-out1) 1.2*max(comp-out1)]);
xlabel('Aika/Hz'); ylabel('Amplitudi/V')
title('Kvantisointivirhe 1 (comp(t)-out1(t))'); grid on; zoom on;

figure %3. arkki

subplot(2,3,1); plot(t,out1);
axis([0 max(t) 1.2*min(x1) 1.2*max(x1)]);
xlabel('Aika/s'); ylabel('Amplitudi/V')
title('Dekompressorin tulosignaali out1(t)'); grid on; zoom on;

title({['Kurssi: S-72.060 K-2005 Laboratoriotyö 4, osa A: PCM, sivu 3'];...
      ['Signaali: ' signaali];...    
      ['Tekijät: ' Teekkari1 ' ja ' Teekkari2];...
      [int2str(h(3)) '.' int2str(h(2)) '.' int2str(h(1)) ' klo '...
          int2str(h(4)) ':' int2str(h(5)) ':' int2str(h(6))]});

subplot(2,3,2); plot(alaki,xa);
axis([1.2*min(alaki) 1.2*max(alaki) 1.2*min(xa) 1.2*max(xa)]);
xlabel('Dekompression tulo/V'); ylabel('Dekompression lähtö/V');
title('A-lain dekompression ominaiskäyrä'); grid on; zoom on;

subplot(2,3,3); plot(t,decomp);
axis([0 max(t) 1.2*min(decomp) 1.2*max(decomp)]);
xlabel('Aika/s'); ylabel('Amplitudi/V');
title('Dekompressoitu signaali decomp(t)'); grid on; zoom on;

subplot(2,3,4); plot(t3,out3);
axis([0 max(t3) 1.2*min(out3) 1.2*max(out3)]);
xlabel('Aika/s'); ylabel('Amplitudi/V');
title('Lähtösignaali out(t)'); grid on; zoom on;

subplot(2,3,5); plot(f3,OUT3f);
axis([0 Fs/2 0 1.2*max(OUT3f)]);
xlabel('Taajuus/Hz'); ylabel('Amplitudi/V');
title('Lähtösignaalin spektri OUT(f)'); grid on; zoom on;

subplot(2,3,6); plot(t3,d);
axis([0 max(t3) 1.2*min(d) 1.2*max(d)]);
xlabel('Taajuus/Hz'); ylabel('Amplitudi/V')
title('Lopullinen kvantisointivirhe out(t) - x1(t)'); grid on; zoom on;

soundsc(out3,Fs); %Halutun signaalin kuuntelu.

