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.