Subtract the mean of the EKG data before calculating the fft. Also, it likely needs to be rescaled, because the amplitude should be between ±1 mV.
The signal needs to be filtered. Some of the noise is broadband, however much can be filtered with a frequency-selective filte —
LD = load('Sebastian Daneli ecg_stress_test.mat')
t = linspace(0, 1, L)/Fs;
[EKGfilt,df] = bandpass(EKG, [1 25], Fs, 'ImpulseResponse','iir');
FTEKGfilt = fft(EKGfilt,NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
[ampmax,idx] = max(abs(FTEKGfilt(Iv))*2);
fprintf('Amplitude = %.2f mV at %.2f Hz',ampmax,Fv(idx))
The pwelch function is probably the mose frequently used method to determine the power spectral density. To see the the signal spectrum and power spectrum as a function of time, use the pspectrum or spectrogram functions.