%This function calculates a fully developed velocity profile for a straight %tube according to Womersley's theory. It is based on knowing the mean %velocity. %For full details see "Ultrasound Imaging Velocimetry with interleaved %images for improved pulsatile arterial flow measurements: a new correction %method, experimental and in vivo validation", J Royal Soc Interface, 2017 %K H Fraser, C Poelma 2016 function vel= womer_Vmean(t, freqs, V_mean, R, nu) %t = time %freqs = number of harmonics to use %V_peak = peak velocity waveform (same number of points as t) %R = radius of artery [p1,q1]=size(t); tot_size=p1*q1; q=length(t); p=tot_size/q; %period time is just the last time in the time vector T=t(end); % main execution loop. calls the function d2 which sets up graphs and calls % function d1 which actually calculates and sums Womersley profile harmonics. n=1; for r=1:q velprof = d2(t(r),freqs, V_mean, R, T, nu); vel(r,:) = velprof; %whole profile end function velprof = d2(t,freqs, V_mean, R, T,nu) % the fourier coefficients. Dividing by the number of samples is necessary % here as the fourier coefficients are calculated with a factor of N. B=fft(V_mean)/length(V_mean); % sets up a complex or real grid, depending on what you are plotting %z = cplxgrid(20); z=0:0.01:1; velprof = 2*B(1)*(1-abs(z).^2) + d1(z,t,B,R, freqs, T,nu); % this function actually calculates the Womersley harmonics. They are summed % in the variable k function k = d1(z,t,B,R, freqs, T,nu) % frequency -in rad/s w=2*pi/T; % Womersley number for mode 0 A0=R*sqrt(w/nu); k=0; % this sums the various components to form the complete Womersley profile for y =2:freqs An=A0*sqrt(y-1); top = besselj(0,An*1i^(3/2))-besselj(0,An*abs(z)*1i^(3/2)); % bottom = besselj(0,An*1i^(3/2)) - 2*besselj(1,An*1i^(3/2))/(An*i^(3/2));% in these 2 lines Holdsworth paper has i^-(3/2) but now fairly sure it should be + k = k+2*B(y)*(top/bottom)*exp(w*i*(y-1)*t); end