Batch FFT operation on EMG data

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.

Answer 1

[Star Strider]

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)

.

[Ines Shekhovtsov]

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]

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]

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]

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]

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.

See the documentation on medfreq for details.

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]

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]

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]

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]

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.

[Ines Shekhovtsov]

Thank you for all of your help!

[Star Strider]

As always, my pleasure!