How to generate smooth filtered envelope on EMG data in Matlab

I'm new to analysing EMG data and would appreciate some carefully explained help.

I would like to generate a smooth, linear enevelope signal of my EMG data (50kHz sampling rate) like the one published in this paper: https://openi.nlm.nih.gov/detailedresult.php?img=PMC3480942_1743-0003-9-29-3&req=4

My end goal is to be able to analyze the relationship between EMG activity (output) and action potentials fired from upstream neurons (putative input) recorded at the same time.

Even though this paper lists the filtering methods out quite clearly, I do not understand what they mean or how to perform them in matlab, which is the analysis tool I have available to me.

In the code I have written so far, I can dc offset as well as rectify my data:

x = EMGtime_data
y = EMGvoltage_data

%dc offset
y2=detrend(y)
% Rectification of the EMG signal
rec_y=abs(y2);
plot(x, rec_y)

But then I am not sure how to proceed. I have tried the envelope function, but it is not as smooth as I would like: For instance, if I used the following:

 envelope(y_rec,2000,'rms')

I get this (which also doesn't seem to care that the data is rectified):

Even if I were to accept the envelope function, I'm not sure how to access just the processed envelope data to adjust the plot (i.e. change the y-range), or analyse the data further for on-set and off-set of the signal since the results of this function seem to be coupled with the original trace.

I have also come across fastrms.m, which seems promising. Unfortunately, I do not understand how to implement this function since the general explanation is over my head and the example code is lacking any defined variable (so I don't know where to integrate my own data!) The example code from fastrms.m file exchange is here

Fs = 200; T = 5; N = T*Fs; t = linspace(0,T,N); 
noise = randn(N,1); 
[a,b] = butter(5, [9 12]/(Fs/2)); 
x = filtfilt(a,b,noise); 
window = gausswin(0.25*Fs); 
rms = fastrms(x,window,[],1); 
plot(t,x,t,rms*[1 -1],'LineWidth',2); 
xlabel('Time (sec)'); ylabel('Signal') 
title('Instantaneous amplitude via RMS')

I will be eternally grateful for help in understanding how to filter and smooth EMG data!

In order to analysis EMG signals in time domain, researcher use The combination of rectification and low pass filtering which is also called finding the “linear envelope” of the signal.
And as mentioned in both the above sentence and your attached article image's explanation, in order to plot overlaid signal, you could simply low pass filter your signal at specific frequency.
In your attached article the said signal was filtered at 8 HZ.

For better understanding the art of EMG signal analysis , i think this document could help you a lot (link)