生体信号の周波数解析 平均周波数について

筋電図の周波数解析をしています。

平均周波数を求めるために、以下のようなcodeで取り組んでみました。

2行n列のcsvデータで、1行目を解析対象としたものになります。

csvFilePath = 'ファイル名';

data = readmatrix(csvFilePath);

firstRowData = data(1, :);

fs = 1000;

fftData = fft(firstRowData);

powerSpectrum = abs(fftData).^2;

N = length(firstRowData);

frequencies = (0:(N-1)) * fs / N;

meanFrequency = sum(frequencies .* powerSpectrum) / sum(powerSpectrum);

varianceFrequency = sum(((frequencies - meanFrequency).^2) .* powerSpectrum) / sum(powerSpectrum);

fprintf('1行目のデータの平均周波数: %.2f Hz\n', meanFrequency);

fprintf('1行目のデータの周波数の分散: %.2f Hz^2\n', varianceFrequency);

上記codeになりますが、

添付している2つの筋電図の波形で、平均周波数が498、499Hz程度と差がでません。

他にも多くのデータを解析しましたが、どのような筋電図の波形を用いても497-499Hzにおさまってしまいます。

見た目がかなり異なる波形においても、平均周波数がほぼ同じになってしまうため、

codeが間違っているのではないかと考えております。

どうすれば解決するかご教授頂けないでしょうか。

Answer 1

[COVAO]

ナイキスト周波数(サンプリング周波数の1/2）を超える周波数のデータで平均を算出しているためと考えられます。

Signal Processing Toolboxのmeanfreqを使うと平均周波数を算出することができます。

上記コードでナイキスト周波数以下のデータを使用すると、meanfreqと近い値になるようです。

% Sample Data
nSamp = 1000;
Fs = 1000;
SNR = 40;
rng default
t = (0:nSamp-1)'/Fs;
x = chirp(t,50,nSamp/Fs,100);
x = x+randn(size(x))*std(x)/db2mag(SNR);
data = x';

%質問のコードで計算
firstRowData = data(1, :);
fs = 1000;
fftData = fft(firstRowData);
powerSpectrum = abs(fftData).^2;
N = length(firstRowData);
frequencies = (0:(N-1)) * fs / N;

meanFrequency = sum(frequencies .* powerSpectrum) / sum(powerSpectrum);
fprintf('1行目のデータの平均周波数: %.2f Hz\n', meanFrequency);

1行目のデータの平均周波数: 500.00 Hz

%ナイキスト周波数以下のデータを使用
powerSpectrum2 = powerSpectrum([1:N/2+1]);
frequencies2 = frequencies([1:N/2+1]);
meanFrequency2 = sum(frequencies2 .* powerSpectrum2) / sum(powerSpectrum2);
fprintf('1行目のデータの平均周波数: %.2f Hz\n', meanFrequency2);

1行目のデータの平均周波数: 75.02 Hz

%meanfreq関数で算出
meanfreq(data,fs)

ans = 75.0220