How to calculate the FWHM from multiple spectra?

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.

Answer 1

[Star Strider]

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.

.

[Sam]

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]

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]

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]

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]

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]

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]

No need to apologise! Thank you so much for your help!

[Star Strider]

As always, my pleasure!