// Questo script caratterizza un sistema lineare a partire dallo stimolo impartitogli
// e dalla risposta del sistema stesso.
// Stimolo e risposta devono essere contenuti in un file wav stereo.
// Il canale sinistro deve contenere lo stimolo,
// il destro la corrispondente risposta.
// Il nome del file wav viene richiesto interattivamente alla console.
// Inoltre compensa per le carattristiche del microfono usato per la presa
// Che si suppone descritto in un file

// Definizione di una funzione.
// dato un valore e un riferimento, restituisce il valore il dB
function y=db(ampl,ref)
	y = 20*log(ampl/ref);
endfunction

// Procedura main:
// Cancella eventuali precedenti grafici
xdel(winsid());
x = []; // definisce un vettore vuoto, per stare nel loop
while length(x)==0
	filein=input("Nome del file wav contenente i dati: ","string"); // chiede il nome del file contenente i dati
	[x,ierr]=fileinfo(filein); // dati sul file appena definito. x contiene la lunghezza del file.
	// se la lunghezza è zero, il file non esiste o è vuoto. Torna a chiedere.
	if length(x)==0 then
		print(%io(2),"Il file non esiste. Abort."); // n modo per mandare messaggi alla console
	end
end
framex=input("Esponente in base 2 della lunghezza della frame (es. 16 per indicare 65536): "); // vogliamo frame lunghe potenze di due
stacksize(100000000); // aumenta lo memoria a disposizione dello scilab.
data = wavread(filein); // leggo il file contenente stimolo e risposta. data è una matrice a due colonne.
in = data(:,2); // lo stimolo (in) nel canale destro, 2 colonna della matrice
framex = floor(framex); // caso mai si fosse dato un esponente decimale, il valore viene troncato.
if (framex < 4) | (framex > 19) then
	print(%io(2),"Esponente fuori dai limiti (4 <= exp <= 18). Abort."); // frame troppo corta o troppo lunga
	return
end
	
framel = 2^framex; // l'effettiva lunghezza della frame dall'esponente.
delete('all');
scf(0); clf(0); // crea figura 0, e la cancella (se era già esistente) 
plot2d(data(1:framel)); // plot di stimolo e risposta
xtitle(filein + " - Prima frame di ingresso e di uscita"); // titolo delal figura 0
xgrid(15); // griglia di colore verde (15)

// Calcolo della funzione di trasferimento mediando diverse frame

nlun = length(data(:,1)); // lunghezza del file (può contenere più di una frame)
nframes = floor(nlun/framel); // quante frame vi sono contenute
in = zeros(framel,1); // ingresso creato come vettore di zeri
out = in; // idem per l'uscita
// media i segnali, sommando le varie frames
for i=1:nframes-1
	start = (nframes - 1) * framel + 1;
	endt   = (nframes * framel);
	in = in   +  data(start:endt,2);
	out = out + data(start:endt,1);
end
// frequenze da 10 Hz a 14000
n0 = floor(10/22050*(framel/2)); // punto corrispondente a 10 Hz
n1 = floor(14000/22050*(framel/2)); //punto corrispondente a 14000 Hz
frq = linspace(10,14000,floor(n1-n0)+1)';
scf(1); clf(1); // figura 1, e pulisci
subplot(2,1,1); // definisce una sezione dell'area di plotaggio
plot2d(in); // plot dell'ingresso nella sezione
xgrid(15); // griglia verde
xtitle(filein + " - Ingresso mediato:"); // Titolo della sezione
subplot(2,1,2); // altra sezione, sotto la precedente
plot2d(out); // Qui plot dell'uscita
xgrid(15); // griglia verde
xtitle(filein + " - Uscita mediata:"); // Titolo della sezione
scf(2); clf(2); // figura 2
ins=fft(in,-1); // Trasformata dell'ingresso (complessa)
outs=fft(out,-1); // idem per l'uscita
// ouF = [db(abs(ins(n0:framel/2)),1) db(abs(outs(n0:framel/2)),1)];
// plot2d(frq,ouF,logflag="ln",leg="Ingresso@Uscita");
xmin = 10; // frequenza minima
ymin = -120; // Ampiezza minima (dB)
xmax = 100000; // frequenza massima
ymax = 0; // massima ampiezza (dB)
outdb = db(abs(outs(n0:n1)),1); // Ampiezza = valore assoluto della trasformata. Espressa in dB
outdb = outdb - max(outdb); // Riferita allo 0dB
plot2d(frq,outdb,logflag="ln",rect=[xmin,ymin,xmax,ymax]); // plot. La frequenza è in scala logaritmica
xtitle(filein + " - Spettro uscita");
xgrid(15); // griglia verde
scf(3); clf(3); // Finestra 3
ins = ins / max(abs(ins)); // normalizzo l'ingresso a 1.
// plot2d(frq,abs(outs(n0:framel/2)./ins(n0:framel/2)),logflag="ll");
// Funzione di trasferimento: rapporto (compelsso) tra trasformata dell'uscita e quella dell'ingresso 
tf = db(abs(outs(n0:n1)./ins(n0:n1)),1); // Dato che c'è una costante arbitraria, solo la forma è valida
// compensazione microfono. Qui il file si riferisce ad uno dei microfoni a disposizione dell'aula.
mic=fscanfMat("akg08638.prn"); // si suppone che questo file contenga i dati di calibrazione del microfono
maxmic = max(mic);
mic = maxmic - mic; // normalizziamo il microfono
tf = tf - mic ; // sottraiamo il microfono
tf = tf - max(tf); // normalizziamo la funzione di trasferimento corretta
plot2d(frq,tf(1:length(tf)),logflag="ln",rect=[xmin,ymin,xmax,ymax]); // Plot della funzione di trasf.
xtitle(filein + " - Funzione di trasferimento"); // Titolo del grafico
xgrid(15); // Griglia verde
// Passiamo al dominio del tempo:
// Calcolo della risposta all'impulso
ft = outs./ins; // Torniamo alla funzionedi trasferimento complessa 
pulser = fft(ft,1); // antitrasformata: otteniamo la risposta all'impulsa.
pulser = real(pulser); // solo la parte reale
scf(4); clf(4); // figura 4
lunp =floor(framel*.9) // Tagliamo la coda che spesso è spuria
pulser = pulser(1:lunp); // idem
fprintfMat(filein+'.txt',pulser); // Scrivo un file di testo con i dati della risposta
xmp = 0;
ymp=min(pulser); // Minimo della risposta
xMp=(lunp*1/44100); // delta dt tempo, per riportarsi a 44100
yMp=max(pulser); // Massimo della risposta
time= linspace(0,xMp,lunp);
plot2d(time,pulser,logflag='nn',rect=[xmp,ymp,xMp,yMp]);// Plot della risposta all'impulso
xgrid(15); // Griglia verde