% Si seleziona un file wav di partenza, e lo si sottocampiona a frequenze ridotte via via di un
% fattore 2 fino ad un minimo di 2 kHz, si interpola alla frequenza fc e si confronta con l'originale
%
% si puo' prefiltrare (se prefiltro = 1) e si possono salvare i file nel formato .wav (se salvataggio = 1)
% default 0 per entrambi

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
[file,path]=uigetfile('*.wav','Scegli il file wav:');
file=file(1:(length(file)-4)); % si toglie .wav dal nome del file
[x,fc,nbit]=wavread([file '.wav']);

N=length(x);
x=x(1:2^(floor(log2(N))))'; % se ne fa una potenza di 2
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);
% if rem(N,2),
%    f=(N-1)/2+df*(0:N-1);
% else
%    f=N/2+df*(0:N-1);
% end   
trasf=abs(fft(x)); %trasf=fftshift(trasf);
semilogy(f,trasf)
xlabel('f'),title(['modulo trasformata del segnale campionato a ' num2str(fc) ' kHz']),figure(1),pause

Ntr=50;Nview=50;
alpha=0.2;         % roll-off

% sottocampionamento a fc/2
Ntr=2*Ntr;Nview=2*Nview;
fattore=2;
if fc/fattore<2000, break; end
fT=fc/2/fattore;    % Frequenza di taglio
fIni=fT/(1+alpha);  % freq. di inizio transizione
Initrans=max(find(f<=fIni));
Finetrans=max(find(f<=fT));
H=ones(1,Initrans-1);
H=[H cos(.5*pi*(0:(Finetrans-Initrans))/(Finetrans-Initrans)).^2];
H=[H zeros(1,N/2-length(H))];
H=[H 0 fliplr(H)];
H=H(1:N);
h=fattore*real(ifft(H)); % N*df*fattore/fc=fattore
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 ' num2str(fT/1000) ' kHz']),figure(1),pause
if prefiltro==1
    xx=conv(x,h)*dt*fc/fattore;
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:fattore:N)=xx(1:fattore:N);
else
    y=zeros(1,N);y(1:fattore:N)=x(1:fattore:N);
end
clf
trasf=abs(fft(fattore*y));semilogy(f(1:N/fattore),trasf(1:N/fattore))
xlabel('f'),title(['modulo trasformata del segnale campionato a ' num2str(fc/fattore/1000) ' kHz']),figure(1),pause
z=conv(y,h); % segnale ricostruito a fc da quello sottocampionato
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,
   if prefiltro==1,
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz_pref.wav']);
   else
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz.wav']);
   end
end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause


% sottocampionamento a fc/4
Ntr=2*Ntr;Nview=2*Nview;
fattore=4;
if fc/fattore<2000, break; end
fT=fc/2/fattore;    % Frequenza di taglio
fIni=fT/(1+alpha);  % freq. di inizio transizione
Initrans=max(find(f<=fIni));
Finetrans=max(find(f<=fT));
H=ones(1,Initrans-1);
H=[H cos(.5*pi*(0:(Finetrans-Initrans))/(Finetrans-Initrans)).^2];
H=[H zeros(1,N/2-length(H))];
H=[H 0 fliplr(H)];
H=H(1:N);
h=fattore*real(ifft(H)); % N*df*fattore/fc=fattore
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 ' num2str(fT/1000) ' kHz']),figure(1),pause
if prefiltro==1
    xx=conv(x,h)*dt*fc/fattore;
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:fattore:N)=xx(1:fattore:N);
else
    y=zeros(1,N);y(1:fattore:N)=x(1:fattore:N);
end
clf
trasf=abs(fft(fattore*y));semilogy(f(1:N/fattore),trasf(1:N/fattore))
xlabel('f'),title(['modulo trasformata del segnale campionato a ' num2str(fc/fattore/1000) ' kHz']),figure(1),pause
z=conv(y,h); % segnale ricostruito a fc da quello sottocampionato
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,
   if prefiltro==1,
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz_pref.wav']);
   else
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz.wav']);
   end
end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause


% sottocampionamento a fc/8
Ntr=2*Ntr;Nview=2*Nview;
fattore=8;
if fc/fattore<2000, break; end
fT=fc/2/fattore;    % Frequenza di taglio
fIni=fT/(1+alpha);  % freq. di inizio transizione
Initrans=max(find(f<=fIni));
Finetrans=max(find(f<=fT));
H=ones(1,Initrans-1);
H=[H cos(.5*pi*(0:(Finetrans-Initrans))/(Finetrans-Initrans)).^2];
H=[H zeros(1,N/2-length(H))];
H=[H 0 fliplr(H)];
H=H(1:N);
h=fattore*real(ifft(H)); % N*df*fattore/fc=fattore
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 ' num2str(fT/1000) ' kHz']),figure(1),pause
if prefiltro==1
    xx=conv(x,h)*dt*fc/fattore;
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:fattore:N)=xx(1:fattore:N);
else
    y=zeros(1,N);y(1:fattore:N)=x(1:fattore:N);
end
clf
trasf=abs(fft(fattore*y));semilogy(f(1:N/fattore),trasf(1:N/fattore))
xlabel('f'),title(['modulo trasformata del segnale campionato a ' num2str(fc/fattore/1000) ' kHz']),figure(1),pause
z=conv(y,h); % segnale ricostruito a fc da quello sottocampionato
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,
   if prefiltro==1,
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz_pref.wav']);
   else
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz.wav']);
   end
end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause

% sottocampionamento a fc/16
Ntr=2*Ntr;Nview=2*Nview;
fattore=16;
if fc/fattore<2000, break; end
fT=fc/2/fattore;    % Frequenza di taglio
fIni=fT/(1+alpha);  % freq. di inizio transizione
Initrans=max(find(f<=fIni));
Finetrans=max(find(f<=fT));
H=ones(1,Initrans-1);
H=[H cos(.5*pi*(0:(Finetrans-Initrans))/(Finetrans-Initrans)).^2];
H=[H zeros(1,N/2-length(H))];
H=[H 0 fliplr(H)];
H=H(1:N);
h=fattore*real(ifft(H)); % N*df*fattore/fc=fattore
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 ' num2str(fT/1000) ' kHz']),figure(1),pause
if prefiltro==1
    xx=conv(x,h)*dt*fc/fattore;
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:fattore:N)=xx(1:fattore:N);
else
    y=zeros(1,N);y(1:fattore:N)=x(1:fattore:N);
end
clf
trasf=abs(fft(fattore*y));semilogy(f(1:N/fattore),trasf(1:N/fattore))
xlabel('f'),title(['modulo trasformata del segnale campionato a ' num2str(fc/fattore/1000) ' kHz']),figure(1),pause
z=conv(y,h); % segnale ricostruito a fc da quello sottocampionato
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,
   if prefiltro==1,
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz_pref.wav']);
   else
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz.wav']);
   end
end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause

% sottocampionamento a fc/32
Ntr=2*Ntr;Nview=2*Nview;
fattore=32;
if fc/fattore<2000, break; end
fT=fc/2/fattore;    % Frequenza di taglio
fIni=fT/(1+alpha);  % freq. di inizio transizione
Initrans=max(find(f<=fIni));
Finetrans=max(find(f<=fT));
H=ones(1,Initrans-1);
H=[H cos(.5*pi*(0:(Finetrans-Initrans))/(Finetrans-Initrans)).^2];
H=[H zeros(1,N/2-length(H))];
H=[H 0 fliplr(H)];
H=H(1:N);
h=fattore*real(ifft(H)); % N*df*fattore/fc=fattore
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 ' num2str(fT/1000) ' kHz']),figure(1),pause
if prefiltro==1
    xx=conv(x,h)*dt*fc/fattore;
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:fattore:N)=xx(1:fattore:N);
else
    y=zeros(1,N);y(1:fattore:N)=x(1:fattore:N);
end
clf
trasf=abs(fft(fattore*y));semilogy(f(1:N/fattore),trasf(1:N/fattore))
xlabel('f'),title(['modulo trasformata del segnale campionato a ' num2str(fc/fattore/1000) ' kHz']),figure(1),pause
z=conv(y,h); % segnale ricostruito a fc da quello sottocampionato
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,
   if prefiltro==1,
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz_pref.wav']);
   else
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz.wav']);
   end
end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause


% sottocampionamento a fc/64
Ntr=2*Ntr;Nview=2*Nview;
fattore=64;
if fc/fattore<2000, break; end
fT=fc/2/fattore;    % Frequenza di taglio
fIni=fT/(1+alpha);  % freq. di inizio transizione
Initrans=max(find(f<=fIni));
Finetrans=max(find(f<=fT));
H=ones(1,Initrans-1);
H=[H cos(.5*pi*(0:(Finetrans-Initrans))/(Finetrans-Initrans)).^2];
H=[H zeros(1,N/2-length(H))];
H=[H 0 fliplr(H)];
H=H(1:N);
h=fattore*real(ifft(H)); % N*df*fattore/fc=fattore
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 ' num2str(fT/1000) ' kHz']),figure(1),pause
if prefiltro==1
    xx=conv(x,h)*dt*fc/fattore;
    xx=xx(Ntr+1:N+Ntr);
    y=zeros(1,N);y(1:fattore:N)=xx(1:fattore:N);
else
    y=zeros(1,N);y(1:fattore:N)=x(1:fattore:N);
end
clf
trasf=abs(fft(fattore*y));semilogy(f(1:N/fattore),trasf(1:N/fattore))
xlabel('f'),title(['modulo trasformata del segnale campionato a ' num2str(fc/fattore/1000) ' kHz']),figure(1),pause
z=conv(y,h); % segnale ricostruito a fc da quello sottocampionato
z=z(Ntr+1:N+Ntr);
E_errore=(x-z)*(x-z)'*dt;
if salvataggio==1,
   if prefiltro==1,
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz_pref.wav']);
   else
      wavwrite(z,fc,nbit,[file num2str(round(fc/fattore/1000)) 'kHz.wav']);
   end
end
plot(dt*(0:N-1),x-z),xlabel('t'),title(['errore      energia = ' num2str(E_errore)]),figure(1),pause