% EE 438 Digital Signal Processing with Applications - Spring 1997
% Application of elementary signal operations to speech signal
%---------------------------------------------------------------
% Initialize
clear all;
close all;
%---------------------------------------------------------------
% Parameters
D = 4
%---------------------------------------------------------------
%Read and play back data sampled at 8192kHz
Fs = 8192;
data=load('erf1s1t0');

%plot(data)
%end1=input('end1?');
%end2=input('end2?');

x=data;
x(8192)=0;	% zero-pad if length(data) < 8192
x=x(1:8192);
N = length(x)

%---------------------------------------------------------------
%Design LPF with cutoff at Fs/D
N_spect = 128;
del_mu = 1./N_spect;
mu = -0.5:del_mu:0.5-del_mu;
N_filt = 64;
delta = 0.02;
bands = [0.0, 1./D-delta, 1./D+delta, 1.0];
H_ideal = [1.0, 1.0, 0.0, 0.0];
h = remez(N_filt, bands, H_ideal);
%---------------------------------------------------------------
%Prepare plots
figure(1);
H = fft(h,N_spect);
H = fftshift(H);
Hmag = abs(H);
Hpha = angle(H);
subplot(3,1,1), stem(h);
ylabel('Filter');
subplot(3,1,2), plot(mu,Hmag);
ylabel('Magnitude');
subplot(3,1,3), plot(mu,Hpha);
ylabel('Phase');
pause;
%===============================================================
%Decimation = LPF followed by downsampler
%===============================================================
%Filter signal
xbl = conv(h, x);
% Downsample by D
n = 1:1:N;
t = (n-1)./Fs;
y = zeros(1,N);
y(1:N/D) = xbl(1:D:N);
y(N/D+1:N) = zeros(1,N-N/D);

figure(2);
subplot(2,1,1), plot(t,x);
xlabel('t sec');
title('Original utterance');
subplot(2,1,2), plot(t,y);
xlabel('t sec');
title('Utterance decimated by D');

X = fft(x);
X = fftshift(X);
Xmag = abs(X);
Y = fft(y);
Y = fftshift(Y);
Ymag = abs(Y);
delta_f = Fs./(N.*1000);
nf = -N./2:1:N/2-1;
f = nf .* delta_f;

figure(3);
subplot(2,1,1), plot(f,Xmag);
xlabel('f kHz');
title('Original utterance');
subplot(2,1,2), plot(f,Ymag);
xlabel('f kHz');
title('Utterance decimated by D');

input('Original utterance')
soundsc(x,Fs);
input('Utterance decimated by D, played at Fs/D')
soundsc(y,Fs./D);
%===============================================================
%Interpolation = upsampler followed by LPF
%===============================================================
% Upsample by D
z = zeros(1,N);
z(1:D:N) = y(1:N/D);
%Filter signal
zbl = conv(h, z);
z = zbl(1:N);

figure(2);
subplot(2,1,1), plot(t,x);
xlabel('t sec');
title('Original utterance');
subplot(2,1,2), plot(t,z);
xlabel('t sec');
title('Utterance interpolated by D');

Z = fft(z);
Z = fftshift(Z);
Zmag = abs(Z);
delta_f = Fs./(N.*1000);
nf = -N./2:1:N/2-1;
f = nf .* delta_f;

figure(3);
subplot(2,1,1), plot(f,Xmag);
xlabel('f kHz');
title('Original utterance');
subplot(2,1,2), plot(f,Zmag);
xlabel('f kHz');
title('Utterance interpolated by D');

input('Original utterance')
soundsc(x,Fs);
input('Utterance interpolated by D')
soundsc(z,Fs);