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

31 ビュー (過去 30 日間)
HelpAStudent
HelpAStudent 2023 年 6 月 17 日
回答済み: Alan 2023 年 9 月 7 日
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;

回答 (1 件)

Alan
Alan 2023 年 9 月 7 日
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

カテゴリ

Help Center および File ExchangeSpectral Analysis についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by