How to run batch of peak finding on data?

I am writing a script in order to locate the exact position of one peak, run this on hundred + sample .dat files and then spit out a .txt file with the peak position for each sample/.dat file. Unfortunately, I am struggling to use fitpeaks.

The data looks like this and I am looking to find the peak position at approx. 0.01.

data = readmatrix('exampledata.dat') ;   % load the data from text file
graph = plot(data(:,1),data(:,2)) 

findpeaks(graph)

It plots the data nicely, but won't find any peaks, so I get this error:

Error using findpeaks
Expected Y to be one of these types:

double, single

Instead its type was matlab.graphics.chart.primitive.Line.
Error in findpeaks>parse_inputs (line 199)
validateattributes(Yin,{'double','single'},{'nonempty','real','vector'},...
Error in findpeaks (line 136)
  = parse_inputs(isInMATLAB,Yin,varargin{:});
Error in PeakFinderS (line 4)
findpeaks(graph) 

I have attached example data.

Answer 1

[Star Strider]

data = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/895415/exampledata.txt')

data = 1268×3

0.0017 0.0103 0.0001 0.0018 0.0730 0.0004 0.0019 0.2810 0.0006 0.0021 0.2447 0.0004 0.0022 0.2747 0.0003 0.0023 0.3600 0.0003 0.0024 0.3246 0.0002 0.0025 0.2888 0.0002 0.0026 0.2204 0.0002 0.0028 0.2311 0.0002

[pks,locs,wdt,prm] = findpeaks(data(:,2), 'MinPeakDistance', 50);

figure
graph = plot(data(:,1),data(:,2));
hold on
plot(data(locs,1), data(locs,2), '^r')
hold off
xlim([0 0.05])
text(data(locs(2),1), data(locs(2),2), sprintf('  \\leftarrow (%.4f, %.4f)',data(locs(2),1), data(locs(2),2)), 'Horiz','left', 'Vert','bottom', 'Rotation',45)

I am not certain what the desired result is. This is one way of getting something close to it.

.

[Sam]

Thank you @Star Strider, always a huge help!

Is there a way to run this on several .dat files and then produce a text file of the resulting peak positions?

[Star Strider]

My pleasure!

See Import Or Export a Sequence Of Files for that part. However, it would be necessary to know if the peak of interest is always at about the same position or if it always the same peak or peak amplitude, or other details about it in order to automate it. Beyond that, it is simply a matter of reading each file (and perhaps plotting its data), identifying the peak, storing the peak coordinates in a matrix, and then saving the matrix to a file at the end.

With only one file and no other information about the characteristics of the desired peak, I cannot be more specific.

[Sam]

Yes, it will always be the same peak. Its position will range from 0.0085 to 0.01.

I have 35,000 .dat files so will definitely need to automate it somehow, aha! Thank you, I will take a look at the link now.

[Star Strider]

My pleasure!

If it is always going to be in the same region, that makes it easier —

data = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/895415/exampledata.txt')

data = 1268×3

0.0017 0.0103 0.0001 0.0018 0.0730 0.0004 0.0019 0.2810 0.0006 0.0021 0.2447 0.0004 0.0022 0.2747 0.0003 0.0023 0.3600 0.0003 0.0024 0.3246 0.0002 0.0025 0.2888 0.0002 0.0026 0.2204 0.0002 0.0028 0.2311 0.0002

[pks,locs,wdt,prm] = findpeaks(data(:,2), 'MinPeakDistance', 50);

Iv = find((data(:,1) >= 0.007) & (data(:,1) <= 0.015));             % Indices Defining Region-Of-Interest
Peak_01 = Iv(ismember(Iv,locs))                                     % The 'locs' Withiin The Region-Of-Interest

Peak_01 = 70

figure
graph = plot(data(:,1),data(:,2));
hold on
plot(data(Peak_01,1), data(Peak_01,2), '^r')
hold off
xlim([0 0.05])
text(data(Peak_01,1), data(Peak_01,2), sprintf('  \\leftarrow (%.4f, %.4f)',data(Peak_01,1), data(Peak_01,2)), 'Horiz','left', 'Vert','bottom', 'Rotation',45)

The changed (updated) code now selects the region-of-interest (the ‘Iv’ vector of indices) and then returns the peak location within the region-of-interest first as a logical vector with the ismember call and then the numeric index as ‘Peak_01’.

The code now finds and plots the desired peak location and amplitude. I have illustrated how to retrieve those, so it will be straightforward to write them to a matrix and then to a file when the file processing loop completes. Processing the files should be straightforward with the procedures in the link in my previous Comment.

.

[Sam]

Thanks again @Star Strider.

I am a unclear as to where the value for peak position appears in the workspace (if it is supposed to)? When I look under locs, the value is not there.

I have my script for exporting the file:

writematrix(locs, 'peak_positions.txt', 'Delimiter', 'tab')
type 'peak_positions.txt'

, but it appears that locs is not what I need to export.

[Star Strider]

The ‘Peak_01’ value is the index into ‘data’ so the position and magnitude are ‘PeakMtx’ in this example code —

for k = 1:numfiles

data = readmatrix(filename{k});

[pks,locs,wdt,prm] = findpeaks(data(:,2), 'MinPeakDistance', 50);

Iv = find((data(:,1) >= 0.007) & (data(:,1) <= 0.015)); % Indices Defining Region-Of-Interest

Peak_01 = Iv(ismember(Iv,locs)); % The 'locs' Withiin The Region-Of-Interest

PeakMtx(k,:) = [data(Peak_01,1), data(Peak_01,2)];

end

writematrix(PeakMtx, 'peak_positions.txt', 'Delimiter', 'tab')

type 'peak_positions.txt'

The ‘PeakMtx’ matrix has the independent variable (location) value in the first column, and the corresponding peak amplitude in the second column. (Add as many other variables as desired, for example ‘Peak_01’, however they all must be numeric in this code.)

The ‘numfiles’ vector are the number of files in the directory you are accessing. Your code needs to return a vector of file names, preferably as a cell array (since that is what my code assumes), and number of files in that vector is the ‘numfiles’ value. (You can use the numel function to return that value.) Each file name is referenced as ‘filename{k}’ in the loop, and may also need to include the full path (directory) information as well as part of that argument, using the fullfile function. That is as specific as I can get, so please refer to the documentation I referred to earlier with respect to processing the sequence of files.

.

[Sam]

Thank you @Star Strider, you've been very helpful!

[Star Strider]

As always, my pleasure!

[Sam]

Sorry to message again, but I've spent all day trying to get this to work, but I have had no luck. I have tried several ways to produce the .txt file with all the relevant data.

This is my current code:

% Specify the folder where the files live.
myFolder = 'E:\sm27583-2\Processed SAXS\DSpacing\i22-569969_dat\i22-569969_ascii';
% Check to make sure that folder actually exists.  Warn user if it doesn't.
if ~isfolder(myFolder)
    errorMessage = sprintf('Error: The following folder does not exist:\n%s\nPlease specify a new folder.', myFolder);
    uiwait(warndlg(errorMessage));
    myFolder = uigetdir(); % Ask for a new one.
    if myFolder == 0
         % User clicked Cancel
         return;
    end
end
% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(myFolder, '*.dat'); % Change to whatever pattern you need.
theFiles = dir(filePattern);
for k = 1 : length(theFiles)
    baseFileName = theFiles(k).name;
    fullFileName = fullfile(theFiles(k).folder, baseFileName);
    fprintf(1, 'Now reading %s\n', fullFileName);
    % Now do whatever you want with this file name,
    data = importdata(fullFileName)

[pks,locs,wdt,prm] = findpeaks(data(:,2), 'MinPeakDistance', 50);
Iv = find((data(:,1) >= 0.007) & (data(:,1) <= 0.015)); % Indices Defining Region-Of-Interest
Peak_01 = Iv(ismember(Iv,locs)); % The 'locs' Withiin The Region-Of-Interest
PeakMtx(k,:) = [data(Peak_01,1), data(Peak_01,2)];

end

writematrix(PeakMtx, 'peak_positions.txt', 'Delimiter', 'tab')
type 'peak_positions.txt'.

Unfortunately, this does not work and I get this error:

Unable to perform assignment because the size of the left side is 1-by-2 and the size of the right side is

0-by-2.

Error in PeakFinderS (line 26)

PeakMtx(k,:) = [data(Peak_01,1), data(Peak_01,2)];

I am unsure why there now appears to be an issue reading the .dat files in, this was not a problem previously.

[Star Strider]

That problem is likely caused by ‘Peak_01’ being empty. That will occur if the ismember call returns an empty vector (meaning that there are no ‘loc’ indices in the ‘Iv’ vector).

The easiest way to deal with that is to check for that and branch around it, however another option is to assign NaN values to ‘PeakMtx’ in case the order is important and the peak information are to be matched to the files —

[pks,locs,wdt,prm] = findpeaks(data(:,2), 'MinPeakDistance', 50);

Iv = find((data(:,1) >= 0.007) & (data(:,1) <= 0.015)); % Indices Defining Region-Of-Interest

Peak_01 = Iv(ismember(Iv,locs)); % The 'locs' Withiin The Region-Of-Interest

if ~isempty(Peak_01)

PeakMtx(k,:) = [data(Peak_01,1), data(Peak_01,2)];

end

alternatively —

[pks,locs,wdt,prm] = findpeaks(data(:,2), 'MinPeakDistance', 50);

Iv = find((data(:,1) >= 0.007) & (data(:,1) <= 0.015)); % Indices Defining Region-Of-Interest

Peak_01 = Iv(ismember(Iv,locs)); % The 'locs' Withiin The Region-Of-Interest

if isempty(Peak_01)

PeakMtx(k,:) = NaN(1,2);

else

PeakMtx(k,:) = [data(Peak_01,1), data(Peak_01,2)];

end

That only changes my previous code. It does not change the loop or the rest of the code.

(I was in the process of writing this when my desktop computer chose to update some of its software. That took a few minutes. My apologies for the delay.)

.

[Sam]

No need to apolgise, you respond very quickly, its really appreciated!

Does this mean that if there are samples/.dat files where there is no peak present, it will break the code?

Unfortunately for both new codes I get this error instead:

Unable to perform assignment because the size of the left side is 1-by-2 and the size of the right side is

2-by-2.

[Star Strider]

‘Does this mean that if there are samples/.dat files where there is no peak present, it will break the code?’

If no provisions are made for that condition, yes. It is not possible to fill a numeric vector with empty values.

A fix for the previous error (empty ‘Peak_01’) and the new error is to save it to a cell array instead of a numeric array. The loop will then complete, however it will be necessary to sort out the resulting cell array afterwards.

That code would be —

[pks,locs,wdt,prm] = findpeaks(data(:,2), 'MinPeakDistance', 50);

Iv = find((data(:,1) >= 0.007) & (data(:,1) <= 0.015)); % Indices Defining Region-Of-Interest

Peak_01 = Iv(ismember(Iv,locs)); % The 'locs' Withiin The Region-Of-Interest

PeakMtx{k,:} = [data(Peak_01,1), data(Peak_01,2)];

The only change is substituting curly braces {} for parentheses () in the ‘PeakMtx’ assignment, denoting cell array indexing. It might still be best to fill the empty array with NaN values, since you can then use the cell2mat function to create a numeric matrix from the cell array. Experiment with and without it to see what works best in your application.

.

[Sam]

This has worked! Thank you so much @Star Strider, you have been such a huge help, I really appreciate your persistance to help me!

[Star Strider]

As always, my pleasure!

This is an interesting problem. I never considered using the technique here to isolate peaks in a specific region of the independent variable before,so I learned something.