Extracting Respiratory Cycle Duration from RPM Waveform Stored in DICOM Files

60 次查看（过去 30 天）
W031997
W031997 2025-9-4，7:12
回答： Star Strider 2025-9-5，11:19
Dear MATLAB Community,
I am currently working with respiratory motion waveform data acquired from a Real-Time Position Management (RPM) system, stored in DICOM files (MW.dcm). The waveform is represented in 8-bit format, with values ranging from 0 to 255. I have written a MATLAB script to plot the raw signal and applied a band-pass filter (removing values <10 and >250, with a passband of 0.2–0.4 Hz) in an effort to denoise the signal and obtain a smooth, sinusoidal-like waveform suitable for identifying individual breathing cycles.
Additionally, I attempted to use the Signal Analyzer tool in MATLAB to visualize the waveform, which is shown in Figure 3. Despite these efforts, the waveform remains significantly noisy, and the calculated mean duration of a single breathing cycle does not align with physiological expectations. For reference, a typical respiratory cycle lasts approximately 3–5 seconds, with each cycle resembling a sinusoidal wave.
I would greatly appreciate any guidance from the community on how to proceed from this point. Specifically, I am seeking reliable methods for extracting the duration of individual breathing cycles from RPM DICOM files. Insights, best practices, or illustrative MATLAB examples from those who have experience with this type of waveform data would be extremely valuable.
Thank you very much for your time and assistance.
.
  5 个评论
显示 3更早的评论
Star Strider
Star Strider 2025-9-4，13:21
@W031997 -- It would help to have a representative sample of your data (raw signal), and the sampling times. (Sampling times are preferable to the sampling frequency, since the time intervals may not be constant.)
Are the 'values ranging from 0 to 255' unsigned integers or signed integers?
W031997
W031997 2025-9-5，6:52
I have a dataset of unsigned values 0–255 from a machine, sampled at 25 Hz.This dataset represents the raw waveform data I obtained from the machine. I’m relatively new to signal processing and was wondering whether it is possible to extract timing information from this data. thank you

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2025-9-5，11:19
Ran in:
Part of the problem is that your instrument is 'saturating', such that the amplitude of the signal exceeds the ability of the instrumentation to accommodate it. This leads to 'flat-topping' and that will distort the spectrum.
I am not certain what you want to do with these data. My approach here filters the data using information from the Fourier Transform (you can of course change the filter passband), then detects the threshold-crossings and plots them. My code also calculates the direction the filtered signal is going at the threshold-crossing. (I set 'thrshld' to zero here because a bandpass filter will filter out the D-C or constant offsets. If you want to use a lowpass filter instead, that will retain the D-C offset, so you will need to set 'thrshld' to a value that makes sense in that context.) It also calculates the signal value at the threshold crossing. This should be close to the value of the threshold, and I include it as a check on that. The threshold crossings with their 'Direction' value can be used to define the elements of the respiratory signal. I am not certain how you want to do that, so I will stop here for now. I will help with subsequent code in order to get all this working to your satisfaction.
Try this --
LD = load('rpm_L01.mat');
s = LD.signal_raw;
Fs = 25; % Hz
L = numel(s); % Signal Length
t = linspace(0, L-1, L).'/Fs; % Time Vector
figure
plot(t, s)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Raw Signal')
xlim([0 10])
[FTs1,Fv] = FFT1(s,t);
[pks,locs] = findpeaks(abs(FTs1)*2, MinPeakProminence=5);
Peaks_Freqs = table(Fv(locs),pks,VariableNames=["Frequency","Magnitude"])
Peaks_Freqs = 7×2 table
Frequency Magnitude _________ _________ 1.5625 132.71 3.125 60.184 4.6875 14.336 6.25 18.04 7.8125 9.2014 9.375 21.796 10.938 18.75
% Harmonics = diff(PeakFreqs{:,1})
figure
plot(Fv, abs(FTs1)*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform')
●
s_filt = bandpass(s, [0.5 2.5], Fs, ImpulseResponse='iir'); % IIR Bandpass Filter
thrshld = 0; % Threshold
zxi = find(diff(sign(s_filt-thrshld))); % Approximate Indices Of Threshold-Crossings
for k = 1:numel(zxi)
idxrng = max(1, zxi(k)-1) : min(L,zxi(k)+1);
tv(k,:) = interp1(s_filt(idxrng), t(idxrng), 0);
sv(k,:) = interp1(t(idxrng), s_filt(idxrng), tv(k));
dsdt(k,:) = median(s_filt(idxrng)./t(idxrng));
end
Threshold_Crossings = table(tv, sv, dsdt, VariableNames=["Time","Amplitude","Direction"])
Midpoint_Crossings = 32430×3 table
Time Amplitude Direction ________ ___________ _________ 0.018357 1.7764e-15 -Inf 0.62802 2.1316e-14 48.163 0.92922 2.3093e-14 -5.6854 1.1601 2.9379e-14 0.058945 1.4918 -6.0396e-14 -6.1388 1.8717 1.1635e-13 15.346 2.2253 -1.7764e-14 -10.033 2.5121 8.8818e-14 11.471 2.8363 -8.8818e-16 -9.9154 3.1799 -1.7941e-13 5.744 3.4868 6.2172e-15 -2.3447 3.8066 -2.6734e-13 1.5522 4.1213 3.3307e-15 -0.19829 4.4459 -8.4377e-14 1.3013 4.7959 5.1692e-13 -8.3874 5.1154 -2.4736e-13 6.0433
figure
plot(t, s_filt, DisplayName="Filtered Signal")
hold on
plot(tv, sv, 'xr', DisplayName='Threshold Crossings')
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Filtered Signal')
legend(Location="best")
xlim([0 10])
function [FTs1,Fv] = FFT1(s,t)
% One-Sided Numerical Fourier Transform
% Arguments:
% s: Signal Vector Or Matrix
% t: Associated Time Vector
t = t(:);
L = numel(t);
if size(s,2) == L
s = s.';
end
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
Fv = Fs*(0:(NFFT/2))/NFFT;
% Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
Fv = Fv(:);
FTs1 = FTs(Iv,:);
end
.

类别

在 Help Center 和 File Exchange 中查找有关 Audio I/O and Waveform Generation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by