function esim1eng
% ESIM1ENG Demo on LTI filtering
%
% Let us first hear a signal with noise, and
% then filter it so that only the frequency band
% of wired phone (< 3300 Hz) is left.

% Jukka Parviainen, parvi@hut.fi, 12.9.2002 / 17.9.2004
% Tiedostot saatavilla 
% http://www.cis.hut.fi/Opinnot/T-61.246/Demo/esim1.shtml
%
% Alkutoimenpiteet: tarkista että esityskoneessasi on 
% kaiuittimet tai salin kyseessä ollen audioliitäntä.
% Lataa tiedostot työhakemistoon tai aseta Matlabissa
% addpath-komennolla polku tunnetuksi.
% Kutista Matlab-työikkuna suhteellisen pieneksi ruudun alareunaan,
% sijoita piirtoikkuna, esimerkiksi ensimmäisen vaiheen jälkeen,
% sopivaan paikkaan, että se näkyy.
%
% Tiedostot: esim1.m, esim1.wav [esim1eng.m, esim1eng.wav]
%
% HUOM! Demoluonteen takia suurin osa _oikeasta_ laskennasta
% tehdään aivan koodin alussa.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;      % shut all Matlab windows
clc;            % clears the screen

fprintf(1,'****************************************\n');
fprintf(1,'****************************************\n');
fprintf(1,'************  esim1eng.m  **************\n');
fprintf(1,'****************************************\n');
fprintf(1,'****       Filtering  Demo           ***\n');
fprintf(1,'****************************************\n');
fprintf(1,'\n\n');

fprintf(1,'** Wait a moment...\n\n');

% Read the wav file into Matlab
%[y, fs, nbits] = wavread('esim1eng.wav'); 
[y, fs, nbits] = wavread('esim1.wav'); 

Wp = 2*3000/fs; % End of passband
Ws = 2*3300/fs; % Start of stopband
Rp = 0.2;       % Ripple in passband
Rs = 40;        % Ripple in stopband
[N,Wn] = ellipord(Wp, Ws, Rp, Rs);
[b,a] = ellip(N, Rp, Rs, Wn);
[hh,ww]=freqz(b,a,1024,fs);   % amplitude- ja phase response

Y = fft(y,2^14);        % Compute spectrum form first 2^14 points

ysuod = filter(b,a,y);  % Filtering for those

Ysuod = fft(ysuod,2^14); % And correspondingly.

fprintf(1,'** Press ENTER...\n\n');
pause;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

fprintf(1,'Let us hear a signal with noise.\n\n');

t = linspace(0, length(y)/fs, length(y));

soundsc(y, fs);  %Kuunnellaan

fprintf(1,'** Press ENTER after listening...\n\n');
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

w = linspace(0, fs, length(Y));
subplot(311);
plot(w, abs(Y));
axis([0 fs/2 0 1.1*max(abs(Y))]);
grid on;
title(['Signal spectrum (temporary) 0..' num2str(fs/2) ' Hertz']);
xlabel('frequency (Hz)');

fprintf(1,'The spectrum shows which frequency components are big in the signal.\n\n');

fprintf(1,'** Press ENTER...\n\n');
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

fprintf(1,'There's a lot of high-frequency noise in the signal.\n\n');

fprintf(1,'Let us realize a lowpass filter ("wired phone").\n\n');
fprintf(1,'Quick changes in the signal in time-domain\ncorrespond high frequencies.\n');
fprintf(1,'Slow changes in the signal in time-domain\ncorrespond low frequencies.\n');

fprintf(1,'** Press ENTER...\n\n');
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Seuraava pelkkää pluffia
h = waitbar(0,'Computing the filter coefficients...');
for i=1:20,
  pause(0.05);
  waitbar(i/20,h)
end
close(h)

subplot(312);
ku = plot(ww,abs(hh).^2);
axis([0 ww(end) 0 1.05]);
set(ku,'LineWidth',2);
grid on;
title('Elliptic lowpass filter H(z)');
subplot(311); xlabel('');
subplot(312); xlabel('frequency (Hz)');

fprintf(1,'An elliptic lowpass filter is implemented (middle).\n');
fprintf(1,'Only low frequencies of the signal will remain.\n\n');

fprintf(1,'** Press ENTER...\n\n');
pause;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

fprintf(1,'Now let us filter the noisy signal...\n\n');

fprintf(1,'\n Convolution theorem: \n convolution in time-domain corresponds product in  frequency-domain\n');
fprintf(1,'\n');
fprintf(1,'                 --------------\n');
fprintf(1,'                 |            |\n');
fprintf(1,'  X(e^jw)  ----> |  H(e^jw)   | ---->  Y(e^jw) = H(e^jw) . X(e^jw)\n');
fprintf(1,'                 |            |\n');
fprintf(1,'                 --------------\n');
fprintf(1,'\n');
fprintf(1,'\n');

% Seuraava pelkkää pluffia
h = waitbar(0,'Filtering...');
for i=1:40,
  pause(0.03);
  waitbar(i/40,h)
end
close(h)

w = linspace(0, fs, length(Ysuod));
subplot(313);
plot(w, abs(Ysuod));
grid on;
axis([0 fs/2 0 1.1*max(abs(Y))]);
title('Spectrum of the filtered signal');
subplot(312), xlabel('');
subplot(313), xlabel('frequency (Hz)');

fprintf(1,'Let us hear the filtered signal.\n\n');

fprintf(1,'** Press ENTER...\n\n');
pause;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

soundsc(ysuod,fs);

fprintf(1,'Does the filtered signal sound "better"?\n\n');

fprintf(1,'** Press ENTER after listening...\n\n');
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Mitä on tapahtunut aikatasossa?!

clc

fprintf(1,'Let us plot the original and filtered\n signal in time-domain\n\n');

figure(2), clf;
alku = 44800;
loppu = 45000;
plot(t(alku:loppu), y(alku:loppu));  
hold on;
ku3 = plot(t(alku:loppu), ysuod(alku:loppu), 'r');
set(ku3,'LineWidth',2);
title('Original and filtered signal x(t)');
xlabel('time (seconds)');
ylabel('amplitude');
legend({'Original','Filtered'});
grid on;
axis([t(alku) t(loppu) 1.1*min(y(alku:loppu)) 1.1*max(y(alku:loppu))]);


fprintf(1,'** THE END. Press ENTER.\n\n');
pause;

% If you need to write the signal back to WAV format:
%wavwrite(ysuod, fs, nbits, 'esim1_output.wav');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
