%% Effect of Input Padding on Perfect Reconstruction
%

% Copyright 2018 The MathWorks, Inc.

%% 
% To enable perfect reconstruction, the |mdct| function zero-pads the front
% and back of the audio input signal. The signal returned from |imdct|
% removes the zero padding added for perfect reconstruction.
%
% Read in an audio file, create a 2048-point Kaiser-Bessel-derived window,
% and then clip the audio signal so that its length is a multiple of 2048.

[x,fs] = audioread('Click-16-44p1-mono-0.2secs.wav');
win = kbdwin(2048);

xClipped = x(1:end - rem(size(x,1),numel(win)));

%%
% Convert the signal to the frequency domain, and then reconstruct it back
% in the time domain. Plot the original and reconstructed signals and
% display the reconstruction error.

C = mdct(xClipped,win);
y = imdct(C,win);

figure(1)
t = (0:size(xClipped,1)-1)'/fs;
plot(t,xClipped,'bo',t,y,'r.')
legend('Original Signal','Reconstructed Signal')
title(strcat("Reconstruction Error = ",num2str(mean((xClipped-y).^2))))
xlabel('Time (s)')
ylabel('Amplitude')

%%
% You can perform the MDCT and IMDCT without input padding using the
% |PadInput| name-value pair. However, there will be a reconstruction error
% in the first half-frame and last half-frame of the signal.

C = mdct(xClipped,win,'PadInput',false);
y = imdct(C,win,'PadInput',false);

figure(2)
t = (0:size(xClipped,1)-1)'/fs;
plot(t,xClipped,'bo',t,y,'r.')
legend('Original Signal','Reconstructed Signal')
title(strcat("Reconstruction Error (Without Input Padding) = ",num2str(mean((xClipped-y).^2))))
xlabel('Time (s)')
ylabel('Amplitude')

%%
% If you specify an input signal to the |mdct| that is not a multiple of
% the window length, then the input signal is padded with zeros. Pass the
% original unclipped signal through the transform pair and compare the
% original signal and the reconstructed signal.

C = mdct(x,win);
y = imdct(C,win);

figure(3)

subplot(2,1,1)
plot(x)
title('Original Signal')
ylabel('Amplitude')
axis([0,max(size(y,1),size(x,1)),-0.5,0.5])

subplot(2,1,2)
plot(y)
title('Reconstructed Signal')
xlabel('Time (s)')
ylabel('Amplitude')
axis([0,max(size(y,1),size(x,1)),-0.5,0.5])

%%
% The reconstructed signal is padded with zeros at the back end. Remove the
% zero-padding from the reconstructed signal, plot the original and
% reconstructed signal, and then display the reconstruction error.

figure(4)
y = y(1:size(x,1));
t = (0:size(x,1)-1)'/fs;
plot(t,x,'bo',t,y,'r.')
legend('Original Signal','Reconstructed Signal')
title(strcat("Reconstruction Error = ",num2str(mean((x-y).^2))))
xlabel('Time (s)')
ylabel('Amplitude')