%% Multirate: Resampling kiisseli.wav
%
% This demo shows how upsampling (ylösnäytteitys) by
% factor L and downsampling (alasnäytteistys) by factor M
% are applied when changing the sampling rate (frequency).
%
%                  x2              x3
%           -----        --------        -----
%  x[n] --> | L | -----> | H_LP | -----> | M | --> y[n]
%           -----        --------        -----
%
%
% The demo gives you four figures. Each figure
% contains a time-domain waveform in seconds and a frequency-domain
% spectrum from 0 .. fs/2 Hz.
%
%   Fig 1  -- original signal x
%   Fig 2  -- upsampled signal x2
%   Fig 3  -- upsampled signal after lowpass filtering x3
%   Fig 4  -- downsampled signal y  <==> result
%
% Anti-imaging (interpolation) filter and anti-aliasing
% filter are combined to a single LP filter with cut-off
% frequency of min(w_M, w_L).
%

% Jukka Parviainen, 28.4.2009 / 25.4.2008

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

%% Variables to be edited
%
% YOU CAN EDIT THESE VARIABLES HERE.
% DO NOT CHANGE ANYTHING ELSE.

% Audiofile to be downloaded and analysed
audiofile = 'kiisseli.wav';

% Upsampling (L) and downsampling (M) factors
L    = 3;
M    = 5;

% Is the lowpass filter in Step 3 in use?
% filterInUse = 0;
filterInUse = 1;

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

%% Shut down figures and remove some variables
%
% NOTE! It is probably useful to say
%   clear all  % clear all variables
%   close all  % shut down all windows
% before starting this demo!

close all;
clear x x2 x3 y 

%% Download the audio file
[x, fs] = wavread(audiofile);
%soundsc(x, fs)
disp(['Audiofile                        ' audiofile]);
disp(['#samples                       = ' num2str(length(x))]);
disp(['Original sampling frequency fs = ' num2str(fs) ' Hz']);
disp(['Original sampling period    Ts = ' num2str(1/fs) ' s']);
disp(['Length in seconds              = ' num2str(length(x)/fs) ' s']);


%% Set the upsampling and downsampling factors
fsL  = L*fs;
fsLM = L*fs/M;
disp(' ');
disp(['Resampling by factor             L/M = ' num2str(L) '/' num2str(M)]);
disp(['Sampling freq. after upsampling  fsL = ' num2str(fsL) ' Hz']);
disp(['Resampled sampling frequency    fsLM = ' num2str(fsLM) ' Hz']);
disp(['Sampling per. after resampling  TsLM = ' num2str(1/fsLM) ' s']);

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

%% Fig 1 -- Original signal
P1  = length(x);
t   = [0 : P1-1]/fs;
xF  = fft(x);
w   = fs*[0 : P1-1]/P1;

figure(1); clf;
subplot(211);
plot(t, x);
grid on;
xlabel('time');
title(['Original signal ' audiofile ' [#samples = ' num2str(P1) ']']);

subplot(212);
plot(w, 20*log10(abs(xF)));
set(gca,'XLim',[0 fs/2]);
grid on;
xlabel(['frequency (fs = ' num2str(fs) ' Hz)']);

% some lines and text:
dummyY = get(gca,'YLim');
pO     = line([fs/2 fs/2], dummyY);
pF     = line([fsLM/2 fsLM/2], dummyY);
set(pO, 'Color', [1 0 0], 'LineWidth', 3);
set(pF, 'Color', [0 1 0], 'LineWidth', 3);

dummyF = dummyY(2) - 0.07*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fs/2), dummyF, '\omega = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

dummyF = dummyY(2) - 0.14*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fsLM/2), dummyF, '\omega_L_M = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

%% BREAKPOINT
%
% - What is sampling frequency?

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

%% Fig 2 -- Upsampling by L
P2  = L*(P1-1)+1;
t2  = [0 : P2-1]/fsL;
x2(1:L:P2) = x;
x2F = fft(x2);
w2  = fsL*[0 : P2-1]/P2;

figure(2); clf;
subplot(211);
plot(t2, x2);
grid on;
xlabel('time');
title(['Upsampled by L = ' num2str(L) ' [#samples = ' num2str(P2) ']']);

subplot(212);
plot(w2, 20*log10(abs(x2F)));
set(gca,'XLim',[0 fsL/2]);
grid on;
xlabel(['frequency (fs = ' num2str(fsL) ' Hz)']);

% some lines and text:
dummyY = get(gca,'YLim');
pO     = line([fs/2 fs/2], dummyY);
pF     = line([fsLM/2 fsLM/2], dummyY);
pT     = line([fsL/2 fsL/2], dummyY);
set(pO, 'Color', [1 0 0], 'LineWidth', 3);
set(pF, 'Color', [0 1 0], 'LineWidth', 3);
set(pT, 'Color', [1 1 0], 'LineWidth', 3);

dummyF = dummyY(2) - 0.07*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fs/2), dummyF, '\omega = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

dummyF = dummyY(2) - 0.14*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fsLM/2), dummyF, '\omega_L_M = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

dummyF = dummyY(2) - 0.21*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fsL/2), dummyF, '\omega_L = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

%% BREAKPOINT
%
% - What is the new sampling frequency f_L? 
% - Confirm that there are L-1 zeros between
%   each original sample of x in waveform

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

if filterInUse

% Fig 3 -- Lowpass filtering
% Removing "images" and anti-aliasing
w_cut   = 2*min([fsL/(2*L) fsL/(2*M)])/fsL; % ???? does it work????
f_cut   = w_cut*fsL/2;
Rp      = 1;
Rs      = 50;
[N, Wn] = ellipord(0.92*w_cut, 0.98*w_cut, Rp, Rs); 
[B, A]  = ellip(N, Rp, Rs, Wn);
x3      = filter(B, A, x2);
x3F     = fft(x3);

figure(3); clf;
subplot(211);
plot(t2, x3);
grid on;
xlabel('time');
title(['Up by L = ' num2str(L) ' and FILTERED with f_c \approx ' num2str(f_cut) ' Hz']);

subplot(212);
plot(w2, 20*log10(abs(x3F)));
set(gca,'XLim',[0 fsL/2]);
grid on;
xlabel(['frequency (fs = ' num2str(fsL) ' Hz)']);

% some lines and text:
dummyY = get(gca,'YLim');
pO     = line([fs/2 fs/2], dummyY);
pF     = line([fsLM/2 fsLM/2], dummyY);
pT     = line([fsL/2 fsL/2], dummyY);
set(pO, 'Color', [1 0 0], 'LineWidth', 3);
set(pF, 'Color', [0 1 0], 'LineWidth', 3);
set(pT, 'Color', [1 1 0], 'LineWidth', 3);

dummyF = dummyY(2) - 0.07*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fs/2), dummyF, '\omega = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

dummyF = dummyY(2) - 0.14*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fsLM/2), dummyF, '\omega_L_M = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

dummyF = dummyY(2) - 0.21*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fsL/2), dummyF, '\omega_L = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

else
  % Not filtering the signal at all! 
  % For (aliasing) demo purposes only.
  warning('No anti-aliasing / anti-imaging filter used!');
  x3 = x2;
end;

%% BREAKPOINT
%
% - What does lowpass filtering affect?
% 

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

%% Fig 4 -- Downsampling by M
y    = x3(1:M:end);
P4   = length(y);
t4   = [0 : P4-1]/fsLM;

yF   = fft(y);
w4   = fsLM*[0 : P4-1]/P4;


figure(4); clf;
subplot(211);
plot(t4, y);
grid on;
xlabel('time');
if filterInUse
  title(['Up L = ' num2str(L) ' AND filtered AND downsampled by M = ' num2str(M)]);
else
  title(['Up L = ' num2str(L) ' (NOT filtered) AND downsampled by M = ' num2str(M)]);
end;

subplot(212);
plot(w4, 20*log10(abs(yF)));
set(gca,'XLim',[0 fsLM/2]);
grid on;
xlabel(['frequency (fs = ' num2str(fsLM) ' Hz)']);

% some lines and text:
dummyY = get(gca,'YLim');
pO     = line([fs/2 fs/2], dummyY);
pF     = line([fsLM/2 fsLM/2], dummyY);
pT     = line([fsL/2 fsL/2], dummyY);
set(pO, 'Color', [1 0 0], 'LineWidth', 3);
set(pF, 'Color', [0 1 0], 'LineWidth', 3);
set(pT, 'Color', [1 1 0], 'LineWidth', 3);

dummyF = dummyY(2) - 0.07*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fs/2), dummyF, '\omega = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

dummyF = dummyY(2) - 0.14*(dummyY(2)-dummyY(1));
tF     = text(0.99*(fsLM/2), dummyF, '\omega_L_M = \pi \rightarrow');
set(tF,'HorizontalAlignment','right');

%% BREAKPOINT
%
% - Listen to the resampled signal y (if possible)

%%
%soundsc(x, fs);     % original
%%
%soundsc(y, L/M*fs); % resampled by L/M

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

%% Summary
%
% In order to change sampling rate
%  (1) upsample
%  (2) apply lowpass filtering
%  (3) downsample
% 
% _All_ this can be done in Matlab with a single command
%
% y = resample(x, L, M);
%