Python和Matlab中平均频率的差异
Python和Matlab中平均频率的差异
提问于 2016-06-20 12:44:17
我有这个肌电信号，我想在这个文章的基础上画出平均工频。我使用以下代码在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中做同样的事情。守则是：
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
df = pd.read_csv('EMG.txt', names=['EMG'])
EMG=df['EMG'].tolist()
sampling_rate=200
N=1000 #my window
FT=np.fft.fft(EMG, axis=0)
psd=FT*np.conj(FT)
NFFT=len(FT)
f =(np.arange(0,NFFT/2)*sampling_rate)/N
NUM=0
DUM=0
MPF=[]
for j in np.arange(1,NFFT/2):
NUM=NUM+f[j]*psd[j]
DUM=DUM+psd[i]
MPF.append(NUM/DUM)
plt.plot(MPF)
plt.show() 强积金的计划是：
为什么不一样？
更新
按照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个样本，结果如下：
回答 2
Stack Overflow用户
发布于 2020-12-30 17:37:14
“下降趋势”选项将导致第一个值(0 Hz或DC)的差异。
默认设置为“常量”，它从原始数据中删除DC组件。
f, Pxx_den = signal.periodogram(x, fs, detrend=False)Stack Overflow用户
发布于 2019-06-07 02:24:37
我正在继续这个问题和答案。
1)在Matlab中，尝试使用内置函数meanfreq，medfreq也可以帮助您找到中值频率。
2)到目前为止，我还没有找到一个包来解决Python的平均/中位频率。
3)我还尝试用Python来求解平均/中位频率。并遵循Matlab meanfreq中的计算步骤。(在Matlab命令窗口中，键入edit meanfreq.m，您可以看到meanfreq的源代码。源代码会有很大帮助。
4)以下是基于python的meanfreq实现
from scipy import signal
def meanfreq(x, fs):
f, Pxx_den = signal.periodogram(x, fs)
Pxx_den = np.reshape( Pxx_den, (1,-1) )
width = np.tile(f[1]-f[0], (1, Pxx_den.shape[1]))
f = np.reshape(f, (1, -1))
P = Pxx_den * width
pwr = np.sum(P)
mnfreq = np.dot(P, f.T)/pwr
return mnfreq5)可以同时调试of和Matlab的源代码。
6)我的代码中也有一个小的'bug‘，或者不是我的。调试它们时，您可能会发现Pxx的第一个值(只有第一个值)完全不同。其余的值是相同的。，我在这里也很困惑。 :)
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/37922928
复制相关文章
相似问题
腾讯云开发者
