How to calculate PSD of a signal?

Hi

i have a signal recorded at 125Hz. I have preprocessed it with zero mean and buterworth filter 0.1 - 0.4Hz as i need to calculate repsiration rate by detcting the highest frequency linked to the highest psd value in it.

I am running the following code to calculate respiratory rate by usingWelch’s method which divides the interbeat intervals (signal) into overlapping segments and computes a modified periodogram for each segment. Then anaverage of all periodograms along with an array of frequencies are returned at the output. The respiratory rate is themaximum frequency within the frequency band (Hz) that is,where PSD is maximum. To calculate the final respiratoryate the frequency is multiplied by 60 to convert therespiratory frequency band from Hz to breaths per minute. But the code returns the number of segments to be -1 and the repiration rate as 0. I am attaching the filtered signal. please help!

samplingFrequency = 125; % Sampling frequency in Hz

% Step 1: Identify the peaks (using the findpeaks function)

[peaks, peakLocations] = findpeaks(filteredData, 'MinPeakProminence',0.15);

% Step 2: Calculate the time difference between consecutive peaks

IBIs = diff(peakLocations) / samplingFrequency; % Calculate the time differences in seconds

% Step 3: Compute periodogram for each IBI segment with overlapping

segmentLength = round(mean(IBIs) * samplingFrequency); % Length of each segment

overlap = 0.5; % Overlap percentage between segments (50% overlap)

segmentSamples = floor(segmentLength * (1 - overlap)); % Number of samples in each segment

numSegments = floor((length(IBIs) - segmentLength) / segmentSamples) + 1; % Calculate the number of segments

% Initialize variables to store periodograms

periodograms = zeros(segmentLength, numSegments);

% Compute periodogram for each IBI segment

for i = 1:numSegments

startIndex = (i - 1) * segmentSamples + 1;

endIndex = startIndex + segmentLength - 1;

segment = filteredData(peakLocations(i):peakLocations(i+1));

[~, freq, periodograms(:, i)] = periodogram(segment, hann(segmentLength), segmentLength, samplingFrequency, 'power', 'psd');

if i == 1

frequencies = freq

end

end

frequencies

% Step 4: Compute average periodogram

averagePeriodogram = mean(periodograms, 2);

% Step 5: Calculate respiratory rate

[~, maxIndex] = max(averagePeriodogram);

maxFrequency = frequencies(maxIndex);

respiratoryRate = maxFrequency * 60; % Convert frequency from Hz to breaths per minute

% Display the respiratory rate

disp(['Respiratory Rate: ', num2str(respiratoryRate), ' breaths per minute']);

[Mathieu NOE]

seems to me the correct way of computing your amount of segments is

numSegments = floor((length(filteredData) - segmentLength) / segmentSamples) + 1; % Calculate the number of segments

[Nyah]

@Mathieu NOE Thanks for the correction. I have to use the window length of 60s. as the overall data is of around 8 minutes. I am using this code to calculate repiration rate now, but i am not sure if i am doing it correctly. an insight will be highly appreciated.

samplingFrequency = 125; % Sampling frequency in Hz

windowLength = 60; % Window length in seconds

% Step 1: Identify the peaks (using the findpeaks function)

[peaks, peakLocations] = findpeaks(filteredData,'MinPeakProminence',0.15);

% Step 2: Calculate the time difference between consecutive peaks

IBIs = diff(peakLocations) / samplingFrequency; % Calculate the time differences in seconds

% Step 3: Compute periodogram for each IBI segment with a 60-second window length

segmentLength = windowLength * samplingFrequency; % Length of each segment

overlap = 0.5; % Overlap percentage between segments (50% overlap)

segmentSamples = floor(segmentLength * (1 - overlap)); % Number of samples in each segment

numSegments = floor((length(filteredData) - segmentLength) / segmentSamples) + 1; % Calculate the number of segments

% Initialize variables to store periodograms

periodograms = [];

% Compute periodogram for each IBI segment

for i = 1:numSegments

startIndex = 1 + (i - 1) * segmentSamples;

endIndex = startIndex + segmentLength - 1;

if endIndex > length(filteredData)

endIndex = length(filteredData);

end

segment = filteredData(startIndex:endIndex);

[~, freq, periodogramData] = periodogram(segment, hann(length(segment)), [], samplingFrequency);

if isempty(periodograms)

periodograms = periodogramData;

frequencies = freq;

else

periodograms = periodograms + periodogramData;

end

end

% Average the periodograms

averagePeriodogram = periodograms / numSegments;

% Step 5: Calculate respiratory rate

[~, maxIndex] = max(averagePeriodogram);

maxFrequency = frequencies(maxIndex);

respiratoryRate = maxFrequency * 60; % Convert frequency from Hz to breaths per minute

% Display the respiratory rate

disp(['Respiratory Rate: ', num2str(respiratoryRate), ' breaths per minute']);

[Mathieu NOE]

hello again

I made a few corrections , seems to work now :

samplingFrequency = 125; % Sampling frequency in Hz
windowLength = 60; % Window length in seconds
% Step 1: Identify the peaks (using the findpeaks function)
[peaks, peakLocations] = findpeaks(filteredData,'MinPeakProminence',0.15);
% Step 2: Calculate the time difference between consecutive peaks
IBIs = diff(peakLocations) / samplingFrequency; % Calculate the time differences in seconds
% Step 3: Compute periodogram for each IBI segment with a 60-second window length
segmentLength = windowLength * samplingFrequency; % Length of each segment
overlap = 0.5; % Overlap percentage between segments (50% overlap)
segmentSamples = floor(segmentLength * (1 - overlap)); % Number of samples in each segment
numSegments = floor((length(filteredData) - segmentLength) / segmentSamples) + 1; % Calculate the number of segments
% Initialize variables to store periodograms
periodograms = [];
samples = numel(filteredData);
periodograms = 0;

% Compute periodogram for each IBI segment
for i = 1:numSegments
    startIndex = 1 + (i - 1) * segmentSamples;
    endIndex = startIndex + segmentLength - 1;
    
    segment = filteredData(startIndex:endIndex);
    
    if endIndex > samples
        % pad data with zeros
        segment= padarray(segment,segmentLength);
    end
    
    [periodogramData, freq] = periodogram(segment, hann(segmentLength), [], samplingFrequency);
    periodograms = periodograms + periodogramData;
    
 end
% Average the periodograms
averagePeriodogram = periodograms / numSegments;
% Step 5: Calculate respiratory rate
[~, maxIndex] = max(averagePeriodogram);
maxFrequency = freq(maxIndex);
respiratoryRate = maxFrequency * 60; % Convert frequency from Hz to breaths per minute
% Display the respiratory rate
disp(['Respiratory Rate: ', num2str(respiratoryRate), ' breaths per minute']);