How to calculate Power Spectral Density using Welch, Periodogram, Yule-Walker, and Burg Methods in ecg?

Hello! I have a signal that has been filtered using a specific code. I'm looking for assistance in calculating the power spectral density using various methods in MATLAB. Specifically, I would like to utilize the Welch method and the periodogram, as well as explore the Parametric estimation of the spectrum using the Yule-Walker and/or Burg method. Could someone please guide me through the process?

load sig_ecg.mat

Fs = 300;
L = numel(sig);
t = linspace(0, L-1, L)/Fs;
% t=linspace(0,length(sig)/fs,length(sig));
figure()
plot(t,sig/max(sig));
grid
title('Original Signal in Time Domain');
xlabel('Time(s)');
ylabel('Amplitude');

figure()
sig2=sig-mean(sig); % Per visualizzare meglio lo spettro
f=linspace(-Fs/2,Fs/2,length(sig2));
plot(f,fftshift(abs(fft(sig2))));
grid
title('Signal in Frequency Domain'); 
xlabel('Frequency (Hz)');
ylabel('Amplitude');

Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
NFFT = 2^nextpow2(L)

sig = bandstop(sig, [48 52], Fs, 'ImpulseResponse','iir');  % Remove 50 Hz Mains Interference
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Ws = 0.5/Fn;                                                % Passband Frequency Vector (Normalised)
Wp = 2.5/Fn;                                                % Stopband Frequency Vector (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs = 150;                                                   % Stopband Attenuation (dB)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                             % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp,'high');                         % Default Here Is A Lowpass Filter
[sos,g] = zp2sos(z,p,k);                                    % Use Second-Order-Section Implementation For Stability
sig_filtered = filtfilt(sos,g,sig);                           % Filter Signal (Here: isrg)
figure
freqz(sos, 2^14, Fs)                                        % Bode Plot Of Filter
set(subplot(2,1,1), 'XLim',[0 50])                          % Optional, Change Limits As Necessary
set(subplot(2,1,2), 'XLim',[0 50])                          % Optional, Change Limits As Necessary

Below there is the code where I tried to calculare the power spectral density and so on:

%       WELCH METHOD POWER SPECTRAL DENSITY
[wel,f_wel] = pwelch(ecg, hamming(1000),250,1000,Fs); % Finestra Hamming di 1000 campioni con overlap di 1/4 e no zero padding

figure()
plot(f_wel,wel)
title('PSD Welch');
%       PERIODOGRAM
[per_rect, f_per_rect] = periodogram(ecg,[],1000,fs);
plot(f_per_rect,per_rect);
title('PSD Periodogram with Rectangular Window')
figure;

[per_ham, f_per_ham] = periodogram(ecg,hamming(length(ecg)),1000,fs);
plot(f_per_ham,per_ham);
title('PSD Periodogram with Hamming Window')
figure;

%       Parametric estimation of the spectrum
% Yule-Walker
order=50; 
ar_yule=zeros(order,order+1);
for i=1:order
    [ar_yule(i,1:i+1),e_ord(i), rc_ord]= aryule(ecg,i);
end 

% a = autoregressive parameters 
% e = variance 
% rc = coefficent of riflection

% Yule
figure()
subplot(1,2,1)
stem(rc_ord);
title('Reflection coefficients')
subplot(1,2,2)
plot(e_ord);
title('Variance');

ordine = 9; % The optimal order is chosen by the behavior of the reflection coefficients (they should go to zero) and of the variance.
[Pxx_y,f_y] = pyulear(ecg,ordine,fs);
figure()
plot(f_y*(fs/2)/pi,Pxx_y); % This is done because the output frequencies from pyulear are normalized frequencies.
title('Yule-Walker')
xlabel('Frequency [Hz]'); 
ylabel('dB/Hz');

% Burg
ord = 50;
ar_burg=zeros(ord,ord+1);
for i=1:ord
    [ar_burg(i,1:i+1),eb_ord(i),rcb_ord]=arburg(ecg,i);
end
figure()
subplot(1,2,1)
stem(rcb_ord);
title('Reflection coefficients');
subplot(1,2,2)
plot(eb_ord);
title('Variance');

ordine_b = 9;
figure()
[Pxx_b,f_b]=pburg(ecg,ordine_b,fs);
plot(f_b*(fs/2)/pi,Pxx_b,'r');
title('Burg')
xlabel('Frequency [Hz]'); 
ylabel('dB/Hz');
hold on;

[HelpAStudent]

Sorry I forgot the signal

Answer 1

[Alan]

Hi,

From what I understand, you are exploring different methods to calculate the Power Spectral Density and conduct Parametric Estimation of the Spectrogram of ECG signals. You can find explanations and examples of Welch and periodogram in MATLAB in the following documentation pages:

1. Welch’s estimate: https://www.mathworks.com/help/signal/ref/pwelch.html
2. Periodogram: https://www.mathworks.com/help/signal/ref/periodogram.html

There are various methods to do parametric estimation of the spectrum including the Yule-Walker and Burg’s methods in the following documentation page: https://www.mathworks.com/help/signal/parametric-spectral-estimation.html

For detailed debugging assistance you could contact MathWorks technical support: https://www.mathworks.com/support/contact_us.html