clear all
close all
plat = 'octave'; %put either 'octave' or 'matlab'

%function y = quantize(x)
%%y = round(x); % arrotondamento
%% y = floor(x); %troncamento
%y = floor(x + (x < 0)); %troncamento simmetrico
%end;

if (plat == 'matlab')
  [s,Fs]=wavread('clarinetto.wav');
else
  [s,Fs]=auload('clarinetto.wav');
end;

%16-bit quantization
nbits=16;
q=2/(2^nbits-1)

s_16=q*quantize(s./q);

viz_frame=[10000:11000];
figure(2)
subplot(3,1,1)
plot(viz_frame,s(viz_frame))
title('Original sample')

nbits=7;
q=2/(2^nbits-1)
s_q=q*quantize(s./q);
hold on
subplot(3,1,2)
plot(viz_frame,s_q(viz_frame),'b')
title('Quantized sample (7 bits)')

%dither
add_dith = -q/2+q*rand(length(s),1);
q=2/(2^nbits-1)
s_qdith=q*quantize((s+add_dith)./q);
subplot(3,1,3)
plot(viz_frame,s_qdith(viz_frame),'b')
title('Additive dither')
hold off;

pause;

figure(3)
N=2*1024;
axis([0, N/2, -20, 20]);
subplot(3,1,1)
S=fft(s(viz_frame),N);
plot(20*log10(abs(S(1:N/2+1))))
title('Original sample')
subplot(3,1,2)
S=fft(s_q(viz_frame),N);
plot(20*log10(abs(S(1:N/2+1))))
title('Quantized sample (7 bits)')
subplot(3,1,3)
S=fft(s_qdith(viz_frame),N);
plot(20*log10(abs(S(1:N/2+1))))
title('Additive dither')


%shaping FIR filter
%ref: Stanley P. Lipshitz, John Vanderkooy, and Robert A. Wannamaker, 
%"Minimally Audible Noise Shaping",  J. Audio Eng. Soc., 
%Vol. 39, No. 11, November 1991.
B=[2.033 -2.165 1.959 -1.590 0.6149];
A=1;

%spectrum shaping
eprec=0; 
eprecf=0;
Z=zeros(length(B)-1,1);
s_qreshape=zeros(1,length(s));
s_qreshapef=zeros(1,length(s));
for (i=1:length(s))
   s_qreshape(i)=q*quantize((s(i)+eprec)./q);
   s_qreshapef(i)=q*quantize((s(i)+eprecf)./q);
   e=s(i)+eprec-s_qreshape(i);
   eprec=e;   % noise shaping without filter
   e=s(i)+eprecf-s_qreshapef(i);
   [eh,Z]=filter(B,A,e,Z); % noise shaping with filter
   eprecf=eh;
end


figure(4)
N=2*1024;
subplot(3,1,1)
S=fft(s(viz_frame),N);
plot(20*log10(abs(S(1:N/2+1))))
title('Original sample')
subplot(3,1,2)
S=fft(s_qreshape(viz_frame),N);
plot(20*log10(abs(S(1:N/2+1))))
title('Noise shaping dither')
subplot(3,1,3)
S=fft(s_qreshapef(viz_frame),N);
plot(20*log10(abs(S(1:N/2+1))))
title('Noise shaping dither with filter')

if (plat == 'matlab') 
  wavwrite(s_q,Fs,16,'s_q.wav');
  wavwrite(s_qdith,Fs,16,'s_qdith.wav');
  wavwrite(s_qreshape,Fs,16,'s_qreshape.wav');
  wavwrite(s_qreshapef,Fs,16,'s_qreshapef.wav');
else
  ausave('s_q.wav',s_q,Fs,'short');
  ausave('s_qdith.wav',s_qdith,Fs,'short');
  ausave('s_qreshape.wav',s_qreshape,Fs,'short');
  ausave('s_qreshapef.wav',s_qreshapef,Fs,'short');
end

%lowpass filter
%[Blp,Alp] = butter(25,0.6,'low');
%freqz(Blp,Alp)
%s_qreshape_lp=filter(Blp,Alp,s_qreshape);