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

10 次查看（过去 30 天）
cho hunseki
cho hunseki 2023-10-8
编辑： covao 2023-11-25
筋電図の周波数解析をしています。
平均周波数を求めるために、以下のような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が間違っているのではないかと考えております。
どうすれば解決するかご教授頂けないでしょうか。

请先登录，再回答此问题。

回答（1 个）

covao
covao 2023-11-25
编辑：covao 2023-11-25
Ran in:
ナイキスト周波数(サンプリング周波数の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

类别

在 Help Center 和 File Exchange 中查找有关 Signal Processing Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!