Batch FFT operation on EMG data
8 次查看(过去 30 天)
显示 更早的评论
I have a basic code to read a .dat file containing EMG data for several muscles. The code scans the file for a specific time and extracts 3 seconds worth of data into a new table. I would like to calculate and plot the FFT spectrum for each EMG channel. The code i have works and gets the job done but i would like to optimize the code and not make it so rudamentary looking. I am very new to matlab coding and wrote the code with a lot of help from google. I have attached the code for reference. A couple of things i would like to streamline with new code:
*Perform the FFT analysis in a for loop or function
*Plot the FFT for each muscle in a for loop or function
*Will also welcome any other code tips for reading the .dat file and getting the table (T) in a more efficient way
Basically, as you will see from my code, it is very repetative and i assume i can get the same output by using a couple of for loops or functions instead of just writing this code in a brute force type method.
Thank you for your help in advance.
采纳的回答
Star Strider
2022-11-15
I am not certain what you want to do.
After you load the data and select the sections you want to analyse, convert the data in ‘Wanted’ to column vectors (if they are not column vectors already), then do something like this (assuming that they are now a matrix of column vectors):
Fs = ...; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
L = size(Wanted,1);
NFFT = 2^nextpow2(L); % For Efficiency
FTW = fft(Wanted, NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn; % Frequency Vector
Iv = 1:numel(Fv); % Index Vector
figure
plot(FV, abs(FTW(Iv,:))*2) % All Channels
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
NrSp = size(FTTW,2); % Number Of Subplots
for k = 1:NrSp
subplot(NrSp, 1, k)
plot(Fv, abs(FTW(Iv,k)*2)
xlabel(Frequency (Hz)')
ylabel('Magnitude')
title(sprintf('Column #%2d', k)
end
That may be a bit more efficient. I leave that to you to determine.
% C = websave('testCODE','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1194378/testCODE.m');
% type(C)
.
11 个评论
Ines Shekhovtsov
2022-11-15
Thank you. I am not sure how to convert the data in 'Wanted' to column vectors. But to clarify, i would like to work with the data in table T. That table starts with a time column, followed by 3 variables that i don't care about in my analysis and then the right and left muscle channels that i wish to perform the FFT analysis on. So basically column 5 to 32 are the billateral muscle channels i want to analyze. I am not even sure if i needed to create this table but if you see from my rudametary code, i was referring to each column (5-32) individually and calculating and plotting the FFT.
When i tried running your code after loading my file and extracting the 3 seconds i get the following error:
File: batchtest.m Line: 16 Column: 30
Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.
I changed the file name so for reference line 16 is:
plot(Fv, abs(FTW(Iv,k)*2)
Thank you!
Ines Shekhovtsov
2022-11-15
I am attaching the figure file that my code generates. I would like something like this to be the output of the new code. I am aware that the spacing is off but just wanted to give idea of what i want my output to be.
Thanks!
Star Strider
2022-11-15
iI did my best to proofread my posted code. I did not catch the missing parenthesis. It should be:
plot(Fv, abs(FTW(Iv,k))*2)
It would help significantly to have the ‘T’ table. Use the writetable function to write it to an .xlsx file or to a .txt file, then upload it here. I will look at the ‘test.fig’ file, however if you are working with the ‘T’ table, it would be best for me to have it to work with as well.
Ines Shekhovtsov
2022-11-16
编辑:Ines Shekhovtsov
2022-11-16
Thank you, I cannot believe that i didn't catch the parenthesis either! Attached is the file table.txt. It doesn't look right in notepad or excel but using readtable function will make it look right in matlab. As a reminder, i don't need to analyze the variables in columns 2-4. The billateral muscle names start with column 5 and end at column 36. And like i mentioned in previous comment, i would like to get a figure output from this code that is similar to my test.fig. I am aware that the margins are off. When i ran my code a week ago, the figure looked fine so i'm sure its an easy fix. And thank you for your help and time!
Star Strider
2022-11-16
As always, my pleasure!
Try this —
F = openfig(websave('test','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1194413/test.fig'));
sgtitle('Original Figure')
T = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1195293/table.txt')
T = 6002×36 table
times Tendonhammer CAP ManualPulse RSCM RIC6 RRA REO REST10 RPS RRF RVL RADD RTA RPL REDL RMH RFHB RMG RSOL LSCM LIC6 LRA LEO LEST10 LPS LRF LVL LADD LTA LPL LEDL LMH LFHB LMG LSOL
______ ____________ ______ ___________ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ ______ ______ ______ ______ ______ ______ ______ ______ ______ ______ ______ ______ ______ ______ ______ ______
30 -20.373 0.011 -0.0061 9.6136 -5.6079 -0.6409 -3.2045 -2.083 -4.1659 0.9614 -5.7682 1.442 -2.7239 -8.8125 -5.4477 -2.083 0 -8.492 1.6023 2.5687 5.7682 8.8125 7.05 21.631 13.139 8.6523 18.266 10.094 7.6909 10.415 7.05 9.1329 13.139 5.1273 7.3704
30 -20.382 0.0104 -0.0089 7.5307 -1.1216 0.4807 -0.6409 0 -3.8454 0.4807 -3.6852 -2.083 -2.2432 -4.967 -6.7295 -0.9614 -2.8841 -3.8454 2.5636 7.2527 11.216 4.1659 10.094 17.625 13.94 6.5693 20.669 7.05 6.7295 11.376 4.0057 11.216 7.3704 9.2932 5.4477
30 -20.382 0.0107 -0.0113 14.26 -6.5693 -7.3704 -12.017 -3.2045 -5.2875 -6.0886 -6.0886 -7.6909 -4.967 -8.1716 -4.967 -5.2875 3.2045 -3.525 1.1216 17.376 4.1659 5.4477 8.9727 20.189 8.8125 12.658 19.548 9.1329 5.4477 12.818 8.9727 11.697 10.895 8.8125 5.7682
30.001 -20.382 0.0101 -0.0128 12.658 -3.0443 -6.0886 -0.6409 2.2432 -5.7682 -0.9614 -4.3261 0.8011 -0.4807 -4.6466 -3.6852 -0.1602 -2.2432 -4.3261 -2.5636 28.255 9.7739 3.525 11.056 20.829 12.658 9.4534 17.945 11.056 7.5307 11.697 9.7739 11.536 12.978 11.857 8.9727
30.002 -20.373 0.0104 -0.0146 -5.6079 -5.2875 -7.3704 -9.9341 -2.083 0.6409 -3.6852 -1.1216 -1.2818 -1.1216 -1.9227 -4.1659 0.9614 -6.2489 -6.7295 -3.8454 33.997 9.2932 7.05 6.0886 18.426 15.222 11.376 24.995 12.658 9.7739 10.575 8.9727 14.26 17.465 9.4534 8.1716
30.002 -20.392 0.0104 -0.0156 -18.106 -8.9727 -3.2045 -10.575 -2.2432 1.1216 -1.9227 -1.6023 -0.4807 -5.4477 -0.4807 -3.3648 1.6023 -7.6909 -4.967 -4.967 30.673 10.735 11.857 14.1 15.542 15.061 13.78 23.714 10.094 10.735 12.818 12.017 12.818 18.426 14.26 12.978
30.003 -20.392 0.0107 -0.0143 -15.061 -2.5636 -5.4477 -5.7682 -5.4477 -0.3205 -5.4477 -1.2818 0.8011 -1.1216 -0.4807 -0.6409 5.1273 -8.1716 -7.6909 -3.2045 31.882 7.3704 11.056 16.183 18.747 13.78 9.9341 18.747 12.978 9.7739 14.741 12.658 12.177 15.222 10.575 11.857
30.003 -20.392 0.0107 -0.0146 3.525 -8.1716 -9.9341 -9.2932 -5.7682 -4.0057 -0.9614 1.7625 8.3318 2.083 0.9614 -2.5636 0.9614 -1.9227 -8.0114 -1.9227 26.895 10.094 12.818 17.144 22.111 16.824 9.6136 24.995 12.498 16.183 10.895 11.056 18.907 17.785 12.978 15.382
30.003 -20.401 0.0107 -0.0134 15.061 -4.1659 1.6023 -6.5693 2.083 -5.4477 -4.3261 -4.3261 2.4034 -4.4864 0 -3.2045 -0.1602 -5.7682 -8.6523 -0.6409 16.47 13.78 15.382 14.26 20.99 17.144 16.984 22.111 13.94 16.664 15.542 9.1329 16.503 25.476 12.498 16.503
30.004 -20.392 0.0104 -0.0116 -2.083 -3.3648 0 -7.6909 -2.2432 -5.7682 -1.9227 -5.7682 1.6023 -4.0057 -2.8841 -2.5636 -1.442 -8.9727 -9.6136 -2.7239 19.189 13.94 15.702 16.664 23.874 17.945 15.061 31.725 12.498 15.702 17.144 9.7739 13.94 21.15 13.619 13.78
30.005 -20.392 0.0107 -0.0085 -6.8898 -2.5636 1.442 -4.3261 0.8011 -4.6466 -2.2432 -4.0057 0.8011 -3.525 -2.5636 -3.8454 -0.6409 -6.2489 -9.2932 -4.0057 17.678 13.94 17.144 15.542 23.874 15.702 12.658 28.52 13.78 12.818 17.465 11.857 15.222 24.515 14.741 15.542
30.005 -20.401 0.0107 -0.0061 -5.9284 -0.6409 -0.9614 -8.8125 -1.6023 -9.6136 -5.7682 -3.8454 -1.442 -1.9227 -0.8011 -2.8841 -2.083 -4.6466 -9.6136 -0.1602 12.239 10.895 12.017 8.492 24.194 15.542 13.619 29.161 15.382 12.338 13.94 12.017 14.26 23.393 15.061 17.144
30.006 -20.382 0.0107 -0.0012 1.1216 -2.5636 -1.1216 -7.6909 -7.3704 -11.697 -2.5636 3.8454 0.3205 -4.1659 -4.967 -7.2102 -2.2432 -5.6079 -15.061 -3.6852 8.1593 8.492 17.625 10.575 20.189 19.067 12.177 31.725 13.139 12.658 14.42 8.8125 18.106 17.625 17.305 14.26
30.006 -20.392 0.011 0.0009 5.2875 -8.3318 0.9614 -4.3261 -0.8011 -14.741 -2.8841 -6.4091 -0.4807 -0.4807 -6.0886 -5.9284 0.6409 -10.575 -8.1716 -3.2045 14.203 9.9341 19.067 10.415 17.785 17.625 14.901 30.283 11.536 11.536 11.216 11.857 12.498 16.984 12.658 13.78
30.006 -20.382 0.0107 0.0052 -4.4864 -1.1216 2.5636 -9.7739 -1.1216 -14.901 -3.3648 -2.4034 -2.083 0.3205 -1.442 -6.7295 0 -6.0886 -6.7295 -2.7239 17.074 12.498 11.376 11.697 19.067 11.376 11.536 23.073 14.741 6.8898 11.216 9.6136 15.542 10.415 12.498 12.658
30.007 -20.401 0.0107 0.0095 5.6079 -7.2102 1.6023 -4.0057 2.2432 -15.061 -2.4034 -4.0057 1.442 -2.2432 -2.5636 -4.3261 0.3205 -10.895 -6.5693 1.1216 24.024 11.056 11.376 12.978 19.227 10.735 6.7295 27.078 13.459 11.536 13.94 9.4534 20.509 9.9341 11.536 11.376
Last32Cols = T(:,5:end);
VN = Last32Cols.Properties.VariableNames;
Wanted = table2array(Last32Cols)
Wanted = 6002×32
9.6136 -5.6079 -0.6409 -3.2045 -2.0830 -4.1659 0.9614 -5.7682 1.4420 -2.7239 -8.8125 -5.4477 -2.0830 0 -8.4920 1.6023 2.5687 5.7682 8.8125 7.0500 21.6307 13.1386 8.6523 18.2659 10.0943 7.6909 10.4148 7.0500 9.1329 13.1386
7.5307 -1.1216 0.4807 -0.6409 0 -3.8454 0.4807 -3.6852 -2.0830 -2.2432 -4.9670 -6.7295 -0.9614 -2.8841 -3.8454 2.5636 7.2527 11.2159 4.1659 10.0943 17.6250 13.9398 6.5693 20.6693 7.0500 6.7295 11.3761 4.0057 11.2159 7.3704
14.2602 -6.5693 -7.3704 -12.0170 -3.2045 -5.2875 -6.0886 -6.0886 -7.6909 -4.9670 -8.1716 -4.9670 -5.2875 3.2045 -3.5250 1.1216 17.3762 4.1659 5.4477 8.9727 20.1886 8.8125 12.6579 19.5477 9.1329 5.4477 12.8182 8.9727 11.6966 10.8954
12.6579 -3.0443 -6.0886 -0.6409 2.2432 -5.7682 -0.9614 -4.3261 0.8011 -0.4807 -4.6466 -3.6852 -0.1602 -2.2432 -4.3261 -2.5636 28.2552 9.7739 3.5250 11.0557 20.8295 12.6579 9.4534 17.9454 11.0557 7.5307 11.6966 9.7739 11.5363 12.9784
-5.6079 -5.2875 -7.3704 -9.9341 -2.0830 0.6409 -3.6852 -1.1216 -1.2818 -1.1216 -1.9227 -4.1659 0.9614 -6.2489 -6.7295 -3.8454 33.9969 9.2932 7.0500 6.0886 18.4261 15.2216 11.3761 24.9954 12.6579 9.7739 10.5750 8.9727 14.2602 17.4647
-18.1057 -8.9727 -3.2045 -10.5750 -2.2432 1.1216 -1.9227 -1.6023 -0.4807 -5.4477 -0.4807 -3.3648 1.6023 -7.6909 -4.9670 -4.9670 30.6728 10.7352 11.8568 14.1000 15.5420 15.0613 13.7795 23.7136 10.0943 10.7352 12.8182 12.0170 12.8182 18.4261
-15.0613 -2.5636 -5.4477 -5.7682 -5.4477 -0.3205 -5.4477 -1.2818 0.8011 -1.1216 -0.4807 -0.6409 5.1273 -8.1716 -7.6909 -3.2045 31.8816 7.3704 11.0557 16.1829 18.7466 13.7795 9.9341 18.7466 12.9784 9.7739 14.7409 12.6579 12.1773 15.2216
3.5250 -8.1716 -9.9341 -9.2932 -5.7682 -4.0057 -0.9614 1.7625 8.3318 2.0830 0.9614 -2.5636 0.9614 -1.9227 -8.0114 -1.9227 26.8954 10.0943 12.8182 17.1443 22.1113 16.8238 9.6136 24.9954 12.4977 16.1829 10.8954 11.0557 18.9068 17.7852
15.0613 -4.1659 1.6023 -6.5693 2.0830 -5.4477 -4.3261 -4.3261 2.4034 -4.4864 0 -3.2045 -0.1602 -5.7682 -8.6523 -0.6409 16.4696 13.7795 15.3818 14.2602 20.9897 17.1443 16.9841 22.1113 13.9398 16.6636 15.5420 9.1329 16.5034 25.4761
-2.0830 -3.3648 0 -7.6909 -2.2432 -5.7682 -1.9227 -5.7682 1.6023 -4.0057 -2.8841 -2.5636 -1.4420 -8.9727 -9.6136 -2.7239 19.1894 13.9398 15.7023 16.6636 23.8738 17.9454 15.0613 31.7250 12.4977 15.7023 17.1443 9.7739 13.9398 21.1500
Ts = 0.001; % Sampling Interval
Fs = 1/Ts; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
L = size(Wanted,1);
NFFT = 2^nextpow2(L); % For Efficiency
Wantedm = Wanted - mean(Wanted); % Subtract Mean To Show Other Peaks
FTW = fft(Wantedm, NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn; % Frequency Vector
Iv = 1:numel(Fv); % Index Vector
figure
plot(Fv, abs(FTW(Iv,:))*2) % All Channels
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
figure
NrSp = size(FTW,2); % Number Of Subplots
Cols = 2; % Number Of Subplot Columns
Rows = fix(NrSp/Cols);
RL = [2:2:NrSp 1:2:NrSp]; % 'Right-Left' Vector
for k = 1:NrSp
subplot(Rows, Cols, RL(k))
plot(Fv, abs(FTW(Iv,k))*2)
[mf(k,:),bp(k,:)] = medfreq(Wantedm(:,k),Fs); % Median Frequency & Bandpower For Each Record
% xlabel('Frequency (Hz)')
% ylabel('Magnitude')
% title(sprintf('Column #%2d', k))
xlim([min(Fv) max(Fv)])
title(VN{k})
% [mv,ix] = max(abs(FTW(Iv,k))); % Spectrum Maxima, Indices, & Associated Frequencies
% Fvix(k,:) = Fv(ix);
end
sgtitle('This Study')
% Fvix
Tfp = table(mf,bp, 'VariableNames',{'MedianFreq','BandPower'}, 'RowNames',VN)
Tfp = 32×2 table
MedianFreq BandPower
__________ _________
RSCM 69.333 247.05
RIC6 24.493 13.644
RRA 157.6 7.5139
REO 29.529 29.824
REST10 80.809 12.586
RPS 60.38 12.457
RRF 110.43 6.7567
RVL 160.02 6.3267
RADD 110.66 8.0695
RTA 36.699 10.957
RPL 126.52 7.909
REDL 234.73 4.1582
RMH 111.76 7.6997
RFHB 60.851 10.009
RMG 119.67 7.5247
RSOL 106.89 8.0257
In the fft call, I first subtracted the mean (D-C offset) from the data before transforming it to the frequency domain in order to show the other peaks. This was apparently not done in the original data, since there is a peak at 0 Hz (the D-C offset) and the other peaks are obscured as the result.
The ‘RL’ vector assigns the ‘right’ subplots to the second column and the corresponding ‘left’ subplots to the first column. There might be other ways to do that, however in my experience, this way is simply easier.
I added:
[mf(k),bp(k)] = medfreq(Wantedm(:,k),Fs); % Median Frequency & Bandpower For Each Record
out of interest, since thse data may be important in such studies, and a table of the results. The median frequency is in Hz since that is the units of ‘Fs’, and the units of bandpower would be whatever the amplitudes units are, squared, so if in millivolts, mV^2.
I also noticed that there appears to be a relatively prominent 30 Hz spike in many of the spectra. I have no idea where that originates, however it could be worth using a notch filter in the time domain signals to get rid of it. One way to do that would be just after ‘Wanted’ is first created to add:
WantedFilt = bandstop(Wanted, [29 31], Fs, 'ImpulseResponse','iir');
Then use ‘WantedFilt’ in place of ‘Wanted’ in the ‘Wantedm’ calculation. The rest of the code is unchanged.
I did my best to proofread this this time. I don’t see any errors.
.
Ines Shekhovtsov
2022-11-16
编辑:Ines Shekhovtsov
2022-11-16
WOW, this is amazing. Thank you so much!!! This code will be a helpful learning tool for me when i need to write other emg analysis codes in the future. A quick side question, is the BandPower similar to Total Power analysis of a signal?? I am also interested in total power analysis and have been using code written by a previous programmer. I've looked up total power analysis and also sometimes see the term power spectrum analysis. I don't fully understand the definition of total power analysis, so any brief clarification in that area will be appreciated! Thank you.
Star Strider
2022-11-16
As always, my pleasure!
I am not certain what ‘total power‘ is or refers to. I might have encountered it by a different name, so if you have a reference to it, I would appreciate seeing it. If it is a PDF file, please attach it here, since if it is behind a paywall, I would likely not have access to it. Looking at the bandpower documentation, the two could be related, with bandpower being the mean of the power of the signal, and total power being the integral of the power in the signal.
Perhaps it would simply be something like taking the integral of the square of a signal, so:
TotlPwr = trapz(t, s.^2);
BandPwr = trapz(t, s.^2)/numel(t);
where ‘t’ is the time vector and ‘s’ is the time domain signal.
This is just a guess.
Ines Shekhovtsov
2022-11-16
We are just starting to use this analysis and i need to get better clarification from our researchers on what they actually want calculated. However i found this snipet of code from a programmer who is no longer here.
Total_Power(j,1) = double(sum(abs(data_int).^2)./numel(data_int));
It looks similar to what you have shown, however i am not familiar with the trapz function.
Thank you, as always!
Star Strider
2022-11-16
As always, my pleasure!
The trapz version, doing trapezoidal integration, is much more robust (and likely more accurate) than using sum. You should have trapz. Also, it is inefficient to calculate the absolute value first, then squaring the result. Just square the vector element-wise within the trapz call, as I did.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)