Remove artefact in signal, and then look at the PSD?

Hello everyone,

I have filtered my signal according to the frequencies I am interested in using butterworth. But there are some large amplitude artefacts in my signal that I want to remove. I have been advised replacing the signal with NaN, but I want to use pwelch to look at the signal's PSD. But pwelch wpn't work if I have a NaN. I've tried an autoregression spectrum, but I am not confident in using that. I don't really want to use a median filter, because it smooths out the data too much (I'm working with EEG data).

What other ways can I use to remove these types of artefacts so that I can look at the signal in the PSD plot?

Example code to replace artefacts with NaN:

N=38400;
M=128;
fs=128;
x=randn(N,1); % Create signal with the same parameters as my data

mi=10;
ma=20;
x(4032:4165)=(ma-mi).*rand(length(x(4032:4165)),1)+mi; % Add in large amplitude artefacts
x(25875:25989)=(ma-mi).*rand(length(x(25875:25989)),1)+mi;
t = (0:N-1)/fs; % Get time vector of x

figure
plot(t,x); % Look at where artefacts are located

% First blip in time domain
x1 = fs*31.5 % X axis units (time domain) (The units change everytime in x because the signal is produced with randn, so wen you do it, you may not get the same units)
x2 = fs*32.53

x3 = t(4032:4165);
x3 = (x3*fs)+1; % Find corresponding indices on the signal (x) that relate to the artefacts on the time vector

% Second blip in time domain
y1 = fs*202.1 % X axis units (time domain)
y2 = fs*203

y3 = t(25875:25989);
y3 = (y3*fs)+1;

x(x3) = NaN; % Replace artefacts with NaN
x(y3) = NaN;