%Demo program for illustrating down and up sampling using different 
%prefilter/reconstruction filters.
%compare the spectrum of the original and down-/up-sampled signals.

%read in a sound file
echo on;
[y,fs]=wavread('chimes');
pause;



%Example 1: decimate by a factor of 2 and then interpolate back by a factor of 2
%both using good filters

%plot waveform and spectrum

figure(1);
sound(y,fs);
pause;
subplot(3,4,1),plot(y);axis tight;title('y');
subplot(3,4,2),psd(y,256,fs);axis([0,fs/2,-60,0]);ylabel('');xlabel('');title('psd-y');
pause;

%decimate by factor of 2, using default FIR filter (order=30)
x21=decimate(y,2,'fir');
pause;

sound(x21,fs/2);
pause;

%plot waveform and spectrum of downsampled signal
subplot(3,4,5),plot(x21);axis tight;title('x21');
subplot(3,4,6),psd(x21,256,fs/2);axis([0,fs/2,-60,0]);ylabel('');xlabel('');title('psd-x21');
pause;

%find the actual  filter used, plot is temporal response and spectrum
%the filter should have cut-off at 0.5 
%note in matlab, 1 corresponds to fs/2 (maximum digital frequency)
%therefore, 0.5 corresponds to fs/4 (half band filter)
h21=fir1(30,0.5);
h21
subplot(3,4,7),bar(-15:1:15,h21,0.02);title('h21');axis([-15,15,-0.1,0.6]);
subplot(3,4,8),[h,f]=freqz(h21,1);
plot(f/pi,20*log10(abs(h)));axis tight;title('freq-h21');ylabel('');xlabel('');

pause;

%interpolate x21 (decimated by a factor of 2 using length 30 filter) 
% by factor of 2, using length 4 filter, with cut-off at 0.5 (half-band)

[z21,b21]=interp(x21,2,4,0.5);
pause;

sound(z21,fs);pause;

%plot waveform and spectrum of interpolated signal
subplot(3,4,9),plot(z21);axis tight;title('z21');
subplot(3,4,10),psd(z21,256,fs);axis tight;ylabel('');xlabel('');title('psd-z21');
pause;

%display b21, plot its temporal response and spectrum
%note that length N in interp() function means that the filter has N non-zero values on each side
%of the origin, the origin  is always 1.
%The actual filter length is 4*N+1, with zeros between non-zero ones.
b21
pause;
subplot(3,4,11),bar(-8:1:8,b21,0.02);title('b21');axis([-8,8,-0.2,1.2]);
subplot(3,4,12),[h,f]=freqz(b21,1);
plot(f/pi,20*log10(abs(h)));axis tight;title('freq-b21');ylabel('');xlabel('');
pause;

%finally compare the three signals in sound quality
sound(y,fs); pause;
sound(x21,fs/2); pause;
sound(z21,fs); pause;

%Example 2: decimate by a factor of 2 and then interpolate back by a factor of 2
%both using simple, not so good filters

figure(2);
sound(y,fs);
pause;
subplot(3,4,1),plot(y);axis tight;title('y');
subplot(3,4,2),psd(y,256,fs);axis([0,fs/2,-60,0]);ylabel('');xlabel('');title('psd-y');
pause;

%decimate by factor of 2, using a FIR filter with order=1
x22=decimate(y,2,1,'fir');
pause;

sound(x22,fs/2);
pause;
%plot waveform and spectrum of downsampled signal
subplot(3,4,5),plot(x22);axis tight;title('x22');
subplot(3,4,6),psd(x22,256,fs/2);axis([0,fs/2,-60,0]);ylabel('');xlabel('');title('psd-x22');


%find the actual  filter used, plot is temporal response and spectrum
h22=fir1(1,0.5);
subplot(3,4,7),bar(0:1,h22,0.02);title('h22');axis([-15,15,-0.1,0.6]);
subplot(3,4,8),[h,f]=freqz(h22,1);
plot(f/pi,20*log10(abs(h)));axis tight;title('freq-h22');ylabel('');xlabel('');

pause;

%interpolate x22 by factor of 2, using length 1 filter, with cut-off at 0.5 (half-band)
[z22,b22]=interp(x22,2,1,0.5);
pause;

sound(z22,fs);pause;

%plot waveform and spectrum of interpolated signal
subplot(3,4,9),plot(z22);axis tight;title('z22');
subplot(3,4,10),psd(z22,256,fs);axis tight;ylabel('');xlabel('');title('psd-z22');
pause;

%display b22, plot its temporal response and spectrum
%note that length N in interp() function means that the filter has N non-zero values on each side
%of the origin, the origin  is always 1.
%The actual filter length is 4*N+1, with zeros between non-zero ones.
b22
pause;
subplot(3,4,11),bar(-2:1:2,b22,0.02);title('b22');axis([-8,8,-0.2,1.2]);
subplot(3,4,12),[h,f]=freqz(b22,1);
plot(f/pi,20*log10(abs(h)));axis tight;title('freq-b22');ylabel('');xlabel('');
pause;

%finally compare the three signals in sound quality
sound(y,fs); pause;
sound(x22,fs/2); pause;
sound(z22,fs); pause;

%Example 3: decimate by a factor of 4, then interpolate by a factor of 4
% both using good default filters


figure(3);
sound(y,fs); pause;
subplot(3,4,1),plot(y);axis tight;title('y');
subplot(3,4,2),psd(y,256,fs);axis([0,fs/2,-60,0]);ylabel('');xlabel('');title('psd-y');
pause;

%decimate using default filter (order=30)
x41=decimate(y,4,'fir');
pause;

sound(x41,fs/4);
pause;
subplot(3,4,5),plot(x41);axis tight;title('x41');
subplot(3,4,6),psd(x41,256,fs/4);axis([0,fs/2,-60,0]);ylabel('');xlabel('');title('psd-x41');

%note that for decimation by factor of 4, the cut-off should be 0.25
%which corresponds to fs/2/4

h41=fir1(30,0.25);
subplot(3,4,7),bar(-15:1:15,h41,0.02);title('h41');axis([-15,15,-0.1,0.6]);
subplot(3,4,8),[h,f]=freqz(h41,1);
plot(f/pi,20*log10(abs(h)));axis tight;title('freq-h41');ylabel('');xlabel('');

pause;

%interpolate x41 by factor of 2, using length 4 filter, with cut-off at 0.25 
%note that for interpolation by factor of 4, the cut-off should be 0.25
%which corresponds to fs/2/4

[z41,b41]=interp(x41,4,4,0.25);
pause;

sound(z41,fs);pause;

%plot waveform and spectrum of interpolated signal
subplot(3,4,9),plot(z41);axis tight;title('z41');
subplot(3,4,10),psd(z41,256,fs);axis tight;ylabel('');xlabel('');title('psd-z41');
pause;

%display b41, plot its temporal response and spectrum
%note that length N in interp() function with interpolation factor R
%means that the filter has N groups of non-zero values on each side
%of the origin, the origin  is always 1.
%each group has $R-1$ non-zero values, followed by a zero
%The actual filter length is 2*R*N+1
b41
pause;
subplot(3,4,11),bar(-16:1:16,b41,0.02);title('b41');axis([-16,16,-0.2,1.2]);
subplot(3,4,12),[h,f]=freqz(b41,1);plot(f/pi,20*log10(abs(h)));axis tight;title('freq-b41');ylabel('');xlabel('');
pause;

%finally compare the three signals in sound quality
sound(y,fs); pause;
sound(x41,fs/4); pause;
sound(z41,fs); pause;




%Example 4: decimate by a factor of 4, then interpolate by a factor of 4
% both using not-so-good filters


figure(4);
sound(y,fs); pause;
subplot(3,4,1),plot(y);axis tight;title('y');
subplot(3,4,2),psd(y,256,fs);axis([0,fs/2,-60,0]);ylabel('');xlabel('');title('psd-y');
pause;

%decimate using default filter (order=3)
x42=decimate(y,4,3,'fir');
pause;

sound(x42,fs/4);
pause;
subplot(3,4,5),plot(x42);axis tight;title('x42');
subplot(3,4,6),psd(x42,256,fs/4);axis([0,fs/2,-60,0]);ylabel('');xlabel('');title('psd-x42');

%note that for decimation by factor of 4, the cut-off should be 0.25
%which corresponds to fs/2/4
h42=fir1(3,0.25);
subplot(3,4,7),bar(-1:1:2,h42,0.02);title('h42');axis([-15,15,-0.1,0.6]);
subplot(3,4,8),[h,f]=freqz(h42,1);
plot(f/pi,20*log10(abs(h)));axis tight;title('freq-h42');ylabel('');xlabel('');

pause;

%interpolate x42 by factor of 4, using length 1 filter, with cut-off at 0.25 
%note that for interpolation by factor of 4, the cut-off should be 0.25
%which corresponds to fs/2/4
[z42,b42]=interp(x42,4,1,0.25);
pause;

sound(z42,fs);pause;

%plot waveform and spectrum of interpolated signal
subplot(3,4,9),plot(z42);axis tight;title('z42');
subplot(3,4,10),psd(z42,256,fs);axis tight;ylabel('');xlabel('');title('psd-z42');
pause;

%display b41, plot its temporal response and spectrum
%note that length N in interp() function with interpolation factor R
%means that the filter has N groups of non-zero values on each side
%of the origin, the origin  is always 1.
%each group has $R-1$ non-zero values, followed by a zero
%The actual filter length is 2*R*N+1
b42
pause;
subplot(3,4,11),bar(-4:1:4,b42,0.02);title('b42');axis([-16,16,-0.2,1.2]);
subplot(3,4,12),[h,f]=freqz(b42,1);
plot(f/pi,20*log10(abs(h)));axis tight;title('freq-b42');ylabel('');xlabel('');
pause;

%finally compare the three signals in sound quality
sound(y,fs); pause;
sound(x42,fs/4); pause;
sound(z42,fs); pause;


%write results to files
wavwrite(x21,fs/2,'x21.wav');
wavwrite(z21,fs,'z21.wav');
wavwrite(x41,fs/4,'x41.wav');
wavwrite(z41,fs,'z41.wav');
wavwrite(x22,fs/2,'x22.wav');
wavwrite(z22,fs,'z22.wav');
wavwrite(x42,fs/4,'x42.wav');
wavwrite(z42,fs,'z42.wav');

echo off;