% JOVE - Blood flow imaging with Ultrafast Doppler
% Jerome Baranger1, Luc Mertens, Olivier Villemain
%
% SUPPLEMENTARY MATERIAL
%
% This script presents a post-processing example of the IQ data obtained 
% with any ultrafast ultrasound scanner.

%% 1    - Data filtering : removing the tissue echoes wiht singular value
% decompositiom (See Demene et al, 2015)
% 1.1   - Reshape IQ matrix
[nz,nx,nt]              = size(IQ); 
IQr                     = reshape(IQ,[nz*nx,nt]);
[U,S,V]                 = svd(IQr,'econ');

% 1.2   - SVD threshold determination : compute Spatial Similarity Matrix C
C                       = corr(abs(U));
figure, imagesc(C); colormap jet(256);

N                       = 30;   % >>ENTER<< the boundary N between the blood/tissue subspaces (See Baranger et al,
                                % 2018)

% 1.2  - Filter IQ
If                      = ones(1,nt);
If(1:N)                 = 0;
If                      = diag(If);
IQf                     = U*S*If*V;
IQf                     = reshape(IQf,[nz,nx,nt]);

%% 2    - Compute and display Power Doppler map
PD                      = mean( abs( IQf ).^2 ,3);
PD_dB                   = 10*log10(PD/max(PD(:)));
figure(122), imagesc(PD_dB); colormap hot(256); colorbar; caxis([-30 0]);
axis equal tight;

%% 3    - Display spectrogram for a Region of Interest

% 3.1  - Define ROI
figure(122);
hROI                    = drawcircle;

%% ---> Move the ROI and reexecute this cell to analyse different areas of 
%       the field of view.

% 3.2  - Average filtered IQ on that ROI
roiIdx                  = find( hROI.createMask );
IQf_r                   = reshape(IQf,[nz*nx,nt]);
IQf_ROI                 = mean( IQf_r(roiIdx,:) ,1);

% 3.3  - Compute spectrogram PSD
Fr                      = 4000;     % >>ENTER<< the framerate of the sequence in Hz.
win                     = 30;       % Window size
[~,f,t,PSD]             = spectrogram(IQf_ROI,hann(win),floor(0.9*win),128,Fr,'psd');
f                       = f - Fr/2; % Center the requency range aorund zero;

% 3.4  - Compute center frequency fc
fc                      = meanfreq(PSD,f);

% 3.5  - Convert frequecy values into blood velocities
TW_f                    = 7.3;      % >>ENTER<< transmit waveform frequnecy in MHz
SoundSpeed              = 1480;     % >>ENTER<< the speed of sound in the medium of the sequence in m/s.
velocityRange           = f * SoundSpeed /(2*TW_f*1e6);
bloodVelocity_meanROI   = fc * SoundSpeed /(2*TW_f*1e6);


% 3.6  - Show results
figure(123), imagesc(t,velocityRange,PSD); colormap gray(256);
hold on, plot(t,bloodVelocity_meanROI,'g--');


