selection of specific data from graph and extract data

%For given graph only peck values between 10 and 20hz is valuble and i like to extract phase value for specific point. 
%Right now my peck function is taking all the small changes, is there any way i can pass some filter or anything?
%Code:
header = 9;
delimiter = '\t';
ind = 1:100;
for i = 1:4
    filname = sprintf('H1, 2_I1sv%05d.txt',i);
    dat(i) =importdata(filname,delimiter,header);
    FRFdata = dat(i).data;
    f = FRFdata(:,1);
    a = FRFdata(:,2);
    b = FRFdata(:,3);
    amp = sqrt(a.^2 + b.^2);
    pha = (180/pi)* atan2(b,a);
    [Ypk,Xpk,Wpk1,Ppk] = findpeaks(amp(ind));
    figure(i)
    subplot(2,1,1),plot(f,amp,f(Xpk),Ypk,'dr') 
    grid on
    ylabel('Amp')
    xlim([0 40])
    ylim([0 2.4]) 
    subplot(2,1,2),plot(f,pha)
    grid on
    xlim([0 40])
    xlabel('Hz')
    ylabel('Phase')
end
%I am doing FRF analysis to find vibration modes

Answer 1

[Star Strider]

The findpeaks function has a number of name-value pair arguments that can be used to return only selected information.

See if:

[Ypk,Xpk,Wpk1,Ppk] = findpeaks(amp(ind), 'MinPeakProminence',0.5);

does what you want.

I also tweaked your code a bit to make it abit mor efficient, specifically using readmatrix and changing ‘amp’ and ‘pha’

header = 9;
delimiter = '\t';
ind = 1:100;
i = 1;
% for i = 1:4
%     filname = sprintf('H1, 2_I1sv%05d.txt',i);
    filname = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323755/H1,%202_I1sv00001.txt';
%     dat(i) =importdata(filname,delimiter,header);
    FRFdata = readmatrix(filname, 'HeaderLines',9);
%     FRFdata = dat;
    f = FRFdata(:,1);
    a = FRFdata(:,2);
    b = FRFdata(:,3);
    amp = hypot(a,b);                                                       % Changed
    pha = atan2d(b,a);                                                      % Changed
    [Ypk,Xpk,Wpk1,Ppk] = findpeaks(amp(ind), 'MinPeakProminence',0.5);      % Changed
    figure(i)
    subplot(2,1,1),plot(f,amp,f(Xpk),Ypk,'dr') 
    grid on
    ylabel('Amp')
    xlim([0 40])
    ylim([0 2.4]) 
    subplot(2,1,2),plot(f,pha, f(Xpk),pha(Xpk),'dr')                        % Changed
    grid on
    xlim([0 40])
    xlabel('Hz')
    ylabel('Phase')

% end

.

[AL]

This code is perfectly highlighting required points; however, variable Ppk, Wpk1, Xpk, Ypk all are emply do you know why this is happening? I wanted to record peak values for Amp and respective Phase angles for respective frequency.

header = 9;

delimiter = '\t';

ind = 1:100;

for i = 1:4

filname = sprintf('H1, 2_I1sv%05d.txt',i);

dat(i) =importdata(filname,delimiter,header);

FRFdata = dat(i).data;

f = FRFdata(:,1);

a = FRFdata(:,2);

b = FRFdata(:,3);

amp = hypot(a,b);

pha = atan2d(b,a);

[Ypk,Xpk,Wpk1,Ppk] = findpeaks(amp(ind), 'MinPeakProminence',0.05); %Changed

figure(i)

subplot(2,1,1),plot(f,amp,f(Xpk),Ypk,'dr')

grid on

ylabel('Amp')

xlim([0 40])

ylim([0 2.4])

subplot(2,1,2),plot(f,pha, f(Xpk),pha(Xpk),'dr')

grid on

xlim([0 40])

xlabel('Hz')

ylabel('Phase')

end

clear a amp b f filname i header delimiter dat ind %added

[Star Strider]

The fourth data file required a different approach, the problem being the ‘ind’ vector.

I adapted the code to work with all of them, and also to save the findpeaks output to individual cell arrays —

filenamesc = {'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323755/H1,%202_I1sv00001.txt'; 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323760/H1,%202_I1sv00002.txt'; 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323765/H1,%202_I1sv00003.txt'; 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323770/H1,%202_I1sv00004.txt'};
header = 9;
delimiter = '\t';
ind = 1:100;
i = 1;
for i = 1:numel(filenamesc)
%     filname = sprintf('H1, 2_I1sv%05d.txt',i);
    filname = filenamesc{i};
%     dat(i) =importdata(filname,delimiter,header);
    FRFdata = readmatrix(filname, 'HeaderLines',9);
    ind = 20:size(FRFdata,1);                                               % Change 'ind'
%     FRFdata = dat;
    f = FRFdata(:,1);
    a = FRFdata(:,2);
    b = FRFdata(:,3);
    amp = hypot(a,b);                                                       % Changed
    pha = atan2d(b,a);                                                      % Changed
    mpp = max(amp(ind))/4;                                                  % Adaptive 'MinPeakProminence' Value 
    [Ypk,Xpk,Wpk,Ppk] = findpeaks(amp(ind), 'MinPeakProminence',mpp, 'WidthReference','halfheight');       % Changed
    Xpk = Xpk+min(ind)-1;
    Ypkc{i} = Ypk;                                                          % Save To Cell Array
    Xpkc{i} = Xpk;                                                          % Save To Cell Array 
    Wpkc{i} = Wpk;                                                          % Save To Cell Array   
    Ppkc{i} = Ppk;                                                          % Save To Cell Array   
    
    figure(i)
    subplot(2,1,1),plot(f,amp,f(Xpk),Ypk,'dr') 
    grid on
    ylabel('Amp')
    xlim([0 40])
    ylim([0 2.4]) 
    subplot(2,1,2),plot(f,pha, f(Xpk),pha(Xpk),'dr')                        % Changed
    grid on
    xlim([0 40])
    xlabel('Hz')
    ylabel('Phase')
    sgtitle("File "+string(i))
end

Results = table([Ypkc{:}],[Xpkc{:}],[Wpkc{:}],[Ppkc{:}], 'VariableNames',{'Ypk','Xpk','Wpk','Ppk'})

Results = 2×4 table
                        Ypk                                  Xpk                              Wpk                                        Ppk                    
    ___________________________________________    _______________________    ____________________________________    __________________________________________

     1.5942      1.5843      1.5803      0.1259    39     39     39    135    1.6724    1.7203    1.6891     3.465    1.5777       1.343      1.3789     0.11466
    0.60418      2.2867      1.9136    0.045885    52     52     52    191    1.4212    1.4796    1.4388    7.6751    0.5718      2.2838      1.8899    0.032798

.

[AL]

I really appricate your effort one last question can you please explain me the logic behind following lines:

ind = 20:size(FRFdata,1); %why start from 20?

mpp = max(amp(ind))/4;

Like I tried playing with them to change the threshold value which can inclue the values for those frequancy as well shown in figure below. whenever, I change values of those it creates an error saying

Error using horzcat

Dimensions of arrays being concatenated are not consistent. (Error will be in Results)

this is my latest code:

header = 9;

delimiter = '\t';

ind = 1:100;

for i = 1:4

filname = sprintf('H1, 2_I1sv%05d.txt',i);

dat(i) =importdata(filname,delimiter,header);

FRFdata = readmatrix(filname, 'HeaderLines',9);

ind =20:size(FRFdata,1);

f = FRFdata(:,1);

a = FRFdata(:,2);

b = FRFdata(:,3);

amp = hypot(a,b);

pha = atan2d(b,a);

mpp = max(amp(ind))/4; % Adaptive 'MinPeakProminence' Value

[Ypk,Xpk,Wpk,Ppk] = findpeaks(amp(ind), 'MinPeakProminence',mpp, 'WidthReference','halfheight'); % Changed

Xpk = Xpk+min(ind)-1;

Ypkc{i} = Ypk; % Save To Cell Array

Xpkc{i} = Xpk; % Save To Cell Array

Wpkc{i} = Wpk; % Save To Cell Array

Ppkc{i} = Ppk; % Save To Cell Array

figure(i)

subplot(2,1,1),plot(f,amp,f(Xpk),Ypk,'dr')

grid on

ylabel('Amp')

xlim([0 40])

ylim([0 2.4])

subplot(2,1,2),plot(f,pha, f(Xpk),pha(Xpk),'dr') % Changed

grid on

xlim([0 40])

xlabel('Hz')

ylabel('Phase')

sgtitle("File "+string(i))

end

Results = table([Ypkc{:}],[Xpkc{:}],[Wpkc{:}],[Ppkc{:}], 'VariableNames',{'Ypk','Xpk','Wpk','Ppk'})

also, I am getting results for file 2,3, and 4 no graphs or results for data1.

any idea why this is happening?

I have tried:

ind =10:size(FRFdata,1); %ind =20:size(FRFdata,1);

mpp = max(amp(ind))/22; %mpp = max(amp(ind))/4;

I like to make this code generalized to find peak points and collect data for different dataset as well.

If you please explain me the logic it whould be helpful.

[Star Strider]

As always, my pleasure!

The choice of indices for ‘ind’ is arbitrary, and appears to work satisfactorily. Beginning at 20 eliminates the low-frequency peak in the beginning, because including it would cause problems for the ‘mpp’ result.

Those peaks other do not appear in all the spectra. Change the ‘mpp’ result to detect all the peaks you want to return.

The ‘Results’ table is simply for convenience. It is not required, and can safely be deleted. The data are saved in the various cell arrays, so nothing is lost by deleting ‘Results’.

This sorts all the returned data by ‘Ppk’ and returns the 4 peaks with the highest prominence, regardless of the 'MinPeakProminence' value, because it is not used here. (This also make the ‘Results’ table work again.)

You could also add ‘Fpk{i} = f(Xpk{i})’ to return the frequencies if you want to. If you add it to ‘Results’ be sure to add a 'VariableNames' entry for it as well.

The latest version —

filenamesc = {'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323755/H1,%202_I1sv00001.txt'; 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323760/H1,%202_I1sv00002.txt'; 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323765/H1,%202_I1sv00003.txt'; 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1323770/H1,%202_I1sv00004.txt'};
header = 9;
delimiter = '\t';
ind = 1:100;
i = 1;
for i = 1:numel(filenamesc)
%     filname = sprintf('H1, 2_I1sv%05d.txt',i);
    filname = filenamesc{i};
%     dat(i) =importdata(filname,delimiter,header);
    FRFdata = readmatrix(filname, 'HeaderLines',9);
    ind = 2:size(FRFdata,1);                                               % Change 'ind'
%     FRFdata = dat;
    f = FRFdata(:,1);
    a = FRFdata(:,2);
    b = FRFdata(:,3);
    amp = hypot(a,b);                                                       % Changed
    pha = atan2d(b,a);                                                      % Changed
    mpp = max(amp(ind))/30;                                                 % Adaptive 'MinPeakProminence' Value 
%     [Ypk,Xpk,Wpk,Ppk] = findpeaks(amp(ind), 'MinPeakProminence',mpp, 'WidthReference','halfheight');       % Changed
    [Ypk,Xpk,Wpk,Ppk] = findpeaks(amp(ind), 'WidthReference','halfheight');       % Changed
    Xpk = Xpk+min(ind)-1;
    [pv,mxidx] = maxk(Ppk,4);
    Ypkc{i} = Ypk(mxidx);                                                          % Save To Cell Array
    Xpkc{i} = Xpk(mxidx);                                                          % Save To Cell Array 
    Wpkc{i} = Wpk(mxidx);                                                          % Save To Cell Array   
    Ppkc{i} = Ppk(mxidx);                                                          % Save To Cell Array   
    
    figure(i)
    subplot(2,1,1),plot(f,amp,f(Xpkc{i}),Ypkc{i},'dr') 
    grid on
    ylabel('Amp')
    xlim([0 40])
    ylim([0 2.4]) 
    subplot(2,1,2),plot(f,pha, f(Xpkc{i}),pha(Xpkc{i}),'dr')                        % Changed
    grid on
    xlim([0 40])
    xlabel('Hz')
    ylabel('Phase')
    sgtitle("File "+string(i))
end

Results = table([Ypkc{:}],[Xpkc{:}],[Wpkc{:}],[Ppkc{:}], 'VariableNames',{'Ypk','Xpk','Wpk','Ppk'})

Results = 4×4 table
                        Ypk                                   Xpk                              Wpk                                           Ppk                      
    ____________________________________________    _______________________    ____________________________________    _______________________________________________

      1.5942      2.2867      1.9136      0.1259    39     52     52    135    1.6724    1.4796    1.4388     3.465      1.5777       2.2838        1.912      0.11466
     0.60418      1.5843      1.5803    0.045885    52     39     39    191    1.4212    1.7203    1.6891    4.7291      0.5718        1.343       1.3789     0.032798
    0.090036     0.29818     0.45969    0.018411    70     69     69     62         4    2.4146    2.0969    2.8026     0.03277      0.29515      0.42127      0.01134
    0.035399     0.10936     0.15162    0.015571     4     17     17     34    1.6269    1.2849    1.2563         3    0.026676      0.10265      0.12791    0.0074241

.

[AL]

Good morning,

Thank you so much for solving my doubts.

This is my latest draft of the code:

header = 9;

delimiter = '\t';

num_files = 4;

Ypkc = cell(num_files, 1);

Xpkc = cell(num_files, 1);

Wpkc = cell(num_files, 1);

Ppkc = cell(num_files, 1);

freq = cell(num_files, 1);

pha_ = cell(num_files, 1);

for i = 1:num_files

filname = fullfile('.', sprintf('H1, 2_I1sv%05d.txt',i));

if exist(filname, 'file') ~= 2

warning('File %s does not exist', filname);

continue;

end

try

dat = importdata(filname, delimiter, header);

FRFdata = readmatrix(filname, 'HeaderLines', header);

f = FRFdata(:,1);

a = FRFdata(:,2);

b = FRFdata(:,3);

amp = hypot(a,b);

pha = atan2d(b,a);

ind = 2:size(FRFdata,1);

mpp = max(amp(ind))/30;

[Ypk,Xpk,Wpk,Ppk] = findpeaks(amp(ind),'MinPeakProminence',mpp, 'WidthReference','halfheight');

Xpk = Xpk + min(ind) - 1;

[pv,mxidx] = maxk(Ppk,4); % to get top 4 peak prominences

Ypkc{i} = Ypk(mxidx); % peak amp

Xpkc{i} = Xpk(mxidx); % peak locations

Wpkc{i} = Wpk(mxidx); % peak widths

Ppkc{i} = Ppk(mxidx); % peak prominences

freq{i} = f(Xpkc{i}); % Frequency

pha_{i} = pha(Xpkc{i}); % Phase

fprintf('File %d processed successfully\n', i);

catch ME

warning('Error processing file %s: %s', filname, ME.message);

end

figure(i)

subplot(2,1,1), plot(f, amp, f(Xpk), Ypk, 'dr')

grid on

ylabel('Amp')

xlim([0 40])

ylim("auto")

subplot(2,1,2), plot(f, pha, f(Xpk), pha(Xpk), 'dr')

grid on

xlim([0 40])

xlabel('Hz')

ylabel('Phase')

sgtitle("File " + string(i))

end

Results = table([freq{:}],[Ypkc{:}],[pha_{:}],[Xpkc{:}], 'VariableNames',{'Frequency','peak amp','Phase','peak locations '})

Result = sortrows(Results, 'Frequency','ascend');

clear a amp b dat delimiter f filname freq FRFdata;

clear header i ind mpp mxidx Ppk Ppkc pv Wpk Wpkc Xpk Xpkc Ypk Ypkc pha pha_ num_files;

I am only facing one problem related to the results.

I tried sorting results with respect to the frequency; however, the code which i have wrote will sort 1st column of the frequncy in ascending order and change all the other raws with respect to that.

do you know how can i make a table which will look something like this?

[AL]

Dear Star Strider,

I have solved the issue.

% Combine the results from all files into a single table

all_results = table();

for i = 1:num_files

% Add the results for the current file to the table

curr_results = table(freq{i}, Ypkc{i}, pha_{i}, Xpkc{i}, ...

'VariableNames', {'Frequency', 'PeakAmplitude', 'Phase', 'PeakLocation'});

curr_results = sortrows(curr_results, 'Frequency','ascend');

% Add a suffix to the variable names to avoid duplicates

suffix = ['_file' num2str(i)];

curr_results.Properties.VariableNames = strcat(curr_results.Properties.VariableNames, suffix);

all_results = [all_results curr_results]; % Concatenate horizontally

end

% Display the final table

disp(all_results);

clear a amp b dat delimiter f filname freq FRFdata suffix curr_results Results;

clear header i ind mpp mxidx Ppk Ppkc pv Wpk Wpkc Xpk Xpkc Ypk Ypkc pha pha_ num_files;

Have a wonderful day. Thank you so much for your effort and it was a good experience learning new things with you.

[Star Strider]

As always, my pleasure!

It appears that you got everything working the way you want!

You can put the file name in the sgtitle if you want to:

sgtitle("File: H1, 2\_I1s" + i + "v.txt")

.The backslant (\) is necessary to ‘escape’ the underscore so that it prints as an underscore and does not subscript the ‘I’ that follows it.

.