Matlab经典功率谱估计——ECG功率谱

---

前言

经典谱估计又称非参数估计，其将数据工作区外的未知数据假设为零，相当于数据加窗。本篇文章介绍了基于ECG（心电信号）的谱估计。主要用到的方法为周期图法、改进的周期图法、Welch法和MTM法。

---

一、相关介绍

周期图法又称直接法，它从信号序列中截取一段，将其视为有限信号序列真实功率谱的估计的抽样。而改进的周期图法在进行快速傅里叶前进行加窗，平滑了边缘，降低了旁瓣的高度。Welch法将信号序列划分成几段，然后采用改进的周期图法，并取平均。MTM法没有使用带通滤波器，而是使用了最优滤波器进行计算估计。

二、功率谱绘制

1.周期图法

代码如下：

clear
clc
load data.txt
ecg=data;
Fs=2000;
N=1024;Nfft=512;
%数据的长度和FFT所用的数据长度
n=0:N-1;t=n/Fs;
Pxx=10*log10(abs(fft(ecg,Nfft).^2)/N);
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);
plot(f,Pxx)
xlabel('频率/Hz');ylabel('功率谱/dB');
title('周期图法');

2.改进的周期图法

代码如下：

clear
clc
load data.txt
ecg=data;
Fs=2000;
N=1024;Nsec=256;
%数据的长度和FFT所用的数据长度
n=0:N-1;t=n/Fs;%采用的时间序列
%四段功率谱Pxx1-Pxx4
Pxx1=abs(fft(ecg(1:256),Nsec).^2)/Nsec;
Pxx2=abs(fft(ecg(257:512),Nsec).^2)/Nsec;
Pxx3=abs(fft(ecg(513:768),Nsec).^2)/Nsec;
Pxx4=abs(fft(ecg(769:1024),Nsec).^2)/Nsec;
Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4/4);%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);
subplot(211),plot(f(1:Nsec/2),Pxx(1:Nsec/2));
xlabel('频率/Hz');ylabel('功率谱/dB');
title('平均周期图（无重叠） N=4*256');grid on;
%运用信号重叠分段估计功率谱
Pxx1=abs(fft(ecg(1:256),Nsec).^2)/Nsec; 
Pxx2=abs(fft(ecg(129:384),Nsec).^2)/Nsec;
Pxx3=abs(fft(ecg(257:512),Nsec).^2)/Nsec;
Pxx4=abs(fft(ecg(385:640),Nsec).^2)/Nsec;
Pxx5=abs(fft(ecg(513:768),Nsec).^2)/Nsec;
Pxx6=abs(fft(ecg(641:896),Nsec).^2)/Nsec;
Pxx7=abs(fft(ecg(769:1024),Nsec).^2)/Nsec;
Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7/7);
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);
subplot(212),plot(f(1:Nsec/2),Pxx(1:Nsec/2));
xlabel('频率/Hz');ylabel('功率谱/dB');
title('平均周期图（重叠1/2） N=1024');

3.Welch法

代码如下：

clear
clc
load data.txt
ecg=data;
Fs=2000;%采样频率
NFFT=1024;
t=0:1/Fs:1;
window=boxcar(100);%采用矩形窗，其他窗可自行验证
noverlap=20;%指定段与段之间的重叠的样本数
[pxx,f]=pwelch(ecg,window,noverlap,NFFT,Fs);
pxx=10*log10(pxx);
plot(f,pxx);
title('矩形窗');

4.MTM法

代码如下：

clear
clc
load data.txt
ecg=data;
Fs=2000;
N=1024;Nfft=256;n=0:N-1;t=n/Fs;
[Pxx,f]=pmtm(ecg,4,Nfft,Fs);
plot(f,10*log10(Pxx));
xlabel('频率/Hz');ylabel('功率谱/dB');
title('多窗口法（MTM）NW=4');

---

三、运行结果

1.周期图法

2.改进周期图法

3.Welch法

4.MTM法

最后

如果你想要进一步了解更多的相关知识，可以关注下面公众号联系~会不定期发布相关设计内容包括但不限于如下内容:信号处理、通信仿真、算法设计、matlab appdesigner，gui设计、simulink仿真…希望能帮到你！