Extracting a respiratory signal from a raw micro-CT-based signal

Hi,

I'm new in signal processing. I want extract a respiratory signal from a raw micro-CT-basedsignal, which is non-uniformly sampled in time domain. The sample rate of the signal is 3.5 Hz and based on Lomb-Scargle periodogram I know that the respiratory frequency is 1.73 Hz. Although the signal is non-uniformly sampled, still with the bandpass function I can extract the respiratory signal and calculate the respiratory rate with a good accuracy. But I need to extract the signal with the actual amplitude in the following form:

While I see my extracted signal in the following form:

I tried to use Chebyshev type-II filter design to bandpass the signal non-uniformly, but I cannot get the same respiratory rate as the bandpass function. I'll attach my code and raw data. I have the following questions:

1. Can I use Chebyshev type-II filter with low sample rate signal (3.5 Hz)? Now when I normalize my passband and stopband frequency ranges, the upper values go over 1. Can I normalize them with my sample rate?
2. I have no idea what values should I use for my passband and stopband ripple.
3. How should I get to the respiratory signal with the frist pic shape? I want to compare the amplitude of different respiratory signals.

this is my code:

%% Design a bandpass Chebyshev type-II filter to extract respiratory signal from raw signals comming from micro-CT images

clear;

clc;

%% Import the raw signals

% The original raw signal: voxel intensities

OrigSig=load('Kaat_m22__RawIntensity.csv');

% The flattend raw signal

FlatSig=load("Kaat_m22__FlattendIntensity.csv");

%% Import the timemarks in millisecond format

% Timemarks (millisec)

t_hat=load ("Kaat_TimeMarks.csv");

N=length(t_hat);

t_hatt=zeros(1,N);

for i=2:N

if t_hat(i)>=t_hat(i-1)

t_hatt(i)=t_hatt(i-1)+(t_hat(i)-t_hat(i-1));

else

t_hatt(i)=t_hatt(i-1)+((1000+t_hat(i))-t_hat(i-1));

end

end

t_hattt=t_hatt(:);

% Final timemarks (sec)

t=t_hattt*0.001;

%% Calculate the periodogram using the Lomb-Scargle method and find the respiratory signal

%Lomb-Scargle method

[Pow,f]=plomb(FlatSig,t);

% Find the dominant frequency corresponding to the respiratory signal

[maxP, index]=max(Pow);

RespFrq=f(index);

%% Chebyshev type-2 filter design

% Sampling frequency (Hz)

Fs=N/t(N);

% Nyquist frequency

Fn=Fs/2;

% Specify the lower cutoff frequency

Lowpass=abs(RespFrq-0.2);

% Specify the upper cutoff frequency

Highpass=abs(RespFrq+0.2);

% Normalize the passband frequncy

Wp=[Lowpass Highpass]/Fs;

% Normalize the stopband frequency

Ws=[Lowpass-0.05 Highpass+0.05]/Fs;

% Passband ripple (dB)

Rp=1;

% Stopband ripple (dB)

Rs=60;

% Filter order

[n,Ws]=cheb2ord(Wp,Ws,Rp,Rs);

% Filter design, specify bandpass

[z,p,k]=cheby2(n,Rs,Ws);

% Convert to second-order-section for stability

[sos,g]=zp2sos(z,p,k);

% Filter Bode plot

figure(1)

freqz(sos, 2^16, Fs);

% Filter the raw signal to find the respiratory signal

RespSig=filtfilt(sos, g, FlatSig);

% Find the respiratory peaks

[pks,locs]=findpeaks(RespSig);

%% Plot the results

figure(2)

% plot the raw signal in time domain

subplot(4,1,1)

plot(t,OrigSig);

xlabel('Time (sec)');

ylabel('Voxel intensities');

grid

% Plot the power spectrum in frequency domain

subplot(4,1,2)

plot(f, Pow)

xlabel('Frequency (Hz)')

ylabel('Power (dB)')

hold on

plot(f(index),maxP, '^',...

'Color',[217 83 25]/255, 'MarkerFaceColor',[217 83 25]/255)

hold on

title(['Respiratory frequency (Hz): ', num2str(f(index))]);

hold off

grid

subplot(4,1,3);

plomb(FlatSig,t);

grid

% Plot the respiratory signal and peaks in time domain

subplot(4,1,4)

plot(t,RespSig);

hold on

plot(t(locs),pks, '^',...

'Color',[217 83 25]/255, 'MarkerFaceColor',[217 83 25]/255)

hold off

xlabel('Time (sec)')

ylabel('Amplitude')

title(['Number of extrema: ', num2str(nnz(pks))]);

grid

Thanks.

Answer 1

[Garmit Pant]

Hello Kaveh Ahookhosh

It is my understanding that you want to extract the respiratory signal from a raw micro-CT-based signal that has been non-uniformly sampled and you are not getting the desired output on applying the bandpass filter designed by you on the input data.

Based on my investigation, I would like to provide you with the following information regarding your queries.:

1. Yes, you can use a “Chebyshev type-II filter” with a low sample rate signal.The method of normalization used in your code is correct.
2. The choice of passband and stopband ripple depends on the requirements of your application. The passband ripple is the amount of variation in the amplitude, within the designated passband of the filter, and stop band attenuation is the minimum attenuation level with the designated rejection band of the filter. Choose your passband and stopband ripple accordingly.
3. Your original data is non-uniformly sampled. You can resample your data to make it uniform before applying the filter. Refer to the following MATLAB Documentation to use “Signal Analyzer” app to resample your data:

https://www.mathworks.com/help/signal/ug/resample-and-filter-signal-with-missing-samples.html

Changing the specifications of your “Chebyshev type-II filter” can also produce different results.

I hope this provides you with enough information to resolve the issues you are facing.

Regards

Garmit