From 2ebea8cffbd05dd4bd71856a53e5ce355acea366 Mon Sep 17 00:00:00 2001 From: Mark Palmeri Date: Mon, 16 Apr 2018 20:41:29 -0400 Subject: [PATCH 1/4] integrate Ned's latest processing updates * All of the hard coded analysis and plotting variables have been moved to the top level wrapper function and are saved in a `par` parameters structure. * The code has been modified to write the plot data to a CSV file. This file can be imported into a spreadsheet if desired. It might be good to add the software version to this file. * I have added two path directories for the saving of (1) intermediate analysis files and (2) the plot and CSV file. These directories are included in the par structure in the wrapper file. For now, the default paths are `pwd`. The code still looks in the current directory to find the raw data, though it would also be easy to change this location. * The maximum frequency included in the plots and CSV output is now determined automatically as the frequency when the phase velocity 95% CI exceeds the 30% of average range. * The last function that makes the plot and saves the text CSV file has been renamed as `MakePlotAndSaveOutputCSVfile.m`. * The last function also generates an output `errorFlag` which is returned as `0` if the results are judged to be valid, and `1` if the results are not valid. Currently, the results are reported as not valid if (1) the gSWS for either displacement or velocity are outside the range 0.5 - 5 m/s (these limits are set in the `par` structure in the wrapper) or (2) if the gSWS 95% CI exceeds the 30% of average for either the displacement or velocity results. --- processing_code/AnalyzeARFIdata.m | 16 +- processing_code/AnalyzeAllAcquisitions.m | 50 +++- .../MakePlotAndSaveOutputCSVfile.m | 218 ++++++++++++++++++ 3 files changed, 264 insertions(+), 20 deletions(-) create mode 100755 processing_code/MakePlotAndSaveOutputCSVfile.m 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 100755 index 0000000..b87064b --- /dev/null +++ b/processing_code/MakePlotAndSaveOutputCSVfile.m @@ -0,0 +1,218 @@ + +%========================================================================== + +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 = nanmean(Vdisp); % find avg,std gSWS results + Vdisp_std = nanstd(Vdisp); + Vvel_avg = nanmean(Vvel); + Vvel_std = nanstd(Vvel); + + [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)') +% title(phantomID) + 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); + if length(fbadidx)>0 + fmax=freqsToAnalyzeHz(min(fbadidx)-1); + end + + fidx = find(kr>=krThreshold & freqsToAnalyzeHz<=fmax); % get plot fidx + ddf=2; + df=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))-df/2; + patchx(2*(nn+1)-ii+1) = freqsToAnalyzeHz(fidx(ii))-df/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))+df/2; + patchx(nn+2) = freqsToAnalyzeHz(fidx(ii-1))+df/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))-df/2; + patchx(2*(nn+1)-ii+1) = freqsToAnalyzeHz(fidx(ii))-df/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))-ddf freqsToAnalyzeHz(fidx(ii))+ddf],[ybot ybot], 'Color',black,'LineWidth',2) + plot([freqsToAnalyzeHz(fidx(ii))-ddf freqsToAnalyzeHz(fidx(ii))+ddf],[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 & temp0 + 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 + +%========================================================================== From 3204912c4caa97db42fdd633b01ae014e55b4196 Mon Sep 17 00:00:00 2001 From: Derek Chan <31013934+derekychan@users.noreply.github.com> Date: Wed, 18 Apr 2018 11:37:35 -0400 Subject: [PATCH 2/4] include kr threshold in determination of max frequency (#12) (#13) --- processing_code/MakePlotAndSaveOutputCSVfile.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing_code/MakePlotAndSaveOutputCSVfile.m b/processing_code/MakePlotAndSaveOutputCSVfile.m index b87064b..2edaec3 100755 --- a/processing_code/MakePlotAndSaveOutputCSVfile.m +++ b/processing_code/MakePlotAndSaveOutputCSVfile.m @@ -96,7 +96,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); + fbadidx = find(CIfactor*phVel_std>fracErrorBar*phVel_avg & kr>=krThreshold); if length(fbadidx)>0 fmax=freqsToAnalyzeHz(min(fbadidx)-1); end From ec7068e61b47526db8727ca2a392e39c95449917 Mon Sep 17 00:00:00 2001 From: Derek Chan <31013934+derekychan@users.noreply.github.com> Date: Thu, 18 Oct 2018 14:24:17 -0400 Subject: [PATCH 3/4] updates to processing code (#14) * omit bad points from group SWS calculation * Remove commented lines * Fix code syntax formatting for `MakePlotAndSaveOutputCSVfile.m` * Move computation of `fidx` to separate function * Rename plotting variables for clarity --- .../MakePlotAndSaveOutputCSVfile.m | 237 ++++++++++-------- 1 file changed, 132 insertions(+), 105 deletions(-) mode change 100755 => 100644 processing_code/MakePlotAndSaveOutputCSVfile.m diff --git a/processing_code/MakePlotAndSaveOutputCSVfile.m b/processing_code/MakePlotAndSaveOutputCSVfile.m old mode 100755 new mode 100644 index 2edaec3..197e6bd --- a/processing_code/MakePlotAndSaveOutputCSVfile.m +++ b/processing_code/MakePlotAndSaveOutputCSVfile.m @@ -23,132 +23,128 @@ 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)); + 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'); + S = load(analysisFiles(ifile).name, 'Vdisp', 'Vvel', 'freqsToAnalyzeHz', 'phVel'); Vdisp(ifile) = S.Vdisp; Vvel(ifile) = S.Vvel; - phVel(ifile,:) = S.phVel; + phVel(ifile, :) = S.phVel; freqsToAnalyzeHz = S.freqsToAnalyzeHz; clear S end - Vdisp_avg = nanmean(Vdisp); % find avg,std gSWS results - Vdisp_std = nanstd(Vdisp); - Vvel_avg = nanmean(Vvel); - Vvel_std = nanstd(Vvel); - - [phVel_avg,phVel_std] = MeanStdOmitBadPtsOneCase_usingPositiveData(phVel,thFactor,maxSpeed); - + [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; + x1 = 1.0; + x2 = 2.0; dx = 0.04; - YTickVals = [0:1:maxPlotSpeed]; + YTickVals = 0:1:maxPlotSpeed; figure(1) clf - set(gcf,'Units','inches','Position',[0.5 2 8 3.5]) + set(gcf, 'Units', 'inches', 'Position', [0.5 2 8 3.5]) - hhh1=subplot(1,2,1); + 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) + 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') + 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)') -% title(phantomID) - set(gca,'Units','inches','Position',[0.5 0.5 1.4 2.7]) + 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 + 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 + 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); + if length(fbadidx) > 0 + fmax = freqsToAnalyzeHz(min(fbadidx) - 1); end - fidx = find(kr>=krThreshold & freqsToAnalyzeHz<=fmax); % get plot fidx - ddf=2; - df=10; - nn=length(fidx); - patchx = NaN(1,2*nn+2); - patchy = NaN(1,2*nn+2); + 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 + for ii = 1:nn+1 switch ii case 1 - patchx(ii) = freqsToAnalyzeHz(fidx(ii))-df/2; - patchx(2*(nn+1)-ii+1) = freqsToAnalyzeHz(fidx(ii))-df/2; - patchy(ii) = (1+fracErrorBar)*phVel_avg(fidx(ii)); - patchy(2*(nn+1)-ii+1) = (1-fracErrorBar)*phVel_avg(fidx(ii)); + 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))+df/2; - patchx(nn+2) = freqsToAnalyzeHz(fidx(ii-1))+df/2; - patchy(ii) = (1+fracErrorBar)*phVel_avg(fidx(ii-1)); - patchy(nn+2) = (1-fracErrorBar)*phVel_avg(fidx(ii-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))-df/2; - patchx(2*(nn+1)-ii+1) = freqsToAnalyzeHz(fidx(ii))-df/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))); + 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); + hhh2 = subplot(1, 2, 2); hold on - fill(patchx,patchy,[0.8 0.8 0.8],'EdgeColor','none'); + fill(patchx, patchy, [0.8 0.8 0.8], 'EdgeColor', 'none'); - for ii=1:length(fidx) + 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))-ddf freqsToAnalyzeHz(fidx(ii))+ddf],[ybot ybot], 'Color',black,'LineWidth',2) - plot([freqsToAnalyzeHz(fidx(ii))-ddf freqsToAnalyzeHz(fidx(ii))+ddf],[ytop ytop], 'Color',black,'LineWidth',2) + 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') + 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]) + set(gca, 'Units', 'inches', 'Position', [2.4 0.5 5.4 2.7]) savefig([saveDir '\gSWS_phVel_figure']) @@ -157,23 +153,23 @@ 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),... + 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),... + 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'); + 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))),... + 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))) ); @@ -181,37 +177,68 @@ 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 ) + 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) +function [out_avg, out_std] = MeanStdOmitBadPtsOneCase_usingPositiveData(data, factor, maxSpeed) - [npts,nfreqs] = size(data); - out_avg = NaN(1,nfreqs); - out_std = NaN(1,nfreqs); + [npts, nfreqs] = size(data); + out_avg = NaN(1, nfreqs); + out_std = NaN(1, nfreqs); - for ifreq=1:nfreqs - temp = data(:,ifreq); + for ifreq = 1:nfreqs + temp = data(:, ifreq); idx = find(temp>0 & temp0 - goodidx = find( abs(temp-temp_avg)0); + if temp_std > 0 + goodidx = find(abs(temp-temp_avg)0); else - goodidx=1:npts; + goodidx = 1:npts; end out_avg(ifreq) = nanmean(temp(goodidx)); - out_std(ifreq) = nanstd(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 From a7871632b3158ffb5b21e096fc87fae97c8877a8 Mon Sep 17 00:00:00 2001 From: Mark Palmeri Date: Thu, 18 Oct 2018 14:28:15 -0400 Subject: [PATCH 4/4] update CHANGELOG --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index adbb508..40c0e9c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ * `routines/save_IQ_data.m` written to save data files that can be post-processed by user. -# v1.3a +# v1.3 * 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