我有这个EMG signal，我想根据这个article绘制平均功率频率。我使用以下代码在Matlab中实现它：clear all;

close all;

EMG=load('EMG.txt');

N=1000; %my window

z=1;

fs=200 %sampling rate

for i=1:length(EMG)-N

DUM=0;

NUM=0;

FT=fft(EMG(i:i+N-1));

psd=FT.*conj(FT);

NFFT=length(fft2);

f = [1:NFFT/2]*fs/N;

for j=1:NFFT/2

NUM=NUM+f(j)*psd(j);

DUM=DUM+psd(j);

end

MPF(z)=NUM/DUM;

z=z+1;

end

强积金的情节是：

下面我尝试在Python中做同样的事情。代码是：

^{pr2}$

强积金的地块是：

为什么不同？在

更新

根据Dan在评论部分的建议，我修改了Python代码如下，结果大致相同，只是Matlab代码比Python快得多，在我的例子中，Python内存不足：sampling_rate=200

N=1000

MPF=[]

for i in range(0,len(EMG)-N):

signal=EMG[i:(i+N)]

FT=np.fft.fft(signal, axis=0)

psd=FT*np.conj(FT)

NFFT=len(FT)

f =(np.arange(0,NFFT/2)*sampling_rate)/N

D_1=0

N_1=0

for j in np.arange(1,NFFT/2):

D_1=D_1+f[j]*psd[j]

N_1=N_1+psd[j]

MPF.append(D_1/N_1)

plt.plot(MPF)

plt.show()

选择前22000个样本，结果如下：