Calculating the aortic flow downstream valve using Womersley and 3-element Windkessel model calculations from my FSI results - Why unusual sine waves at the end of flow curve?

So I have some aortic valve FSI results and I want to examine these results by modelling Womersley flow and using a 3-element Windkessel model in series to it (both downstream of the valve).

I am using the pressure (at the STJ) just downstream of my valve for the calculations (P_STJ, provided in the attached csv file).

I calculate the Womersley Impedance as according to literature (https://leifh.folk.ntnu.no/teaching/tkt4150/._main023.html) and it is connected to a 3 element windkessel model in series (as you can see in the code, with the values provided).

I perform Fourier transform on the pressure data and then divide it by total impedance to get the flow. However, as you can see when running the code (and attached image), that there is a discrepancy in values between the FSI values, that is approximately equal to the amplitude of that unusual sine wave at the end of the flow curve.

Am I doing something wrong? Should I attempt another method than fourier transform and series? Any recommendations?

Thank you in advance.

[Hussam]

If any additional information is needed, please let me know.

Answer 1

[Star Strider]

I’m not familiar with your abbreviations.

FSI = ?

STJ = ?

Beyond that, your code runs without error.

Please explain what your code does, since that is not obvious to me.

(1.) What specific problem are you having?

(2.) Is there a reason you are using your own Fourier transform routine rather than the fft function?

The two curves in the supplied image appear to me to be reasonable approximations to each other. If one is the inverse Fourier transform of the original time-domain signal, it may be necessary to examine the functions you used to do the Fourier and inverse Fourier transforms, since they should be the same if you are using all the available frequencies. The imaginary part of the inverse transform should be zero (or close to it). In that instance, see the documentation for the ifft function, and note the difference between the 'symmetric' and 'nonsymmetric' results.

imshow(imread('Flow_Comp.PNG'))

% clear all
% close all
% clc

%% Input Data

data = readtable('ES_N_00_data_even.csv');
P_STJ_ref = data.P_STJ;
P_STJ_mmHg_ref = 0.00750062*P_STJ_ref;
t_ref = data.t;
Q_FSI_ref = data.Q;

% Interpolating and optional smoothing

dt = 1e-4;
t=(0:dt:t_ref(end))';

P_STJ_interp = interp1(t_ref, P_STJ_ref, t);
P_STJ_smooth = smooth(P_STJ_interp, 50);
P_STJ = P_STJ_interp;
P_STJ_mmHg = 0.00750062*P_STJ;

Q_FSI_interp = interp1(t_ref, Q_FSI_ref, t);
Q_FSI_smooth = smooth(Q_FSI_interp, 50);
Q_FSI = Q_FSI_interp;

% Aorta 3-Element Windkessel Parameters
R1 = 1.1224e+06;    % Proximal Resistance
R2 = 1.6052e+08;    % Distal Resistance
C = 9.2237e-09;     % Capacitance

r_inlet = 0.0215/2;     % Radius - Inlet
r_00 = 0.0255/2;        % Radius - Baseline Outlet

l = 120e-3;             % Inlet to Outlet Length
l_asc = 80e-3;          % Ascending Aorta Tubular Region Length
l_LVOT = 20e-3;         % LVOT Region Length

mu = 0.0035;            % Viscosity
rho = 1060;             % Density
t_cycle = 0.854;        % Cycle Time

%% Poiseuille Flow Calculation for Baseline

R_asc_00 = (8*mu*l_asc) / (pi*(r_00^4));
R_AoWK = R1 + R2;

% Total resistance for Poiseuille flow
R_total = R_asc_00 + R_AoWK;

%% Plotting FT for P and Q_00

P_STJ_FT = fft(P_STJ);
P_STJ_FT_shift = fftshift(P_STJ_FT);
P_STJ_FT_abs = abs(P_STJ_FT_shift);
P_STJ_FT_angle = angle(P_STJ_FT_shift);


Q_FSI_FT = fft(Q_FSI);
Q_FSI_FT_shift = fftshift(Q_FSI_FT);
Q_FSI_FT_abs = abs(Q_FSI_FT_shift);
Q_FSI_FT_angle = angle(Q_FSI_FT_shift);

Ts = mean(diff(t));
Ts_std = std(Ts);
fs = 1/Ts;
N = length(t);
f_axis = [0:N-1]'*fs/N;
f_shift = (-N/2:N/2-1)*(fs/N);

f_range = 100;
f_0 = N/2 + 1;
f_range_idx = f_0 + f_range;

figure(1)
subplot(2,2,1)
% stem(f_axis, P_STJ_FT_abs, 'r.', 'MarkerSize', 15)
stem(f_shift, P_STJ_FT_abs, 'r.', 'MarkerSize', 15)
title('STJ Pressure Fourier Transform Magnitude vs Frequency')
xlabel('Frequency (Hz)')
ylabel('|FT|')

subplot(2,2,3)
stem(f_shift(f_0:f_range_idx), P_STJ_FT_abs(f_0:f_range_idx), 'r.', 'MarkerSize', 15)
title('STJ Pressure Fourier Transform Magnitude vs Frequency')
xlabel('Frequency (Hz)')
ylabel('|FT|')
ylim([0 1e7])

subplot(2,2,2)
stem(f_shift, P_STJ_FT_angle, 'r.', 'MarkerSize', 15)
title('STJ Pressure Fourier Transform Angle vs Frequency')
xlabel('Frequency (Hz)')
ylabel('Phase Angle (rad)')

subplot(2,2,4)
stem(f_shift(f_0:f_range_idx), P_STJ_FT_angle(f_0:f_range_idx), 'r.', 'MarkerSize', 15)
title('STJ Pressure Fourier Transform Angle vs Frequency')
xlabel('Frequency (Hz)')
ylabel('Phase Angle (rad)')
ylim([-pi pi])

figure(2)
subplot(2,2,1)
% stem(f_axis, Q_FSI_FT_abs, 'r.', 'MarkerSize', 15)
stem(f_shift, Q_FSI_FT_abs, 'r.', 'MarkerSize', 15)
title('Q 00 Flow Fourier Transform Magnitude vs Frequency')
xlabel('Frequency (Hz)')
ylabel('|FT|')

subplot(2,2,3)
stem(f_shift(f_0:f_range_idx), Q_FSI_FT_abs(f_0:f_range_idx), 'r.', 'MarkerSize', 15)
title('Q 00 Flow Fourier Transform Magnitude vs Frequency')
xlabel('Frequency (Hz)')
ylabel('|FT|')
% ylim([0 1e7])

subplot(2,2,2)
stem(f_shift, Q_FSI_FT_angle, 'r.', 'MarkerSize', 15)
title('Q 00 Flow Fourier Transform Angle vs Frequency')
xlabel('Frequency (Hz)')
ylabel('Phase Angle (rad)')

subplot(2,2,4)
stem(f_shift(f_0:f_range_idx), Q_FSI_FT_angle(f_0:f_range_idx), 'r.', 'MarkerSize', 15)
title('Q 00 Flow Fourier Transform Angle vs Frequency')
xlabel('Frequency (Hz)')
ylabel('Phase Angle (rad)')
ylim([-pi pi])

%% Plotting Impedances

% Womersley for freq range
f_set = 0:0.001:100;
for i=1:length(f_set)

    Z_asc_set(i) = l_asc*womersley_impedance(r_00, rho, mu, f_set(i));

    Z_AoWK_set(i) = R1 + ( R2 / ( 1 + 1i*C*R2*(2*pi*f_set(i)) ) );

    % Q_set(i) = real( (P_Inlet - P_Outlet) / Z_set(i));
    % Q_set(i) = real( (dP_systole_avg) / Z_set(i));
    % Q_new_set(i) = real( (P_Inlet - P_Outlet) / Z_new_set(i));
    % Q_new_set(i) = real( (dP_systole_avg) / Z_new_set(i));
end

Z_total_set = Z_asc_set + Z_AoWK_set;

figure(3)

subplot(2,2,1)
hold on
plot(0,R_asc_00,'o','LineWidth',3,'DisplayName',strcat('AA Poiseuille Resistance, d_{AA} = 25.5 mm'))
plot(f_set(2:end),abs(Z_asc_set(2:end)),'DisplayName',strcat('AA Womersley Impedance, d_{AA} = 25.5 mm'))
lgd = legend('Location','southeast');
% fontsize(lgd,'decrease')
title('AA Womersley Impedance Magnitude vs Frequency')
xlabel('Frequency (Hz)')
ylabel('|Z|')
grid on
hold off

subplot(2,2,2)
hold on
plot(0,R_AoWK,'o','LineWidth',3,'DisplayName',strcat('Aortic 3-Element WK Poiseuille Resistance'))
plot(f_set(2:end),abs(Z_AoWK_set(2:end)),'DisplayName',strcat('Aortic 3-Element WK Impedance'))
lgd = legend('Location','northeast');
% fontsize(lgd,'decrease')
title('Aortic WK Impedance Magnitude vs Frequency')
xlabel('Frequency (Hz)')
ylabel('|Z|')
grid on
hold off

subplot(2,2,3)
hold on
plot(0,R_total,'o','LineWidth',3,'DisplayName',strcat('Total Poiseuille Resistance, d_{AA} = 25.5 mm'))
plot(f_set(2:end),abs(Z_total_set(2:end)),'DisplayName',strcat('Total Impedance, d_{AA} = 25.5 mm'))
lgd = legend('Location','northeast');
% fontsize(lgd,'decrease')
title('Total Impedance Magnitude vs Frequency')
xlabel('Frequency (Hz)')
ylabel('|Z|')
grid on
hold off

subplot(2,2,4)
hold on
plot(0,R_total,'o','LineWidth',3,'DisplayName',strcat('Total Poiseuille Resistance, d_{AA} = 25.5 mm'))
plot(f_set(2:end),abs(Z_total_set(2:end)),'DisplayName',strcat('Total Impedance, d_{AA} = 25.5 mm'))
lgd = legend('Location','northeast');
% fontsize(lgd,'decrease')
title('Total Impedance Magnitude vs Frequency')
xlabel('Frequency (Hz)')
ylabel('|Z|')
ylim([1.2e7 1.5e7])
xlim([1.1 1.3])
grid on
hold off

figure(4)

subplot(2,2,1)
hold on
% plot(0,R_asc_00,'o','LineWidth',3,'DisplayName',strcat('AA Poiseuille Resistance, d_{AA} = 25.5 mm'))
% plot(0,R_asc_03,'o','LineWidth',3,'DisplayName',strcat('AA Poiseuille Resistance, d_{AA} = 21.5 mm'))
plot(f_set(2:end),angle(Z_asc_set(2:end)),'DisplayName',strcat('AA Womersley Impedance, d_{AA} = 25.5 mm'))
lgd = legend('Location','southeast');
% fontsize(lgd,'decrease')
title('AA Womersley Impedance Phase Angle vs Frequency')
xlabel('Frequency (Hz)')
ylabel('Phase Angle (rad)')
grid on
hold off

subplot(2,2,2)
hold on
% plot(0,R_AoWK,'o','LineWidth',3,'DisplayName',strcat('Aortic 3-Element WK Poiseuille Resistance'))
plot(f_set(2:end),angle(Z_AoWK_set(2:end)),'DisplayName',strcat('Aortic 3-Element WK Impedance'))
lgd = legend('Location','northeast');
% fontsize(lgd,'decrease')
title('Aortic WK Impedance Phase Angle vs Frequency')
xlabel('Frequency (Hz)')
ylabel('Phase Angle (rad)')
grid on
hold off

subplot(2,2,3)
hold on
% plot(0,R_total_sys,'o','LineWidth',3,'DisplayName',strcat('Total Poiseuille Resistance, d_{AA} = 25.5 mm'))
% plot(0,R_total_03_sys,'o','LineWidth',3,'DisplayName',strcat('Total Poiseuille Resistance, d_{AA} = 21.5 mm'))
plot(f_set(2:end),angle(Z_total_set(2:end)),'DisplayName',strcat('Total Impedance, d_{AA} = 25.5 mm'))
lgd = legend('Location','northeast');
% fontsize(lgd,'decrease')
title('Total Impedance Phase Angle vs Frequency')
xlabel('Frequency (Hz)')
ylabel('Phase Angle (rad)')
grid on
hold off

subplot(2,2,4)
hold on
% plot(0,R_total_sys,'o','LineWidth',3,'DisplayName',strcat('Total Poiseuille Resistance, d_{AA} = 25.5 mm'))
% plot(0,R_total_03_sys,'o','LineWidth',3,'DisplayName',strcat('Total Poiseuille Resistance, d_{AA} = 21.5 mm'))
plot(f_set(2:end),angle(Z_total_set(2:end)),'DisplayName',strcat('Total Impedance, d_{AA} = 25.5 mm'))
lgd = legend('Location','northeast');
% fontsize(lgd,'decrease')
title('Total Impedance Phase Angle vs Frequency')
xlabel('Frequency (Hz)')
ylabel('Phase Angle (rad)')
ylim([-1.4 -1.3])
% xlim([1.1 1.3])
grid on
hold off

%% Fourier Series for STJ Pressure

[t_FS_Full, P_STJ_FS_Full] = compute_fourier_series(t, P_STJ_FT_shift, f_shift, N, f_0, fs/2);
P_STJ_FS_Full_mmHg =  0.00750062*P_STJ_FS_Full;

[t_FS_30, P_STJ_FS_30] = compute_fourier_series(t, P_STJ_FT_shift, f_shift, N, f_0, 30);
P_STJ_FS_30_mmHg =  0.00750062*P_STJ_FS_30;

[t_FS_40, P_STJ_FS_40] = compute_fourier_series(t, P_STJ_FT_shift, f_shift, N, f_0, 40);
P_STJ_FS_40_mmHg =  0.00750062*P_STJ_FS_40;

[t_FS_85, P_STJ_FS_85] = compute_fourier_series(t, P_STJ_FT_shift, f_shift, N, f_0, 85);
P_STJ_FS_85_mmHg =  0.00750062*P_STJ_FS_85;

figure
subplot(2,2,1)
plot(t,P_STJ_mmHg,'DisplayName',strcat('FSI Data'))
title('STJ Pressure from FSI')
xlabel('Time (s)')
ylabel('Pressure (mmHg)')
ylim([0 180])

subplot(2,2,2)
hold on
plot(t,P_STJ_mmHg,'DisplayName',strcat('FSI Data'))
plot(t_FS_Full,real(P_STJ_FS_Full_mmHg(1:length(t_FS_Full))),'DisplayName',strcat('Fourier Series, All Frequencies'))
title('STJ Pressure from Fourier Series (All Frequencies)')
xlabel('Time (s)')
ylabel('Pressure (mmHg)')
ylim([0 180])
legend
hold off

subplot(2,2,3)
hold on
plot(t_FS_30,real(P_STJ_FS_30_mmHg(1:length(t_FS_30))),'DisplayName',strcat('Fourier Series, f_{max} = 30 Hz'))
plot(t_FS_40,real(P_STJ_FS_40_mmHg(1:length(t_FS_40))),'DisplayName',strcat('Fourier Series, f_{max} = 40 Hz'))
plot(t_FS_85,real(P_STJ_FS_85_mmHg(1:length(t_FS_85))),'DisplayName',strcat('Fourier Series, f_{max} = 85 Hz'))
title('STJ Pressure from Fourier Series')
xlabel('Time (s)')
ylabel('Pressure (mmHg)')
legend
ylim([0 180])
hold off

subplot(2,2,4)
hold on
plot(t,P_STJ_mmHg,'DisplayName',strcat('FSI Data'))
plot(t_FS_85,real(P_STJ_FS_85_mmHg(1:length(t_FS_85))),'DisplayName',strcat('Fourier Series, f_{max} = 85 Hz'))
title('STJ Pressure from Fourier Series (Up to Freq = 85 Hz)')
xlabel('Time (s)')
ylabel('Pressure (mmHg)')
ylim([0 180])
legend
hold off

%% Fourier Series for Q 00 Flow (from FSI)

[t_FS_Full, Q_FSI_FS_Full] = compute_fourier_series(t, Q_FSI_FT_shift, f_shift, N, f_0, fs/2);
% Q_FSI_FS_Full_mmHg =  0.00750062*Q_FSI_FS_Full;

[t_FS_32, Q_FSI_FS_32] = compute_fourier_series(t, Q_FSI_FT_shift, f_shift, N, f_0, 32);
% Q_FSI_FS_30_mmHg =  0.00750062*Q_FSI_FS_30;

[t_FS_50, Q_FSI_FS_50] = compute_fourier_series(t, Q_FSI_FT_shift, f_shift, N, f_0, 50);
% Q_FSI_FS_40_mmHg =  0.00750062*Q_FSI_FS_40;

[t_FS_85, Q_FSI_FS_85] = compute_fourier_series(t, Q_FSI_FT_shift, f_shift, N, f_0, 85);
% Q_FSI_FS_85_mmHg =  0.00750062*Q_FSI_FS_85;

figure
subplot(2,2,1)
plot(t,Q_FSI)
title('Q Flow from FSI')
xlabel('Time (s)')
ylabel('Flow (m^3/s)')
% ylim([0 180])

subplot(2,2,2)
hold on
plot(t,Q_FSI,'DisplayName',strcat('FSI Data'))
plot(t_FS_Full,real(Q_FSI_FS_Full(1:length(t_FS_Full))),'DisplayName',strcat('Fourier Series, All Frequencies'))
title('Q Flow from Fourier Series (All Frequencies)')
xlabel('Time (s)')
ylabel('Flow (m^3/s)')
% ylim([0 180])
hold off

subplot(2,2,3)
hold on
plot(t_FS_32,real(Q_FSI_FS_32(1:length(t_FS_32))),'DisplayName',strcat('Fourier Series, f_{max} = 32 Hz'))
plot(t_FS_50,real(Q_FSI_FS_50(1:length(t_FS_50))),'DisplayName',strcat('Fourier Series, f_{max} = 50 Hz'))
plot(t_FS_85,real(Q_FSI_FS_85(1:length(t_FS_85))),'DisplayName',strcat('Fourier Series, f_{max} = 85 Hz'))
title('Q Flow from Fourier Series')
xlabel('Time (s)')
ylabel('Flow (m^3/s)')
legend
% ylim([0 180])
hold off

subplot(2,2,4)
hold on
plot(t,Q_FSI,'DisplayName',strcat('FSI'))
plot(t_FS_85,real(Q_FSI_FS_85(1:length(t_FS_85))),'DisplayName',strcat('Fourier Series, f_{max} = 85 Hz'))
title('Q Flow from Fourier Series (Up to Freq = 85 Hz)')
xlabel('Time (s)')
ylabel('Flow (m^3/s)')
% ylim([0 180])
hold off

%% Flow Calculation based on STJ Pressure Fourier Series, and AA Womersley and Ao WK Impedance

for i=1:length(f_shift)

    Z_asc_shift(i) = l_asc*womersley_impedance(r_00, rho, mu, f_shift(i));

    Z_AoWK_shift(i) = R1 + ( R2 / ( 1 + 1i*C*R2*(2*pi*f_shift(i)) ) );

end

Z_total_shift = Z_asc_shift + Z_AoWK_shift;

Q_FT_shift = -P_STJ_FT_shift./Z_total_shift';
Q_FT_shift(f_0) = P_STJ_FT_shift(f_0)/R_total;

[t_Q_FS, Q_FS] = compute_fourier_series(t, Q_FT_shift, f_shift, N, f_0, 85);                 %%f_shift(7034)
Q_FS_real = real(Q_FS);

figure
hold on
plot(t_Q_FS, Q_FS_real,'DisplayName',strcat('d_{AA} = 25.5 mm'))
plot(t_ref, Q_FSI_ref, 'DisplayName',strcat('FSI, d_{AA} = 25.5 mm'))
yline(0)
title('Flow Rate from Fourier Series Using Same Baseline STJ Pressure from FSI')
xlabel('Time (s)')
ylabel('Flow Rate (m^3/s)')
legend('Location','northeast')
grid on
grid minor
hold off

Q_FT_diff = Q_FT_shift - Q_FSI_FT_shift;
Q_FT_shift_fixed = Q_FT_shift - Q_FT_diff;

[t_Q_FS_fixed, Q_FS_fixed] = compute_fourier_series(t, Q_FT_shift_fixed, f_shift, N, f_0, 85);                 %%f_shift(7034)
Q_FS_fixed_real = real(Q_FS_fixed);

figure
subplot(2,1,1)
stem(f_shift(f_0:f_0 + 20), abs(Q_FT_diff(f_0:f_0 + 20)), 'r.', 'MarkerSize', 15)
title('Magnitude Difference between Calculated Q and Q from FSI (Frequency Domain)')
xlabel('Frequency (Hz)')
ylabel('|\DeltaQ|')

subplot(2,1,2)
stem(f_shift(f_0:f_0 + 20), angle(Q_FT_diff(f_0:f_0 + 20)), 'r.', 'MarkerSize', 15)
title('Phase Difference between Calculated Q and Q from FSI (Frequency Domain)')
xlabel('Frequency (Hz)')
ylabel('Angle(\DeltaQ)')

figure
subplot(2,1,2)
hold on
plot(t_Q_FS_fixed, Q_FS_fixed,'DisplayName',strcat('Adjusted, d_{AA} = 25.5 mm'))
plot(t_ref, Q_FSI_ref, 'DisplayName',strcat('FSI, d_{AA} = 25.5 mm'))
yline(0)
title('Adjusted Flow Rate from Fourier Series Using Same Baseline STJ Pressure from FSI')
xlabel('Time (s)')
ylabel('Flow Rate (m^3/s)')
legend('Location','northeast')
grid on
grid minor
hold off

subplot(2,1,1)
hold on
plot(t_Q_FS, Q_FS_real,'DisplayName',strcat('d_{AA} = 25.5 mm'))
plot(t_ref, Q_FSI_ref, 'DisplayName',strcat('FSI, d_{AA} = 25.5 mm'))
yline(0)
title('Flow Rate from Fourier Series Using Same Baseline STJ Pressure from FSI')
xlabel('Time (s)')
ylabel('Flow Rate (m^3/s)')
legend('Location','northeast')
grid on
grid minor
hold off

%% Womersley Impedance Function
function Z_w = womersley_impedance(r, rho, mu, f)
    % Inputs:
    %   r - Vessel radius (m)
    %   rho - Blood density (kg/m^3)
    %   mu - Blood dynamic viscosity (Pa.s)
    %   f - Frequency of pulsatile flow (Hz)
    
    omega = 2 * pi * f;  % Angular frequency (rad/s)
    nu = mu / rho;       % Kinematic viscosity (m^2/s)
    
    % Compute Womersley number
    alpha = r * sqrt(omega / nu);

    % ALPHA = alpha * (1i^(3/2)) ;
    ALPHA = alpha * exp(1i * (3 * pi / 4));
    
    % Womersley function F_10(alpha)
    F_alpha = (2 * besselj(1, ALPHA)) / (ALPHA * besselj(0, ALPHA));
    
    % Calculate longitudinal impedance
    Z_w = ( 1i * omega * rho ) / (1 * pi * (r^2) * (1 - F_alpha));
end

%% Compute Fourier Series Function

function [t_FS, X_FS] = compute_fourier_series(t, X_FT_shift, f_FT_shift, N, f_0, f_max)

t_idx = 0;
dt = 1e-4;

% Pre-allocate arrays for efficiency
num_points = length(0:dt:t(end));
t_FS = zeros(1, num_points);
X_FS = zeros(1, num_points);

for time=0:dt:t(end)

    t_idx =  t_idx + 1;
    t_FS(t_idx) = time;

    % Using exp:
    % X_FS(t_idx) = (X_FT_shift(f_0)/N)*exp(1i*2*pi*f_FT_shift(f_0)*time);

    % Real part represents cosine, imaginary part represents sine
    % Central frequency component (DC)
    X_FS(t_idx) = (real(X_FT_shift(f_0)) / N) * cos(2 * pi * f_FT_shift(f_0) * time) - ...
                  (imag(X_FT_shift(f_0)) / N) * sin(2 * pi * f_FT_shift(f_0) * time);

    i = 1;
    while f_FT_shift(f_0 + i) <= f_max              %fs/2       

        % Using exp:
        % X_FS(t_idx) = X_FS(t_idx) + (X_FT_shift(f_0 + i)/N)*exp(1i*2*pi*f_FT_shift(f_0 + i)*time) + (X_FT_shift(f_0 - i)/N)*exp(1i*2*pi*f_FT_shift(f_0 - i)*time);
       
        % Using sine and cosine:
        % Positive frequency component
            X_FS(t_idx) = X_FS(t_idx) + ...
                (real(X_FT_shift(f_0 + i)) / N) * cos(2 * pi * f_FT_shift(f_0 + i) * time) - ...
                (imag(X_FT_shift(f_0 + i)) / N) * sin(2 * pi * f_FT_shift(f_0 + i) * time);

        % Negative frequency component
            X_FS(t_idx) = X_FS(t_idx) + ...
                (real(X_FT_shift(f_0 - i)) / N) * cos(2 * pi * f_FT_shift(f_0 - i) * time) - ...
                (imag(X_FT_shift(f_0 - i)) / N) * sin(2 * pi * f_FT_shift(f_0 - i) * time);
        
        i = i + 1;
        if (f_0 + i) == length(f_FT_shift)
            break
        end
    end
end

% P_STJ_FS_mmHg = 0.00750062*P_STJ_FS;

end

.

[Hussam]

Thank you very much for your reply.

So FSI, stands for Fluid-Structure Interaction simulation. So basically, the plot for FSI is my reference plot data from my simulations that I am trying to achieve using analytical mathematical modelling.

This leads to STJ, which is the sinotubular junction, which is just above the aortic valve. So, when you see P_STJ, then that is the pressure at the STJ.

So, basically, my FSI simulations include the aortic valve, ascending aorta, and then a 3-element Windkessel model. I can see changes in my FSI simulations when I change geometric parameters of the ascending aorta. We are trying to model analytically the region beyond the valve to see if these changes are due to that region downstream of the valve, or is there more happening at the valve/because of the valve.

In order to model the region downstream of the valve, we have introduced the Womersley model connected in series to a 3-element Windkessel model (1 proximal resistance connected in series to another (distal) resistance and capacitor in parallel).

End result, is that the Womersley + Windkessel model should not give that sine wave at the end of the flow curve - it should approach and stay at zero. Not sure where mathematically I have went wrong.

Regarding the use of fft() - I am using it for the transform, but I was not aware of the time that I can use ifft() to do the inverse. Regardless, I tested ifft() and it gives the same results.

Basically, the idea of calculating flow from pressure is: Q = P/Z in the frequency domain, where Z is the combined impedance of the Womersley and Windkessel models. Z at 0 frequency though is calculated for the Womersley using Poiseuille flow, and for Windkessel is just the summation of the (real) resistance components. Then the results are converted back to the time domain, and Q takes the real part of that.

I hope that answers your questions. Please let me know if you have any further questions.

[Star Strider]

My pleasure!

It’s very difficult for me to follow your code.

If you’re using the Fourier transform to filter your time-domain signal (selecting a specific range of frequencies and then inverting that result), that could account for the oscillation. I note that the last argument (‘fmax’) in your ‘compute_fourier_series’ calls is not always the same, so I suspect this could be the problem. Using ‘boxcar’ filters in this approach frequently causes oscillations (‘Gibbs phenomenon’) near the transitiion.

A significantly better filteriing approach would use the Signal Proceessing Toolbox filter functions. Filter transients might still occur, however they’re much easier to deal with in that instance.

Answer 2

[William Rose]

@Hussam,

[Edit: fix spelling errors.]

I have not tried to understand all of your code. There are many plots above, and I cannot always tell what is what, due to overlapping titles, legends that cover the data, etc.

Are you using measured P_stj to predict (forward) aortic flow at the level of the stj? Are you comparing measured flow to predicted flow, and you have a significant mismatch in diastole? Or are you concerned about a mismatch betweem FSI flow and flow predicted by passing the FSI pressure through a Womersley+Windkessel impedance?

Does the figure below (taken form the post above) illustrate your most significant concern?

Is your principal concern the mismatch between blue flow and FSI flow in late diastole? Where is the "data1" trace in the plot above? Is it important?

The flow magnitudes above are not realistic. Consider FSI systolic flow (orange in the plot above), from t=0 to t=0.2 s. This is period is systolic ejection. Therefore the volume moved in this period is the stroke volume. The mean flow indicated during this period is about 2.5 m^3/s. The stroke volume is the mean flow during ejection times the duration of ejection. Using the mean flow in the plot, this is 2.5 m^3/s x (0.2 s) = 0.5 m^3 = 500 L. That is too much by a factor of about 10^4.

The non-zero blue flow in late diastole (t=0.5 to 0.84) in the plot above may be due to the inability of the Womersley+Windkessel model to account for the nearby aortic valve.

Waves are reflected at branches and bifurcations. Are you accounting for wave reflection where the Womersley tube is terminated at the Windkessel?

Are the units MKS, for the resistance and capacitance values?

You say in your comment "Q takes the real part of that". If the input pressure and flow are real, and one does the impedance calculations correctly in the frequency domain, then, after doing the inverse FFT, the result should be real. The imaginary part should be vanishingly small. If this is not the case, then there is a mistake somewhere. Therefore I recommend that you check to make sure the imaginary part is very small compared to the real part, at all times. If it is not, then check for errors.

I have used Womersley and imedance models to estimate pressures and flows in the circulation.

Johnson, D.A., Rose, W.C., Edwards, J.W., Naik, U.P., Beris, A.N., Application of 1D blood flow models of the human arterial network to differential pressure predictions. J. Biomech. 44: 869-876, 2011.

Johnson, D.A., Spaeth, J.R., Rose, W.C., Edwards, J.W., Naik, U.P., Beris, A.N., An impedance model for blood flow in the human arterial system. Part I: Model development and MATLAB implementation. Computers and Chemical Engineering 35: 1304-1316, 2011.

[Hussam]

Hi William,

Thank you very much for your detailed answer.

Here are my replies to your questions:

Are you using measured P_stj to predict (forward) aortic flow at the level of the stj?

Yes, just above the aortic valve.

Are you comparing measured flow to predicted flow, and you have a significant mismatch in diastole? Or are you concerned about a mismatch betweem FSI flow and flow predicted by passing the FSI pressure through a Womersley+Windkessel impedance?

2nd statement.

Does the figure below (taken form the post above) illustrate your most significant concern?

Yes.

Is your principal concern the mismatch between blue flow and FSI flow in late diastole? Where is the "data1" trace in the plot above? Is it important?

Yes, correct. Also, "data1" is just the horizontal axis (used for visualization) - not important.

The flow magnitudes above are not realistic. Consider FSI systolic flow (orange in the plot above), from t=0 to t=0.2 s. This is period is systolic ejection. Therefore the volume moved in this period is the stroke volume. The mean flow indicated during this period is about 2.5 m^3/s. The stroke volume is the mean flow during ejection times the duration of ejection. Using the mean flow in the plot, this is 2.5 m^3/s x (0.2 s) = 0.5 m^3 = 500 L. That is too much by a factor of about 10^4.

Behind the title you can see that the vertical axis is multiplied by 10^-4. This aligns with what you have stated - to be in the expected magnitude range.

The non-zero blue flow in late diastole (t=0.5 to 0.84) in the plot above may be due to the inability of the Womersley+Windkessel model to account for the nearby aortic valve.

Can you recommend a method to overcome this issue? Or is there a way to include the valve in my calculations?

Waves are reflected at branches and bifurcations. Are you accounting for wave reflection where the Womersley tube is terminated at the Windkessel?

No, how can I do that?

Are the units MKS, for the resistance and capacitance values?

Yes.

You say in your comment "Q takes the real part of that". If the input pressure and flow are real, and one does the impedance calculations correctly in the frequency domain, then, after doing the inverse FFT, the result should be real. The imaginary part should be vanishingly small. If this is not the case, then there is a mistake somewhere. Therefore I recommend that you check to make sure the imaginary part is very small compared to the real part, at all times. If it is not, then check for errors.

Verified that the imaginary part of the calculated flow is very small (negligible).

Thank you for sharing your papers as well. I will read them. I hope these answer your questions. Looking forward to your reply.

Many thanks,

Hussam

[William Rose]

@Hussam,

This an interesting discussion which I woudl like to continue. Please send me a secure message by clicking on the "WR" next to my posts, and then click on the envelope icon. I will reply.

[Hussam]

Sent :)