% camp96 = segnale con banda di 20 kHz campionato a 96 kHz (in totale 2^17=131072 campioni)
%
% si campiona a frequenza ridotta, si interpola a 96 kHz e si confronta con l'originale
%
% si puo' prefiltrare, nei casi di campionamento a frequenza molto bassa (se prefiltro = 1)
%
% si possono salvare i file nel formato .wav (se salvataggio = 1)
%
% COMMENTO: tipico esempio di programma scritto di fretta per un solo caso e poi riutilizzato (male) per N casi!!
%           si dovrebbero usare due o tre semplici funzioni per le operazioni fondamentali (e ripetitive);
%           sarebbe bene introdurre anche variabili ausiliarie, per semplificare la lettura
%
disp(' '),disp(' DOPO OGNI FIGURA PREMERE RETURN PER PROSEGUIRE'),drawnow
close all
if exist('prefiltro')~=1,prefiltro=0;end
if exist('salvataggio')~=1,salvataggio=0;end
disp(' ')
if prefiltro==1,disp(' prefiltro'),disp(' '),end
if salvataggio==1,disp(' salvataggio'),disp(' '),end
[x fc nbit]=wavread('camp96');x=x';
N=length(x);
dt=1/fc;               % passo di campionamento
E_segnale=x*x'*dt;     % energia del segnale
plot(dt*(0:N-1),x),xlabel('t')
title(['segnale originale      energia = ' num2str(E_segnale)]),figure(1),pause
T=N*dt;                % durata complessiva
df=1/T;                % passo in frequenza
f=df*(0:N-1);
trasf=abs(fft(x));semilogy(f,trasf)
xlabel('f'),title('modulo trasformata del segnale campionato a 96 kHz'),figure(1),pause

if prefiltro==0
    
% filtro passo basso ideale con banda 24 kHz
% troncato e ritardato di Ntr campioni (in figura sono mostrati Nview campioni)
Ntr=300;Nview=50;
h=sin(.5*pi*(-Ntr:Ntr)+1e-6)./(.5*pi*(-Ntr:Ntr)+1e-6);
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'.'),hold on
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'r'),grid
xlabel('t'),title('risposta impulsiva (troncata) passo basso ideale 24 kHz'),figure(1),pause
clf
y=zeros(1,N);y(1:2:N)=x(1:2:N);
trasf=abs(fft(2*y));semilogy(f(1:N/2),trasf(1:N/2))
xlabel('f'),title('modulo trasformata del segnale campionato a 48 kHz'),figure(1),pause
z=conv(y,h);
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,wavwrite(z,fc,'camp48_ideale'),end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause
Ntr=100;Nview=50;
H=zeros(1,N);
n20=round(20000/df);
H(1:n20)=1;
H(n20+1:N/4+1)=cos(.5*pi*(1:(N/4-n20)+1)/(N/4-n20)).^2;
H(N:-1:3*N/4+1)=H(2:N/4+1);
plot(f(1:N/2),H(1:N/2)),xlabel('f'),title('passa basso con taglio a 20-24 kHz'),figure(1),pause
h=2*real(ifft(H));
h=[h(N+1-Ntr:N) h(1:Ntr+1)];
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'.'),hold on
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'r'),grid
xlabel('t'),title('risposta impulsiva passo basso con taglio a 20-24 kHz'),figure(1),pause
clf
y=zeros(1,N);y(1:2:N)=x(1:2:N);
z=conv(y,h);
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,wavwrite(z,fc,'camp48'),end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause
Ntr=100;Nview=100;
H=zeros(1,N);
n10=round(10000/df);
H(1:n10)=1;
H(n10+1:N/8+1)=cos(.5*pi*(1:(N/8-n10)+1)/(N/8-n10)).^2;
H(N:-1:7*N/8+1)=H(2:N/8+1);
h=4*real(ifft(H));
h=[h(N+1-Ntr:N) h(1:Ntr+1)];
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'.'),hold on
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'r'),grid
xlabel('t'),title('risposta impulsiva passo basso con taglio a 10-12 kHz'),figure(1),pause
clf
y=zeros(1,N);y(1:4:N)=x(1:4:N);
trasf=abs(fft(4*y));semilogy(f(1:N/4),trasf(1:N/4))
xlabel('f'),title('modulo trasformata del segnale campionato a 24 kHz'),figure(1),pause
z=conv(y,h);
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,wavwrite(z,fc,'camp24'),end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause

end

Ntr=200;Nview=200;
H=zeros(1,N);
n5=round(5000/df);
H(1:n5)=1;
H(n5+1:N/16+1)=cos(.5*pi*(1:(N/16-n5)+1)/(N/16-n5)).^2;
H(N:-1:15*N/16+1)=H(2:N/16+1);
h=8*real(ifft(H));
h=[h(N+1-Ntr:N) h(1:Ntr+1)];
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'.'),hold on
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'r'),grid
xlabel('t'),title('risposta impulsiva passo basso con taglio a 5-6 kHz'),figure(1),pause
clf
if prefiltro==1
    xx=.125*conv(x,h);
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:8:N)=xx(1:8:N);
else
    y=zeros(1,N);y(1:8:N)=x(1:8:N);
end
trasf=abs(fft(8*y));semilogy(f(1:N/8),trasf(1:N/8))
if prefiltro==0
    xlabel('f'),title('modulo trasformata del segnale campionato a 12 kHz'),figure(1),pause
else
    xlabel('f'),title('modulo trasformata del segnale prefiltrato e campionato a 12 kHz'),figure(1),pause
end
z=conv(y,h);
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1 & prefiltro==0,wavwrite(z,fc,'camp12'),end
if salvataggio==1 & prefiltro==1,wavwrite(z,fc,'camp12_prefiltrato'),end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause
Ntr=400;Nview=400;
H=zeros(1,N);
n2_5=round(2500/df);
H(1:n2_5)=1;
H(n2_5+1:N/32+1)=cos(.5*pi*(1:(N/32-n2_5)+1)/(N/32-n2_5)).^2;
H(N:-1:31*N/32+1)=H(2:N/32+1);
h=16*real(ifft(H));
h=[h(N+1-Ntr:N) h(1:Ntr+1)];
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'.'),hold on
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'r'),grid
xlabel('t'),title('risposta impulsiva passo basso con taglio a 2.5-3 kHz'),figure(1),pause
clf
if prefiltro==1
    xx=.0625*conv(x,h);
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:16:N)=xx(1:16:N);
else
    y=zeros(1,N);y(1:16:N)=x(1:16:N);
end
trasf=abs(fft(16*y));semilogy(f(1:N/16),trasf(1:N/16))
if prefiltro==0
    xlabel('f'),title('modulo trasformata del segnale campionato a 6 kHz'),figure(1),pause
else
    xlabel('f'),title('modulo trasformata del segnale prefiltrato e campionato a 6 kHz'),figure(1),pause
end
z=conv(y,h);
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1 & prefiltro==0,wavwrite(z,fc,'camp6'),end
if salvataggio==1 & prefiltro==1,wavwrite(z,fc,'camp6_prefiltrato'),end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause
Ntr=800;Nview=800;
H=zeros(1,N);
n1_25=round(1250/df);
H(1:n1_25)=1;
H(n1_25+1:N/64+1)=cos(.5*pi*(1:(N/64-n1_25)+1)/(N/64-n1_25)).^2;
H(N:-1:63*N/64+1)=H(2:N/64+1);
h=32*real(ifft(H));
h=[h(N+1-Ntr:N) h(1:Ntr+1)];
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'.'),hold on
plot(dt*(-Nview:Nview),h(Ntr+1-Nview:Ntr+1+Nview),'r'),grid
xlabel('t'),title('risposta impulsiva passo basso con taglio a 1.25-1.5 kHz'),figure(1),pause
clf
if prefiltro==1
    xx=.03125*conv(x,h);
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:32:N)=xx(1:32:N);
else
    y=zeros(1,N);y(1:32:N)=x(1:32:N);
end
trasf=abs(fft(32*y));semilogy(f(1:N/32),trasf(1:N/32))
if prefiltro==0
    xlabel('f'),title('modulo trasformata del segnale campionato a 3 kHz'),figure(1),pause
else
    xlabel('f'),title('modulo trasformata del segnale prefiltrato e campionato a 3 kHz'),figure(1),pause
end
z=conv(y,h);
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1 & prefiltro==0,wavwrite(z,fc,'camp3'),end
if salvataggio==1 & prefiltro==1,wavwrite(z,fc,'camp3_prefiltrato'),end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1)