14 views (last 30 days)
Show older comments
Sam on 20 Jun 2022
Commented: Star Strider on 21 Jun 2022
Accepted Answer: Star Strider
- example.csv
Open in MATLAB Online
Hi,
I have been looking over the documentation for findpeaks and several FWHM scripts, but I am not having any luck.
I have multiple .csv files of a single peak of Raman data, with the magnitude and wavelength. I would like to have a script that reads each .csv file and calculate the FWHM of the peak.
I currently have a script that reads my .csv files from a folder and spits out information such as peak postion, area etc., so maybe its easier to just tag it onto that, however it is using peakfit, not fitpeak. I have inserted below the section of the script which fits the peaks.
% Code for reading in files is above here!
% Now go through all the files in that folder and apply signal processing.
% Skip any files that generate and error and keep going.
for FileNum=1:NumFiles
try
datamatrix=load([DataDirectory '\' fn(FileNum,1:NameLength)]);
DataArray=[datamatrix(:,1) datamatrix(:,2)]; % Change to suit format of data
disp([num2str(FileNum) ': ' fn(FileNum,1:NameLength)])% Prints out the file number and name
disp(' Peak# Position Height Width Area')
% Curve fitting settings:
windowcenter=0; % Change this to fit other regions of the spectrum
windowwidth=0; % Change to zoom in or out to smaller or larger region.
NumPeaks=15; % Change this to the expected number of peaks.
Shape= 4; % 1-Gaussian, 2-Lorentzian, etc. Type "help peakfit" for list.
NumTrials=10; % Usually 1 to 10. Lower numbers faster but less stable.
startvector=0; % Difficult cases may require start vector, see help file.
BaselineMode=0; % Can be 0, 1, 2,or 3 (1=linear baseline subtraction)
[FitResults,GoodnessOfFit]=peakfit(DataArray,windowcenter,windowwidth,NumPeaks,Shape,1,NumTrials,startvector,BaselineMode);
disp(FitResults)
disp(' % fitting error R2')
disp(GoodnessOfFit)
drawnow
disp(' ')
% Add code here if you want to perform another set of
% operations on the same data file with different settings, i.e.,
% different data region, number of peaks, peak shape, etc.
catch me
disp([ num2str(FileNum) ': Error encounted with file: ' me.identifier])
disp(' ')
end
(I cannot remember who to credit with the above code, as I found this many months ago!)
Would it be best to just tag a FWHM code into this script, or write something completely separate?
Thanks in advance.
EDIT: Addition of example data.
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Sign in to answer this question.
Accepted Answer
Star Strider on 20 Jun 2022
Edited: Star Strider on 20 Jun 2022
Open in MATLAB Online
It might be easier if we had one or two typical .csv files to experiment with.
The findpeaks WidthReference option is 'halfprom' by default, although it can be set to 'halfheight' to calculate the FWHM value. If baseline offset or an inconstant baseline is a problem, there are several ways to deal with that.
EDIT — (20 Jun 2022 at 14:12)
Do you want to detrend the baseline first? If so, that is relatively straightforward —
D = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1039090/example.csv');
[~,ixu] = unique(D(:,1));
x = D(ixu,1);
y = D(ixu,2);
figure
plot(x, y)
grid
B = [x([1 end]) ones(2,1)] \ y([1 end]);
ydt = [x ones(size(y))] * B;
ybc = y-ydt; % Baselilne Corrected
figure
plot(x, ybc)
grid
It would be difficult to fit Gaussian peaks to this because the peak contributions are not well-defined. Calculating the FWHM of the single peak however would be relatively straightforward.
[pkmax,ixp] = max(ybc);
ixv = find(diff(sign(ybc-(pkmax/2))));
for k = 1:numel(ixv)
idxrng = max(1,ixv(k)-1) : min(numel(x),ixv(k)+1);
xv(k) = interp1(ybc(idxrng), x(idxrng), pkmax/2);
yv(k) = interp1(x(idxrng), ybc(idxrng), xv(k));
end
figure
plot(x, ybc)
hold on
plot(xv, yv, '+r')
hold off
text(x(ixp),pkmax/2, sprintf('\\leftarrow FWHM = %.1f \\rightarrow', diff(xv)), 'Horiz','center', 'Vert','middle', 'FontSize',7)
Make appropriate changes to get the result you want.
.
8 Comments Show 6 older commentsHide 6 older comments
Show 6 older commentsHide 6 older comments
Sam on 20 Jun 2022
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1744025-how-to-calculate-the-fwhm-from-multiple-spectra#comment_2224475
So sorry, I knew I was forgetting something. I have added the example data in now.
Is it possible to use both findpeaks and peakfit simultaneously or will this cause issues?
Sam on 20 Jun 2022
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1744025-how-to-calculate-the-fwhm-from-multiple-spectra#comment_2224600
Open in MATLAB Online
I have managed to get a fit from my data now using:
data = importdata('example.csv');
wavelengths = data(:, 1); % Wavelengths in column 1
magnitude = data(:, 2); % Signal magnitude in column 2
findpeaks(data(:,2),data(:,1),...
'MinPeakHeight',3,'MinPeakProminence',1,'Annotate','Extents'...
);
And from here, using the graph I can use the x values from the width line to calculate the FWHM
Sam on 20 Jun 2022
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1744025-how-to-calculate-the-fwhm-from-multiple-spectra#comment_2224780
That is great, thank you so much @Star Strider!
Looking back at my data, I think it would be better to fit a single peak from the whole spectrum, rather than "crop" the spectrum to the single peak and then fit it.
If I was to load in the whole spectrum, how would I "trim" it, so there is a specific window to select the correct peak from?
Star Strider on 20 Jun 2022
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1744025-how-to-calculate-the-fwhm-from-multiple-spectra#comment_2224875
Open in MATLAB Online
As always, my pleasure!
‘If I was to load in the whole spectrum, how would I "trim" it, so there is a specific window to select the correct peak from?’
That depends on the data. It would be necessary to select a specific peak manually or programmatically first. To get the ‘valleys’ on either side, use findpeaks (or islocalmin with find) to get the valley coordinate indices.
That would go something like this —
x = 920 : 980; % Create Data
y = sinc((x-mean(x))/5); % Create Data
figure
plot(x, y)
grid
[pks,plocs] = findpeaks(y, 'MinPeakProminence',1); % Desired Peak & Index
[vys,vlocs] = findpeaks(-y, 'MinPeakProminence',0.1); % Valleys & Indices
[vlc,ixv] = mink(abs(vlocs-plocs),2); % Choose Two Closest Valleys To Desired Peak
vlcs = plocs + [-1 1].*vlc; % Calculate Indices
idxr = vlcs(1) : vlcs(2); % Initial Index Range
figure
plot(x(idxr), y(idxr)) % Plot Desired Peak Region
grid
This would be a bit more difficult with a wandering baseline, so that would likely need to be addressed first, if it was a problem.
After that, my code simply analyses the peak it is presented with.
Sam on 20 Jun 2022
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1744025-how-to-calculate-the-fwhm-from-multiple-spectra#comment_2224920
This is great, thank you!
There seems to be an issue when I run this, the following error comes ups:
Matrix dimensions must agree.
Error in FWHMTrim (line 13)
[vlc,ixv] = mink(abs(vlocs - plocs),2); % Choose Two Closest Valleys To Desired Peak
Star Strider on 20 Jun 2022
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1744025-how-to-calculate-the-fwhm-from-multiple-spectra#comment_2225075
Open in MATLAB Online
As always, my pleasure!
I forgot that your data are column vectors. My data are row vectors, so transposing them to give a column vector requires:
vlcs = plocs + [-1; 1].*vlc; % Calculate Indices
for my code to work correctly. The mink (introduced in R2017b) call works correctly for the column vector.
This now appears to work correctly with column vector data (in R2022a) —
x = (920 : 980).'; % Create Data (Column)
y = sinc((x-mean(x))/5); % Create Data (Column)
figure
plot(x, y)
grid
[pks,plocs] = findpeaks(y, 'MinPeakProminence',1); % Desired Peak & Index
[vys,vlocs] = findpeaks(-y, 'MinPeakProminence',0.1); % Valleys & Indices
[vlc,ixv] = mink(abs(vlocs-plocs),2); % Choose Two Closest Valleys To Desired Peak
vlcs = plocs + [-1; 1].*vlc; % Calculate Indices
idxr = vlcs(1) : vlcs(2); % Initial Index Range
figure
plot(x(idxr), y(idxr)) % Plot Desired Peak Region
grid
My apologies for overlooking that detail.
.
Sam on 21 Jun 2022
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1744025-how-to-calculate-the-fwhm-from-multiple-spectra#comment_2225955
No need to apologise! Thank you so much for your help!
Star Strider on 21 Jun 2022
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1744025-how-to-calculate-the-fwhm-from-multiple-spectra#comment_2226120
As always, my pleasure!
Sign in to comment.
More Answers (0)
Sign in to answer this question.
See Also
Categories
Signal ProcessingSignal Processing ToolboxMeasurements and Feature ExtractionDescriptive Statistics
Find more on Descriptive Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- 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
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office
