%% Read in sound sample
[y,FSamp] = audioread('Ex1_Mary.wav');
%% Extract wholeband Hilbert envelope, low-pass filter and plot
env = abs(hilbert(y));
cutoff = 50; % cutoff for low-pass filter, typically 40 Hz or 50 Hz
[a,b] = bilinear([0 0 0 15], [1 6 15 15], FSamp/(2*pi*cutoff)); % filter coefficients
LPenv = filtfilt(a,b,env); % perform zero-phase filtering (no phase distortion)
figure,
plot([1:length(y)]/FSamp,y,'k'),hold on
plot([1:length(y)]/FSamp,LPenv,'r','LineWidth',2)
axis tight
xlabel('Time (s)','FontName','Times','FontSize',18,'FontWeight','bold')
ylabel('Amplitude','FontName','Times','FontSize',18,'FontWeight','bold')
title('Wholeband Amplitude Envelope','FontName','Times','FontSize',24,'FontWeight','bold')