How to plot a frequency of respiration signal on y-axis while time of measurement is on x-axis?

I have a respiration (breathing) signal and a sampling frequency of 25 Hz and need to detect where is the lowest breathing frequency on a time scale, which should tell me actually when the person became sleepy. Fourier transform in its classical form doesn't give me much useful information. So, to clarify: the time of measurement should be on the x-axis and the breathing frequency should be on the y-axis. Then, I suppose, lower amplitudes of the signal will show the slower breathing. What should be done with the signal to plot it the way I need?

[Birdman]

Can you share your data in a mat file?

[miki90]

Here is also a vector in a txt file, if it is ok.

Answer 1

[Star Strider]

This should get you started:

D = load('respiratory.txt');
Fs = 25;                                                    % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency
Ts = 1/Fs;                                                  % Sampling Time (sec)
L = numel(D);
t = linspace(0, L, L)*Ts;                                   % Time Vector (sec)

figure(1)
plot(t, D)
grid
% axis([0  60    -850  -750])
axis([xlim    -850  -750])
xlabel('Time')
ylabel('Amplitude')

FTD = fft(D-mean(D))/L;                                     % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                         % Frequency Vector
Iv = 1:numel(Fv);                                           % Index Vector

figure(2)
plot(Fv, abs(FTD(Iv))*2)
grid
axis([0  2.5    ylim])
xlabel('Frequency (Hz)')
ylabel('Amplitude')

Wp = [0.35 0.65]/Fn;                                        % Passband Frequency (Normalised)
Ws = [0.30 0.75]/Fn;                                        % Stopband Frequency (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs =  50;                                                   % Stopband Ripple (dB)
[n,Ws]  = cheb2ord(Wp,Ws,Rp,Rs);                            % Filter Order
[z,p,k] = cheby2(n,Rs,Ws);                                  % Filter Design, Sepcify Bandpass
[sos,g] = zp2sos(z,p,k);                                    % Convert To Second-Order-Section For Stability

figure(3)
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot

D_filtered = filtfilt(sos, g, D);                           % Filter Signal

[pks,locs] = findpeaks(D_filtered, 'MinPeakDist',40);

figure(4)
plot(t, D_filtered)
hold on
plot(t(locs), pks, '^r')
hold off
grid
% axis([0  60    ylim])
axis([0  60    -15  15])
xlabel('Time')
ylabel('Amplitude')

tdif = diff([0 t(locs)]);                                   % Time Difference Between Peaks (sec)
Dfrq = 60./tdif;                                            % Frequency (Respirations/Minute)

figure(5)
plot(t(locs), Dfrq)
grid
axis([xlim    10  40])
xlabel('Time (sec)')
ylabel('Frequency (minute^{-1})')

I do not intend this to be a complete solution. It does some signal processing and preliminary data analysis. Make necessary changes to get the result you want. See the documentation on the various functions to understand how the code works.

Have fun!

[miki90]

This is more than I needed. Thank you so much for a detailed answer. I cited you on Stack Overflow.

[Star Strider]

As always, my pleasure!

Thank you!

[Thomas]

Hi, I hope you still get to read this thread: I want to apply the same code you wrote, but instead of to the loaded respiratory.txt data, to a real time signal I am reading from an Arduino:

clc;

clear all;

close all;

a = arduino();

clear a;

a = arduino("COM4", "Uno");

writeDigitalPin(a, "D12", 1); %LED is turned on when "1" is written

Fs = 10;

T = 1/Fs;

L=1000;

t=(0:L-1)*T;

x = 0;

for k = 1:1000

b = readVoltage(a, "A0")

x = [x, b]

plot (x);

ylabel('Voltage in volts')

xlabel('Samples')

ylim([3 4])

grid;

drawnow;

end

How do I have to change the first lines of code then?

[miki90]

I don't have Matlab on my PC right now but I will try to answer like this. What do you get when you just plot your real-time signal? If you want to get the result as from the file, you may try to put your real time data into an array (make a buffer) and to plot the chunks, one after another.

[Star Strider]

Thomas — I don’t have Arduino, so I can’t help you with this. You need to post this as a new Question, then delete it here.