diff --git a/CHANGELOG.md b/CHANGELOG.md index 0ba0d7f..eab6add 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,3 +8,8 @@ * Change tracking time from 20 ms -> 40 ms for more robust dispersion analysis * Change from Apache v2.0 to MIT license * Fix calculation of C5-2 B-mode aperture (Issue #16) + +# v1.4.0 +* change tracking time from 20 ms -> 40 ms for more robust dispersion analysis +* change from Apache v2.0 to MIT license +* omit bad points from group SWS calculation diff --git a/processing_code/AnalyzeARFIdata.m b/processing_code/AnalyzeARFIdata.m index 784d8f0..36d6b50 100755 --- a/processing_code/AnalyzeARFIdata.m +++ b/processing_code/AnalyzeARFIdata.m @@ -1,21 +1,10 @@ %========================================================================== -function AnalyzeARFIdata(filestamp) +function AnalyzeARFIdata(filestamp,par) % - par.DOFmm = 2; % parameters for processing arfidata - par.minLatMM = 4; - par.maxLatMM = 15; - par.maxTimeMS = 40; - par.minTimeMS = -15; - par.rolloffTimeMS = 15; - par.nStepsToRemove = 2; - par.LPFcutoffKHz = 1; - par.desirecPRFkHz = 20; - par.freqsToAnalyzeHz = 20:10:800; - filename = [filestamp '_fromIQ_arfidata.mat']; S = load(filename,'arfidata','axial','lat','numrefs','push_focus','T'); arfidata = double(S.arfidata); % arfidata saved in um @@ -46,9 +35,8 @@ function AnalyzeARFIdata(filestamp) freqsToAnalyzeHz = par.freqsToAnalyzeHz; phVel = Construct2DFTandCalcPhVel(wtVelPlane,t,lat,freqsToAnalyzeHz); - saveFile = [filestamp '_FL' num2str(fliplat) '_phVel_gSWS_data.mat']; + saveFile = [par.analysisDir '\' filestamp '_FL' num2str(fliplat) '_phVel_gSWS_data.mat']; save(saveFile,'dispPlane','tms','latmm','velPlane','veltms','vellatmm','TTPdisp','TTPvel','Vdisp','Vvel','freqsToAnalyzeHz','phVel') - end end diff --git a/processing_code/AnalyzeAllAcquisitions.m b/processing_code/AnalyzeAllAcquisitions.m index 0aa31d1..3b8c3de 100755 --- a/processing_code/AnalyzeAllAcquisitions.m +++ b/processing_code/AnalyzeAllAcquisitions.m @@ -19,22 +19,60 @@ % AnalyzeARFIdata % SetupARIFdataPlane % Construct2DFTandCalcPhVel -% MakePlot +% MakePlotAndSaveCSVfile - [filepath, name, ext] = fileparts(mfilename('fullpath')); - addpath(filepath) + + + % analysis and save directories + + par.analysisDir = pwd; % directory for saving analysis data + par.saveDir = pwd; % directory to save plot and CSV file + + % parameters for processing arfidata + + par.DOFmm = 2; % depth of field to average (mm) + par.minLatMM = 4; % min lateral position (mm) to start analysis + par.maxLatMM = 15; % max lateral position (mm) to end analysis + par.maxTimeMS = 40; % maximum time (ms) to use in analysis + par.minTimeMS = -15; % minimum time (ms) to zero-pad before push + par.rolloffTimeMS = 15; % time (ms) for rolloff at early and late times + par.nStepsToRemove = 2; % number of reverb steps to remove + par.LPFcutoffKHz = 1; % low-pass filter cutoff frequency (kHz) + par.desirecPRFkHz = 20; % desired PRF (kHz) after upsampling in time + par.freqsToAnalyzeHz = 20:10:800; % frequencies (Hz) for phase velocity measurements + + % parameters for plotting and saving output - paramFiles = dir('*_parameters.mat'); % get parameter files in directions + par.phantomID = 'E2297-A3'; % phantom ID string + par.maxPlotSpeed = 4.5; % maximum speed for gSWS and phVel plots + par.maxplotfreq = 820; % max frequency (Hz) for plot + par.fracErrorBar = 0.3; % fixed fraction for error bars + par.CIfactor = 1.96; % confidence interval for error bars + par.fmax = 800; % maximum frequency (Hz) to plot + par.rmin = 0.004; % minimum lateral position (m) for low freq cutoff + par.krThreshold = 1.5; % kr threshold for low freq cutoff + par.maxSpeed = 8; % max speed for "good" result + par.thFactor = 2; % phVel accept range = +/- thfactor * std + par.maxValidGSWS = 5; % maximum group speed where results are valid + par.minValidGSWS = 0.5; % minimum group speed where results are valid + + +% [filepath, name, ext] = fileparts(mfilename('fullpath')); +% addpath(filepath) + + paramFiles = dir('*_parameters.mat'); % get parameter files for i = 1:length(paramFiles) % for each file, get filestamp and compute displacements disp(['Processing ' num2str(i) ' of ' num2str(length(paramFiles))]) % comment to suppress output filestamp = paramFiles(i).name(1:14); genDispMTL(filestamp) - AnalyzeARFIdata(filestamp); + AnalyzeARFIdata(filestamp,par); end - MakePlot + + errorFlag = MakePlotAndSaveOutputCSVfile(par); + end %========================================================================== diff --git a/processing_code/MakePlotAndSaveOutputCSVfile.m b/processing_code/MakePlotAndSaveOutputCSVfile.m new file mode 100644 index 0000000..197e6bd --- /dev/null +++ b/processing_code/MakePlotAndSaveOutputCSVfile.m @@ -0,0 +1,245 @@ + +%========================================================================== + +function errorFlag = MakePlotAndSaveOutputCSVfile(par) + + % + + analysisDir = par.analysisDir; % directory with result files saved in analysis + saveDir = par.saveDir; % directory to save output plot and CSV data + + freqsToAnalyzeHz = par.freqsToAnalyzeHz; % analysis frequencies (Hz) + phantomID = par.phantomID; + maxPlotSpeed = par.maxPlotSpeed; + maxplotfreq = par.maxplotfreq; + fracErrorBar = par.fracErrorBar; + CIfactor = par.CIfactor; + fmax = par.fmax; + rmin = par.rmin; + krThreshold = par.krThreshold; + maxSpeed = par.maxSpeed; + thFactor = par.thFactor; + + + analysisFiles = dir([analysisDir '\*_FL*_phVel_gSWS_data.mat']); % get analysis files + + Vdisp = NaN(1, length(analysisFiles)); + Vvel = NaN(1, length(analysisFiles)); + phVel = NaN(2 * length(analysisFiles), length(freqsToAnalyzeHz)); + + for ifile = 1:length(analysisFiles) % load Vdisp, Vvel, phVel data + + S = load(analysisFiles(ifile).name, 'Vdisp', 'Vvel', 'freqsToAnalyzeHz', 'phVel'); + Vdisp(ifile) = S.Vdisp; + Vvel(ifile) = S.Vvel; + phVel(ifile, :) = S.phVel; + freqsToAnalyzeHz = S.freqsToAnalyzeHz; + clear S + end + + [Vdisp_avg, Vdisp_std] = MeanStdOmitBadPtsOneCase_gSWS_usingPositiveData(Vdisp, thFactor, maxSpeed); + [Vvel_avg, Vvel_std] = MeanStdOmitBadPtsOneCase_gSWS_usingPositiveData(Vvel, thFactor, maxSpeed); + + [phVel_avg, phVel_std] = MeanStdOmitBadPtsOneCase_usingPositiveData(phVel, thFactor, maxSpeed); + + black = [0.0 0.0 0.0]; % make plot + gray = [0.4 0.4 0.4]; + x1 = 1.0; + x2 = 2.0; + dx = 0.04; + YTickVals = 0:1:maxPlotSpeed; + + figure(1) + clf + set(gcf, 'Units', 'inches', 'Position', [0.5 2 8 3.5]) + + hhh1 = subplot(1, 2, 1); + hold on + ybot = (1 - fracErrorBar) * Vdisp_avg; % +/- fraction of avg displacement gSWS + ytop = (1 + fracErrorBar) * Vdisp_avg; + plot([x1 x1], [ybot ytop], 'Color',gray, 'LineWidth',1) + plot([x1-dx x1+dx], [ybot ybot], 'Color',gray, 'LineWidth',1) + plot([x1-dx x1+dx], [ytop ytop], 'Color',gray, 'LineWidth',1) + + ybot = Vdisp_avg - CIfactor * Vdisp_std; % +/- CI + ytop = Vdisp_avg + CIfactor * Vdisp_std; + plot([x1 x1], [ybot ytop], 'Color',black,'LineWidth',2) + plot([x1-dx x1+dx], [ybot ybot], 'Color',black,'LineWidth',2) + plot([x1-dx x1+dx], [ytop ytop], 'Color',black,'LineWidth',2) + + ybot = (1 - fracErrorBar) * Vvel_avg; % +/- 30% of avg velocity gSWS + ytop = (1 + fracErrorBar) * Vvel_avg; + plot([x2 x2], [ybot ytop], 'Color',gray, 'LineWidth',1) + plot([x2-dx x2+dx], [ybot ybot], 'Color',gray, 'LineWidth',1) + plot([x2-dx x2+dx], [ytop ytop], 'Color',gray, 'LineWidth',1) + + ybot = Vvel_avg - CIfactor * Vvel_std; % +/- CI + ytop = Vvel_avg + CIfactor * Vvel_std; + plot([x2 x2], [ybot ytop], 'Color',black,'LineWidth',2) + plot([x2-dx x2+dx], [ybot ybot], 'Color',black,'LineWidth',2) + plot([x2-dx x2+dx], [ytop ytop], 'Color',black,'LineWidth',2) + hold off + + set(gca, 'Units', 'inches', 'YTick', YTickVals) + set(gca, 'XTick', [1 2], 'XTickLabel', {'disp' 'vel'}) + set(gca, 'Box', 'on') + ylim([0 maxPlotSpeed]) + xlim([0 3]) + ylabel('group SWS (m/s)') + set(gca, 'Units', 'inches', 'Position', [0.5 0.5 1.4 2.7]) + + kr = 2*pi*freqsToAnalyzeHz ./ phVel_avg .* rmin; % kr values to find min freq in plot + + fmax = max(freqsToAnalyzeHz); % max freq based on size of error bars + fbadidx = find(CIfactor*phVel_std>fracErrorBar*phVel_avg & kr>=krThreshold); + if length(fbadidx) > 0 + fmax = freqsToAnalyzeHz(min(fbadidx) - 1); + end + + fidx = GetPlotFreqIdx(kr, krThreshold, freqsToAnalyzeHz, fmax); + + capsize = 2; % size of caps on black error bars + patchwidth = 10; + nn = length(fidx); + patchx = NaN(1, 2*nn+2); + patchy = NaN(1, 2*nn+2); + + for ii = 1:nn+1 + switch ii + case 1 + patchx(ii) = freqsToAnalyzeHz(fidx(ii)) - patchwidth/2; + patchx(2*(nn+1)-ii+1) = freqsToAnalyzeHz(fidx(ii)) - patchwidth/2; + patchy(ii) = (1 + fracErrorBar) * phVel_avg(fidx(ii)); + patchy(2*(nn+1)-ii+1) = (1 - fracErrorBar) * phVel_avg(fidx(ii)); + case nn+1 + patchx(ii) = freqsToAnalyzeHz(fidx(ii-1)) + patchwidth/2; + patchx(nn+2) = freqsToAnalyzeHz(fidx(ii-1)) + patchwidth/2; + patchy(ii) = (1 + fracErrorBar) * phVel_avg(fidx(ii-1)); + patchy(nn+2) = (1 - fracErrorBar) * phVel_avg(fidx(ii-1)); + otherwise + patchx(ii) = freqsToAnalyzeHz(fidx(ii)) - patchwidth/2; + patchx(2*(nn+1)-ii+1) = freqsToAnalyzeHz(fidx(ii)) - patchwidth/2; + patchy(ii) = (1 + fracErrorBar) / 2 * (phVel_avg(fidx(ii-1)) + phVel_avg(fidx(ii))); + patchy(2*(nn+1)-ii+1) = (1 - fracErrorBar) / 2 * (phVel_avg(fidx(ii-1)) + phVel_avg(fidx(ii))); + end + end + + hhh2 = subplot(1, 2, 2); + hold on + fill(patchx, patchy, [0.8 0.8 0.8], 'EdgeColor', 'none'); + + for ii = 1:length(fidx) + + ybot = phVel_avg(fidx(ii)) - CIfactor * phVel_std(fidx(ii)); % +/- CI + ytop = phVel_avg(fidx(ii)) + CIfactor * phVel_std(fidx(ii)); + plot([freqsToAnalyzeHz(fidx(ii)) freqsToAnalyzeHz(fidx(ii)) ], [ybot ytop], 'Color', black, 'LineWidth', 2) + plot([freqsToAnalyzeHz(fidx(ii))-capsize freqsToAnalyzeHz(fidx(ii))+capsize], [ybot ybot], 'Color', black, 'LineWidth', 2) + plot([freqsToAnalyzeHz(fidx(ii))-capsize freqsToAnalyzeHz(fidx(ii))+capsize], [ytop ytop], 'Color', black, 'LineWidth', 2) + end + hold off + set(gca, 'Units', 'inches', 'YTick', YTickVals) + set(gca, 'Box', 'on') + ylim([0 maxPlotSpeed]) + xlim([0 maxplotfreq]) + ylabel('phase velocity (m/s)') + xlabel('frequency (Hz)') + title(phantomID) + set(gca, 'Units', 'inches', 'Position', [2.4 0.5 5.4 2.7]) + + savefig([saveDir '\gSWS_phVel_figure']) + + + % make and save CSV file + + saveFile = [saveDir '\gSWS_phVel_data.txt']; + + fid = fopen(saveFile, 'w'); + fprintf(fid, '\n'); + fprintf(fid, 'phantom %s\n', phantomID); + fprintf(fid, '\n'); + fprintf(fid, ',group SWS (m/s), 95%% CI, 30%% mean\n'); + fprintf(fid, 'displacement,%s,%s,%s\n', sprintf('%7.3f', Vdisp_avg),... + sprintf('%7.3f', 1.95*Vdisp_std),... + sprintf('%7.3f', 0.3*Vdisp_avg) ); + fprintf(fid, 'velocity,%s,%s,%s\n', sprintf('%7.3f', Vvel_avg),... + sprintf('%7.3f', 1.95*Vvel_std),... + sprintf('%7.3f', 0.3*Vvel_avg) ); + fprintf(fid, '\n'); + fprintf(fid, '\n'); + fprintf(fid, 'frequency (Hz),phase velocity (m/s), 95%% CI, 30%% mean\n'); + + for ii = 1:length(fidx) + fprintf(fid, '%s,%s,%s,%s\n', sprintf('%7.0f', freqsToAnalyzeHz(fidx(ii))),... + sprintf('%7.3f', phVel_avg(fidx(ii))),... + sprintf('%7.3f', 1.96*phVel_std(fidx(ii))),... + sprintf('%7.3f', 0.30*phVel_avg(fidx(ii))) ); + end + fclose(fid); + + errorFlag = 0; + if (Vdisp_std*CIfactor > fracErrorBar*Vdisp_avg || Vvel_std*CIfactor > fracErrorBar*Vvel_avg ... + || Vdisp_avg > par.maxValidGSWS || Vdisp_avg < par.minValidGSWS ... + || Vvel_avg > par.maxValidGSWS || Vvel_avg < par.minValidGSWS) + errorFlag=1; + end +end + +%========================================================================== + +function [out_avg, out_std] = MeanStdOmitBadPtsOneCase_usingPositiveData(data, factor, maxSpeed) + + [npts, nfreqs] = size(data); + out_avg = NaN(1, nfreqs); + out_std = NaN(1, nfreqs); + + for ifreq = 1:nfreqs + temp = data(:, ifreq); + idx = find(temp>0 & temp 0 + goodidx = find(abs(temp-temp_avg)0); + else + goodidx = 1:npts; + end + out_avg(ifreq) = nanmean(temp(goodidx)); + out_std(ifreq) = nanstd(temp(goodidx)); + end +end + +%========================================================================== + +function [gSWS_avg, gSWS_std] = MeanStdOmitBadPtsOneCase_gSWS_usingPositiveData(data, factor, maxSpeed) + + idx = find(data>0 & data 0 + goodidx = find(abs(data-temp_avg)0); + else + goodidx = 1:length(data); + end + gSWS_avg = nanmean(data(goodidx)); + gSWS_std = nanstd(data(goodidx)); +end + +%========================================================================== + +function freqIdx = GetPlotFreqIdx(kr, krThreshold, freqsToAnalyzeHz, fmax) + + freqIdx = find(kr>=krThreshold & freqsToAnalyzeHz<=fmax); + if length(freqIdx) == 0 + fmax = max(freqsToAnalyzeHz); + freqIdx = find(kr>=krThreshold & freqsToAnalyzeHz<=fmax); + end +end + +%==========================================================================