-
Notifications
You must be signed in to change notification settings - Fork 0
/
HiFiPDV2.m
2032 lines (1718 loc) · 75.4 KB
/
HiFiPDV2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
%% **************** High Fidelity PDV Analysis Version 2.0 ****************
% ************************************************************************
% ************************************************************************
%
% This code investigates a broad array of reasonable spectrogram inputs
% and peak extraction methods, then analyzes the distribution of output
% histories to determine the best fit & uncertainty.
%
% Author: Travis J. Voorhees, [email: [email protected]]
% Most recent release: October 2, 2023
% Most recent bug fix: May 6, 2024
%
% Original HiFiPDV release was on December 12, 2020, embedded in PhD
% thesis and included in the associated copyright (Travis Voorhees 2020).
% =========================================================================
%% CODE DESCRIPTION & NOTES:
% This code is intended to be run with parallel CPU computing only
%
% Other Versions of this Code That Are Available:
% - GPU computing version (deprecated)
% -- Requires large amounts of GPU VRAM
% -- Latest GPU code update is December 2020
%
% Planned Updates/Upgrades (also see latest release notes):
% - Ability to load in ROI bounds from previous run.
% - Choose your favorite spectrogram to put in final plot.
% =========================================================================
%% INPUT VARIABLES INVESTIGATED:
% - Duration Length (length of timebin)
% -- 5 to 100 ns, in 5 ns steps
% -- Duration Overlap is scaled to produce a constant 2.5 ns time-step
% - Window function
% -- Hamming
% -- Hahn
% -- Blackman
% - Peak determination method
% -- Max of power spectrum at time
% -- Centroid (center of mass) of power spectrum at time
% -- Peak of Gaussian fit to power spectrum at time
%
% Total iterations to be performed: [DurLength]*[WinFun]*[PeakEx]=180
% NumSpec = 20*3;
% NumHist = NumSpec*3; % number of histories per window function
% =========================================================================
%% INPUT VARIABLES THAT ARE HELD CONSTANT
% - Number of frequency bins
% -- 2^15 = 32,768 bins
% -- Results in slighter smaller than 1 m/s velocity bin size for max
% frequency on HSR lab high speed scope
% - Time-step / point-to-point time resolution
% -- ~2.5 ns skip between DFT calculations
% -- Overlap is scaled with duration length accordingly
% -- NOTE: Time-step needs to be a factor of all window lengths, so
% the time-step may expand beyond 2.5ns to accomodate this.
% Code output will tell you if this happens, and explicitly
% state the new time-step.
% =========================================================================
%% ASSUMPTIONS MADE IN THIS CALCULATION
% (1) A single velocity signal is present within the user's ROI bounds
% (2) Signal spectral width does not collide with neighboring signals
% (3) No window corrections are applied (can be added if requested)
% (4) No PDV probe angle corrections are applied (assumed normal to vector)
% =========================================================================
%% LATEST RELEASE NOTES:
% ============================ (May 6, 2024) ===========================
% Bug Fixes:
% - Replaced range function with max-min to avoid mix-up of the standard
% MATLAB range function with the Communications toolbox's range
% function.
% -- Thank you Colton Cagle for finding this bug!
% - Added a logic statement prior to calculating the mean of final
% histories. Versions of MATLAB prior to R2020a use 'nanmean',
% versions R2020b or newer use 'mean' with the added flag of
% 'omitmissing'
% -- Also thank you to Colton Cagle for finding this bug!
%
% ============================ (October 2, 2023) ===========================
% Added:
% - Advanced user options for colormap choice and baseline filtering
% - More descriptive comments/notes
% - More consistent output file naming
%
% ============================ (September 1, 2023) ===========================
% Fixed:
% - Gaussian fitting procedure for Heatmap's history extraction
% -- Original method was throwing spurious results occasionally
%
% ============================ (August 16, 2023) ===========================
% Changed:
% - More efficient Gaussian fitting process -- thank you Dan Champion!
% - Spectrum outside region of interest set to zero (instead of -60 dB)
% -- Helps with new Gaussian fitting process
% - Moved Cropping parameters outside of Spectrogram parfor loop
% -- Makes cropping vectors permanent variables
% - Removed negative velocities from Heatmap calculation
% - Memory reduction via FP32 or FP16 data storage (laptop friendly).
%
% Fixed:
% - Calculation of baseline frequency
% -- Previous version had an typo in conversion to real space for
% Blackman and Hann windows
%
% Beta Testing:
% - Reduced spectrogram bit precision to FP16 instead of FP64
% -- Significant RAM reduction, negligible precision reduction
% - Faster Gaussian fitting program
% - Better cropping of signal during peak fitting
%
% Remaining Updates Before Full Version Release:
% - Updates for down-shifted signals
% -- Not enough urgency to add this functionality
%
% ============================ (March 6, 2023) ============================
% Added:
% - Ability to read in HDF5 (*.h5) files
% -- Borrows "agilend_h5read.m" from SIRHEN; credit to Dan Dolan
% - Baseline filtering/removal with a bandstop filter
%
% Changed:
% - Renaming "RequiredFunctions" to "private"
% -- This should make all enclosed scripts automatically in path
% -- Thanks to Dan Dolan for advice!
%
% ============================ (March 1, 2023) ============================
% Added:
% - Option to correct up-shifted signals
% -- Additional yes/no option in first pop-up window
% -- Early step to fit baseline with very long durations
% - All prompts are now case insensitive
%
% ============================ (August 17, 2022) ==========================
% Added:
% - Method to read .dig files
% -- Borrowed from SIRHEN; credit to Dan Dolan
%
% ============================ (December 9, 2021) =========================
% Added:
% - Better method to read in .dat and .csv files
% -- any columns containing NaN data are removed
%
% ============================ (October 25, 2021) =========================
% Added:
% - Way to read *.isf files
%
% ============================ (December 12, 2020) ========================
% HiFiPDV first released as an embedded tarball on Page 181 of PhD thesis.
% Copyright Travis J Voorhees 2020
% See license file in /private/associated_licenses/ for details
%
% ============================ (August 7, 2020) ===========================
% Changed:
% - Example spectrogram is now made with longest duration length (100ns)
% - Example spectrogram bin-to-bin resolution is 10 ns
% - Increased linewidth of output plots to 2
% Added:
% - Window correction factor options: LiF, PMMA, more coming 'someday'
% ============================ (June 16, 2020) ============================
% Added:
% - Option to output Rise Times (10-50-90) of autodetected front
% - You can now select whether or not to run Gaussian fitting process
% - Improved method for extrapolation of user inputted ROI bounds
% -- No longer extrapolates beyond greatest or lowest user input
% - User ROI bounds are now shown in final spectrogram plot
%
% Removed:
% - Checking for "MonteCarlo" on the path (more user naming freedom)
% -- Code still looks for "RequiredFunctions" directory
%
% Considering working on:
% - Moving user prompts from a pop-up window to being inside the actual
% figure UI (not sure how this would affect performance).
% - Lowering RAM usage (laptop friendly)
% =========================================================================
% *************************************************************************
% ****************************** BEGIN CODE *******************************
% *************************************************************************
%% Start fresh
close all; clear; clc
%% ADVANCED USER OPTIONS:
% These are advanced options that users can vary if they wish.
% Default values are:
% >> SPECTROGRAM_COLORMAP = jet(256); % 'jet' colormap with 8-bit color
% >> HEATMAP_COLORMAP = hot(256); % 'hot' colormap with 8-bit color
% >> FILT_AMT = 0.01; % ±1 percent baseline filter
SPECTROGRAM_COLORMAP = jet(256); % colormap for spectrograms
HEATMAP_COLORMAP = hot(256); % colormap for heatmap
FILT_AMT = 0.01; % ± percentage of baseline filter, 0.01 = ±1 % filter
%% Initiate a parallel CPU pool if one hasn't already been setup
Pool = gcp; % Get Current Parallel pool
%% Read in the signal
% Use a GUI to select the PDV data file
% Filetypes supported: *.trc, *.csv, *.txt
[File,Path] = uigetfile({'*.trc;*.csv;*.txt;*.wfm;*.isf;*.dig;*.h5;*.H5'},'MultiSelect','off');
if File==0 % User did not select anything or hit cancel
clear;
error('User selected cancel or closed window. Process terminating now.');
else
[~,~,FileFormat] = fileparts([Path File]); % updated from just '*.xxx'
if ~FileFormat(1) == '.'
error('File format could not be determined.')
else
end
disp(['User selected file: ' File]);
end
% Set up a GUI for user to input quick cropping parameters & Gauss select
prompt = {'Time Units (s,us,ns):','Approximate Start Time (μs):',...
'Approximate End Time (μs):','Laser Wavelength (nm):',...
'Run Gaussian Fit? (yes/no)', 'Calculate 10-50-90 Rise Times (yes/no)?',...
'Window Material:', 'Correct for frequency-shifting? (up/down/no)',...
'Filter out signal baseline? (yes/no)'};
dlgtitle = 'Data Import Settings';
dims = [1 45];
definput = {'s','0','300','1550','yes','no','none','no','no','no'};
% Launch a GUI for the user to input cropping parameters
C = inputdlg(prompt,dlgtitle,dims,definput);
% Determine Conversion Factor to Convert Time-Scale to Nanoseconds
switch lower(C{1})
case {'s','seconds','sec','1'}
ConvFactor = 1e9;
case {'ns','nanoseconds','nano','1e-9'}
ConvFactor = 1;
case {'\mus','mus','microseconds','micro','1e-6'}
ConvFactor = 1e3;
case {'ms','mili','milli','miliseconds','milliseconds','1e-3'}
FirstTimePoint = Data(1,1)*1e3;
Time = Data(:,1).*1e3;
disp('Are you sure your time format is milliseconds?');
otherwise
error('foo:bar',['========== USER DID NOT ENTER USEABLE TIMESCALE ==========' ...
'\nACCEPTABLE TIMESCALES: s, ns, us, \mus, seconds, nanoseconds, microseconds, nano, micro' ...
'\nNOTE: TIME-SCALE ENTRIES ARE CASE-SENSITIVE, USE LOWER-CASE']);
end
% Determine if user wants to run Gaussian fitting or not
switch lower(C{5})
case {'yes','y','1','yes I want to','yes please','true'}
RunGauss = 1;
disp('Gaussian fitting process selected.');
case {'no','n','0','no thank you','no I do not want to','false'}
RunGauss = 0;
disp('Gaussian fitting process not selected.');
case {''}
RunGauss = 0;
disp('User left Gaussian fitting question blank.');
disp('Gaussian fitting process will not occur.');
otherwise
error('foo:bar',['========== GAUSSIAN FITTING ANSWER NOT ACCEPTED ==========' ...
'\nACCEPTABLE ANSWERS: y, n, yes, no, Yes, No, YES, NO, etc.' ...
'\nUser entered the following answer: ' C{5}]);
end
% Determine the window correction factor (if needed)
switch lower(C{7})
case {'yes','y','not sure'}
WFdlg = inputdlg('What is the window correction factor?','Enter Window Correction Factor',[1 45],{'0.7895'});
if WFdlg == ''
error('User did not enter a window correction factor when asked');
else
end
WF.a = str2double(WFdlg);
WF.b = 1;
case {'no','NO','none','No','N/A'}
disp('No window correction factor applied')
WF.a = 1;
WF.b = 1;
case {'LiF'}
WF.a = 0.7895;
WF.b = 0.9918;
case {'PMMA'}
disp(' ** PMMA requires a negligible window correction factor.')
WF.a = 1; % negligible
WF.b = 1; % negligible
otherwise
WFdlg = inputdlg('What is the window correction factor?','Enter Window Correction Factor',[1 45],'0.7895');
if WFdlg == ''
error('User did not enter a window correction factor when asked');
else
end
WF.a = str2double(WF);
WF.b = 1;
end
% Determine if the user wants to correct for frequency shifting or not:
switch lower(C{8})
% The user might input 'yes', which is not descriptive enough
% Ask them to specify if the signal is "up" or "down" shifted:
case {'yes','y','ya'}
FSdlg = inputdlg('Is the signal up-shifted or down-shifted?','Is the signal up-shifted or down-shifted?',[1 45],{'up'});
if strcmpi(FSdlg,'')
warning('User did not state if signal was up- or down-shifted');
warning('Proceeding with analysis, assuming signal is not frequency shifted');
warning('RESTART ANALYSIS IF THIS IS INCORRECT');
C{8} = 'no';
else
C{8} = FSdlg{1};
end
if strcmpi(FSdlg,'yes')
error('Please retry. "Yes" is not a valid answer for direction of frequency shifting.')
end
end
switch lower(C{8})
case {'no','none','n/a'}
disp('Signal is *NOT* frequency-shifted.')
case {'up','up-shifted','UP','Up','upshifted'}
SHIFT_SIGN = 1; % will be used in later section
case {'down','downshifted','down-shifted','DOWN','Down'}
SHIFT_SIGN = -1; % will be used in later section
warning('Code is not well suited for down-shifted signals yet');
otherwise
warning(['Frequency-shifting answer could not be interpreted. User input: ' C{8}])
error('Please retry.')
end
% Determine if the user wants to filter out the baseline
switch lower(C{9})
case {'yes','y','ya','true'}
disp('User selected to filter baseline.')
disp('Bandstop filter will be applied from 99-101% of baseline')
BaseRemoval = 1;
case {'no','n','false','none','n/a'}
disp('User selected not to filter baseline from signal')
BaseRemoval = 0;
otherwise
warning('Bad input for baseline filtering question')
error(['When asked if signal should be filtered, user input: ' C{9}])
end
% Read in the user inputs
Start = str2double(C{2})*1e3; % convert to nanoseconds
End = str2double(C{3})*1e3; % convert to nanoseconds
Wavelength = str2double(C{4});
% Read in the data
switch FileFormat
case '.trc'
data = ReadLeCroyBinaryWaveform([Path File]);
Signal = [data.x data.y];
case '.txt'
Signal = textscan([Path File],'%f%f');
Signal(any(isnan(Signal), 2), :) = []; % remove any rows with NaN
case {'.csv','.dat'}
Signal = readmatrix([Path File]);
Signal(:,isnan(mean(Signal))) = [];
warning("CSV data read in such that column 1 is time, column 2 is voltage:");
disp(Signal(1:10,:));
disp("...")
warning("If this is incorrect, please reformat input data")
if width(Signal) == 1
error("Only one column of data was read in. Time and Voltage is needed.")
else
end
case '.wfm'
[data.y,data.x] = wfm2read([Path File]);
Signal = [data.x data.y];
case '.isf'
data = isfread([Path File]);
Signal = [data.x data.y];
case '.dig'
[data.y,data.x] = digread([Path File]);
Signal = [data.x data.y];
case {'.h5','.H5'}
% HDF5 can hold a lot of channels, have the user select one
H5dlg = inputdlg('HDF5 Detected. Which Channel are you interested in?',...
'HDF5 Detected, Choose Channel.',[1 55],{'1'});
% If user does not input a selection, default to channel 1
if strcmpi(H5dlg,'')
warning('User did not select a channel. Default is 1.');
warning('Proceeding analysis with Channel 1 selected.');
warning('** RESTART ANALYSIS IF THIS IS INCORRECT **');
Channel = 1;
else
% Try to evaluate user's channel selection
% Throw an error if non-numeric value is entered
try
Channel = eval(H5dlg{1});
catch
error(['User input "' H5dlg{1} '" as channel selection']);
end
end
try
% Try to read the data
[data.y,data.x] = agilent_h5read([Path File],Channel);
Signal = [data.x data.y];
warning('Reading HDF5 format is a *NEW* feature in HiFiPDV.');
warning('Contact [email protected] if bugs arise.');
catch
% Throw an error if data cannot be read.
error(['Cannot read Channel ' num2str(Channel) ' in ' Path File])
end
otherwise
error('foo:bar',['===== USER DID NOT CHOOSE COMPATIBLE FILE FORMAT =====' ...
'\nACCEPTABLE FORMATS: .trc, .csv, .txt, .wfm, .h5']);
end
% Convert data timescale to nanoseconds
Signal(:,1) = Signal(:,1).*ConvFactor;
% Determine Signal Data Sampling Rate (must be done prior to cropping)
SampleRate = 1e9/(mean(diff(Signal(:,1)))); %1e9 b/c time is in nanoseconds
% if Sample Rate is negative, then flip data across the X-axis
if sign(SampleRate)==-1
Signal = flipud(Signal);
SampleRate = 1e9/(mean(diff(Signal(:,1)))); %recalculate
else
end
% Spit out an error if sample rate is slower than 2 points/ns
if SampleRate<2.49e9
error('foo:bar',['Signal Sample Rate is ' num2str(SampleRate./1e9,'%.2e') ...
' points/nanosecond\nSample Rate must be 2.00 points/nanosecond or' ...
' greater for this code to run correctly.'])
else
end
% Crop the data
[~,FirstInd] = min(abs(Signal(:,1)-Start)); % Finds closest value to start limit
[~,LastInd]=min(abs(Signal(:,1)-End)); % Finds closest value to end limit
Signal = Signal(FirstInd:LastInd,:);
% Time-shift signal
FirstTimePoint = Signal(1,1);
Signal(:,1) = Signal(:,1)-FirstTimePoint; % nanoseconds, starting at zero
%% Create STFT inputs and pre-allocate arrays
% REMEMBER: time axis is now in nanoseconds, starting at zero
% First time point is "FirstTimePoint" (also in nanoseconds)
% Detector voltage is second column of signal, time is first column
DetectorVoltage = Signal(:,2);
% Determine the number of data points in a 1 nanosecond duration
Num1ns = SampleRate/1e9;
NumDataPointsFor1ns = round(Num1ns); % needs to be integer value
% NumDataPointsFor1ns needs to be even
% if this value is odd, then add 1 and let the user know the new time-resolution
if mod(NumDataPointsFor1ns,2) == 1 % odd case
disp(['Number of data points per 1 ns is: ' num2str(Num1ns,'%.1f')]);
NumDataPointsFor1ns = NumDataPointsFor1ns+1;
disp(['Point-to-point time resolution set to: ' num2str(2.5*(NumDataPointsFor1ns...
/Num1ns),'%.4f') 'ns']);
else % Even case
disp(['Number of data points per 1 ns is: ' num2str(Num1ns,'%.1f')]);
disp(['Point-to-point time resolution set to: ' num2str(2.5*(NumDataPointsFor1ns...
/Num1ns),'%.4f') 'ns']);
end
% Create a range of duration lengths: ~5 ns to ~100 ns
RangeOfDurationLengths = [5:5:100].*NumDataPointsFor1ns;
% Define overlap as duration - 2 ns, this creates a 2 ns time-step
RangeOfOverlaps = RangeOfDurationLengths-round(2.5*NumDataPointsFor1ns);
% Create a variable for number of frequency bins
nFreqs = 2^15; % Reduce this value if RAM is being maxed out
%% Determine the reference baseline (for frequency shifted signals)
% If the signal is frequency-shifted, have the user draw lines around the
% baseline region (the region that exists before the region of interest)
% Only perform this process if the SHIFT_SIGN variable exists
if exist('SHIFT_SIGN')
BaseClr = 'm';
% Make temporary inputs for spectrogram
TmpWin = hamming(RangeOfDurationLengths(end)); % longest window
TmpOvr = RangeOfOverlaps(end); % corresponding overlap
% Calculate a temporary spectrogram
[zexample,f,texample] = spectrogram(DetectorVoltage,TmpWin,TmpOvr,nFreqs,SampleRate);
zexample = real(zexample.*conj(zexample)); % get rid of imaginary numbers
zexample = zexample./max(zexample(:)); % normalize globally
zexample = 10.*log10(zexample); % convert to dB for plot
texample = texample.*1e6 + FirstTimePoint/1e3; % convert time to microseconds and shift
f = f./1e9; % convert to GHz
% Show the spectrogram, get user input for cropping prior to lineout calcs
figure
hold on
imagesc(texample,f.*(Wavelength/2),zexample)
%colormap(jet(16)) % option for less color depth if X11 forwarding
colormap(SPECTROGRAM_COLORMAP)
colorbar
caxis([-60 0])
axis([-inf +inf -25 +inf]) % we don't usually see velocities over 5,000 m/s
xlabel('Detector Time (\mus)')
ylabel('Signal Velocity (m/s)')
title('Reference/Baseline Selection')
hold off
% Have the user zoom into their region of interest
zoom on
Prompt{1} = 'ZOOM to your BASELINE/REFERENCE REGION';
Prompt{2} = '';
Prompt{3} = ' - LEFT-CLICK and DRAG a box around this region.';
Prompt{4} = '';
Prompt{5} = ' - PRESS ENTER when you are finished.';
Prompt{6} = '';
Prompt{7} = ' - DOUBLE LEFT-CLICK to undo/reset zoom level.';
mb = msgbox(Prompt,'Zoom to REFERENCE region');
pause
zoom off
delete(mb)
clear mb
% Upper Bound of Baseline
PromptUPb{1} = 'Draw an UPPER BOUND for your BASELINE.';
PromptUPb{2} = '';
PromptUPb{3} = ' - Full screen viewing is recommended for higher accuracy.';
PromptUPb{4} = '';
PromptUPb{5} = ' - LEFT-CLICK to place points, RIGHT-CLICK when you are finished.';
mb = msgbox(PromptUPb,'Draw an UPPER BOUND for your BASELINE');
roiUPb = images.roi.Polyline('color',BaseClr);
draw(roiUPb);
ROIUPb = roiUPb.Position;
ROIUPb = sortrows(ROIUPb); % must be in order for future interpolation
% Clean up the ROI and replot it
% Don't allow multiple y-values per x-value
ROIxstepUPb = ROIUPb(2:end,1)-ROIUPb(1:(end-1),1);
if mean(ROIxstepUPb)>0 % moving rightwards
ROIxstepUPb = [1;ROIxstepUPb];
ROIUPb(ROIxstepUPb<=0,:)=[];
else
ROIxstepUPb = [-1;ROIxstepUPb];
ROIUPb(ROIxstepUPb>=0,:)=[]; % moving leftwards
end
roiUPb.Position=ROIUPb; %update the drawing on the figure
ROIUPxb = ROIUPb(:,1).*1e3 -FirstTimePoint; % convert to nanoseconds
ROIUPyb = ROIUPb(:,2); % frequency in GHz
delete(mb)
clear mb
% Lower Bound of Baseline
PromptLOWb = PromptUPb;
PromptLOWb{1} = 'Draw a LOWER BOUND for your BASELINE.';
mb = msgbox(PromptLOWb,'Draw LOWER BOUND for BASELINE');
roiLOWb = images.roi.Polyline('color',BaseClr);
draw(roiLOWb);
ROILOWb = roiLOWb.Position;
ROILOWb = sortrows(ROILOWb); % must be in order for future interpolation
% Clean up the ROI and replot it
% Don't allow multiple y-values per x-value
ROIxstepLOWb = ROILOWb(2:end,1)-ROILOWb(1:(end-1),1);
if mean(ROIxstepLOWb)>0 % moving rightwards
ROIxstepLOWb = [1;ROIxstepLOWb];
ROILOWb(ROIxstepLOWb<=0,:)=[];
else
ROIxstepLOWb = [-1;ROIxstepLOWb];
ROILOWb(ROIxstepLOWb>=0,:)=[]; % moving leftwards
end
roiLOWb.Position=ROILOWb; %update the drawing on the figure
ROILOWxb = ROILOWb(:,1).*1e3 - FirstTimePoint; % convert to nanoseconds, shift
ROILOWyb = ROILOWb(:,2); % frequency in GHz
delete(mb)
clear mb
% Check if user is okay with their frequency bounds before proceeding
UserPrompt2 = questdlg('Are you satisfied with the current BASELINE frequency bounds?',...
'Use These Bounds?','Yes','No','Cancel','No');
switch UserPrompt2
case 'Yes'
disp('Baseline freqency bounds entered, beginning baseline extraction process.');
disp(' ');
case 'No'
while strcmpi(UserPrompt2,'No')
delete(roiUPb);
roiUPb = images.roi.Polyline('color',BaseClr);
delete(roiLOWb);
roiLOWb = images.roi.Polyline('color',BaseClr);
% Redo the upper bound
mb = msgbox(PromptUPb,'Draw UPPER BOUND for your BASELINE');
draw(roiUPb);
ROIUPb = roiUPb.Position;
ROIUPb = sortrows(ROIUPb); % must be in order for future interpolation
% Clean up the ROI and replot it
% Don't allow multiple y-values per x-value
ROIxstepUPb = ROIUPb(2:end,1)-ROIUPb(1:(end-1),1);
if mean(ROIxstepUPb)>0 % moving rightwards
ROIxstepUPb = [1;ROIxstepUPb];
ROIUPb(ROIxstepUPb<=0,:)=[];
else
ROIxstepUPb = [-1;ROIxstepUPb];
ROIUPb(ROIxstepUPb>=0,:)=[]; % moving leftwards
end
roiUPb.Position=ROIUPb; %update the drawing on the figure
ROIUPxb = ROIUPb(:,1).*1e3 -FirstTimePoint; % convert to nanoseconds
ROIUPyb = ROIUPb(:,2); % frequency in GHz
delete(mb)
clear mb
% Redo the lower bound
mb = msgbox(PromptLOWb,'Draw LOWER BOUND for your BASELINE');
roiLOWb = images.roi.Polyline('color',BaseClr);
draw(roiLOWb);
ROILOWb = roiLOWb.Position;
ROILOWb = sortrows(ROILOWb); % must be in order for future interpolation
% Clean up the ROI and replot it
% Don't allow multiple y-values per x-value
ROIxstepLOWb = ROILOWb(2:end,1)-ROILOWb(1:(end-1),1);
if mean(ROIxstepLOWb)>0 % moving rightwards
ROIxstepLOWb = [1;ROIxstepLOWb];
ROILOWb(ROIxstepLOWb<=0,:)=[];
else
ROIxstepLOWb = [-1;ROIxstepLOWb];
ROILOWb(ROIxstepLOWb>=0,:)=[]; % moving leftwards
end
roiLOWb.Position=ROILOWb; %update the drawing on the figure
ROILOWxb = ROILOWb(:,1).*1e3 - FirstTimePoint; % convert to nanoseconds, shift
ROILOWyb = ROILOWb(:,2); % frequency in GHz
delete(mb)
clear mb
% Ask if the new bounds are satisfactory
UserPrompt2 = questdlg('Are you satisfied with the new BASELINE frequency bounds?',...
'Use These Bounds?','Yes','No','Cancel','No');
end
case 'Cancel'
close all; clear; clc;
error('User selected cancel. Terminating processes now.')
end
% Convert y-axis values to GHz
ROIUPyb = ROIUPyb./(Wavelength/2);
ROILOWyb = ROILOWyb./(Wavelength/2);
close all % close the figure
% Crop out the baseline for further analysis
% Determine minimum and maximum time in drawn ROI limits
BasexMin = min([ROIUPxb;ROILOWxb]);
BasexMax = max([ROIUPxb;ROILOWxb]);
BaseyMin = min([ROIUPyb;ROILOWyb]);
BaseyMax = max([ROIUPyb;ROILOWyb]);
% Crop signal to ROI left and right limits
[~,BaseMinInd] = min(abs(Signal(:,1)-BasexMin));
[~,BaseMaxInd] = min(abs(Signal(:,1)-BasexMax));
BASELINE=Signal(BaseMinInd:BaseMaxInd,:);
% Redefine DetectorVoltage with new time limits
BaseTime = BASELINE(:,1);
BaseVoltage = BASELINE(:,2);
% Calculate a SINGLE velocity within the baseline region of interest
nFreqsB = nFreqs*10; % 10 fold increase for the baseline determination
BASEfft.x = (1:nFreqs).*SampleRate./nFreqsB./1e9; % Convert to GHz
BASEfft.ham = fft(BaseVoltage.*hamming(length(BaseVoltage)),nFreqsB); % FFT with hamming window
BASEfft.blk = fft(BaseVoltage.*blackman(length(BaseVoltage)),nFreqsB); % FFT with blackman window
BASEfft.han = fft(BaseVoltage.*hann(length(BaseVoltage)),nFreqsB); % FFT with hann window
% Convert to real space
BASEfft.ham = real(BASEfft.ham.*conj(BASEfft.ham)); % abs needs to be applied FIRST
BASEfft.blk = real(BASEfft.blk.*conj(BASEfft.blk));
BASEfft.han = real(BASEfft.han.*conj(BASEfft.han));
% Crop each of the FFTs to just the region of interest (user input)
[~,BASEfft.Low] = min(abs(BASEfft.x-BaseyMin)); % find the low index
[~,BASEfft.High] = min(abs(BASEfft.x-BaseyMax)); % find the high index
BASEfft.ham = BASEfft.ham(BASEfft.Low:BASEfft.High); % crop
BASEfft.blk = BASEfft.blk(BASEfft.Low:BASEfft.High); % crop
BASEfft.han = BASEfft.han(BASEfft.Low:BASEfft.High); % crop
BASEfft.x = BASEfft.x(BASEfft.Low:BASEfft.High); % crop
% Normalize after cropping
BASEfft.ham = BASEfft.ham/max(max(BASEfft.ham));
BASEfft.blk = BASEfft.blk/max(max(BASEfft.blk));
BASEfft.han = BASEfft.han/max(max(BASEfft.han));
% Extract the peak location from each FFT
BASEfft.MAX = zeros(1,3); % max extraction
BASEfft.CEN = zeros(1,3); % centroid extraction
BASEfft.GAS = zeros(1,3); % Gaussian extraction
BASEfft.ALL = zeros(1,9); % 3 windows, 3 extractions, 9 answers
% MAXIMUM
[~,tmp] = max(BASEfft.ham);
BASEfft.MAX(1) = BASEfft.x(tmp); clear tmp
[~,tmp] = max(BASEfft.blk);
BASEfft.MAX(2) = BASEfft.x(tmp); clear tmp
[~,tmp] = max(BASEfft.han);
BASEfft.MAX(3) = BASEfft.x(tmp); clear tmp
% CENTROID
% f_centerofmass is equal to {int(f*Z)df / int(Z)df}
% In other words: integral of curve*x / area under curve
BASEfft.CEN(1) = trapz(BASEfft.x,BASEfft.x.*BASEfft.ham')./...
trapz(BASEfft.x,BASEfft.ham);
BASEfft.CEN(2) = trapz(BASEfft.x,BASEfft.x.*BASEfft.blk')./...
trapz(BASEfft.x,BASEfft.blk);
BASEfft.CEN(3) = trapz(BASEfft.x,BASEfft.x.*BASEfft.han')./...
trapz(BASEfft.x,BASEfft.han);
% Gaussian fit
% NOTE: Uses a different Gaussian fit than Spectrogram and Heatmap extraction
% Why? Old code that hasn't proven to be an issue (yet). May update in next release.
BaseFit.Xf = BASEfft.x(1):(0.01*2/Wavelength):BASEfft.x(end);
BaseFit.Xv = BaseFit.Xf.*(Wavelength/2);
tmp = fit(BASEfft.x',BASEfft.ham,'gauss1');
BaseFit.ham = tmp(BaseFit.Xf); % do this for plots later
BASEfft.GAS(1) = tmp.b1; clear tmp
tmp = fit(BASEfft.x',BASEfft.blk,'gauss1');
BaseFit.blk = tmp(BaseFit.Xf); % do this for plots later
BASEfft.GAS(2) = tmp.b1; clear tmp
tmp = fit(BASEfft.x',BASEfft.han,'gauss1');
BaseFit.han = tmp(BaseFit.Xf); % do this for plots later
BASEfft.GAS(3) = tmp.b1; clear tmp
BASEfft.ALL = [BASEfft.MAX BASEfft.CEN BASEfft.GAS];
SHIFT_AMT = SHIFT_SIGN*mean(BASEfft.ALL)*Wavelength/2; % converted to velocity
SHIFT_UNC = 1.95*std(BASEfft.ALL)*Wavelength/2;
Basefilter = [1-FILT_AMT 1+FILT_AMT].*(SHIFT_AMT*2/Wavelength).*1e9; % Bandstop filter limits
else
disp(''); % signal is not shifted
Basefilter = [0 0.5e8]; % bandstop filter limits (0 to 50 MHz/39 m/s)
% NOTE TO USERS:
% - feel free to change this value, it is arbitrary
% - Remember: Velocity = Signal_Frequency * (Laser_Wavelength/2)
end
%% Apply a bandstop filter to the filter out the baseline
if BaseRemoval == 1
DetectorVoltage = bandstop(DetectorVoltage,Basefilter,SampleRate); % simple bandstop filter
else
end
%% User input Region of Interest for next step in calculations
% Also have them select the baseline if data is frequency shifted
% Define line colors for these steps
UpClr = 'w'; % line color for signal region of interest
LowClr = 'w'; % line color for signal region of interest
% Make temporary inputs for spectrogram
IndexForExample = 15; % default is end, but you can select a different spectrogram for the example
% NOTE: Consider allowing users to choose their own spectrogram for ROI
TmpWin = hamming(RangeOfDurationLengths(IndexForExample)); % Use Hamming for example
TmpOvr = RangeOfOverlaps(IndexForExample); % Corresponding overlap
% Calculate a temporary spectrogram
[zexample,f,texample] = spectrogram(DetectorVoltage,TmpWin,TmpOvr,nFreqs,SampleRate);
zexample = real(zexample.*conj(zexample)); % get rid of imaginary numbers
zexample = zexample./max(zexample(:)); % normalize globally
zexample = 10.*log10(zexample); % convert to dB for plot
texample = texample.*1e6 + FirstTimePoint/1e3; % convert time to microseconds and shift
f = f./1e9; % convert to GHz
% Show the spectrogram, get user input for cropping prior to lineout calcs
figure
hold on
imagesc(texample,f.*(Wavelength/2),zexample)
colormap(SPECTROGRAM_COLORMAP)
%colormap(jet(16)) % option for less color depth if X11 forwarding
colorbar
caxis([-60 0])
axis([-inf +inf -25 +inf]) % we don't usually see velocities over 5,000 m/s
xlabel('Detector Time (\mus)')
ylabel('Signal Velocity (m/s)')
hold off
% Have the user zoom into their region of interest
zoom on
Prompt{1} = 'ZOOM to your SIGNAL region of interest';
Prompt{2} = '';
Prompt{3} = ' - LEFT-CLICK and DRAG a box around your SIGNAL region of interest.';
Prompt{4} = '';
Prompt{5} = ' - PRESS ENTER when you are finished.';
Prompt{6} = '';
Prompt{7} = ' - DOUBLE LEFT-CLICK to undo/reset zoom level.';
Prompt{8} = '';
Prompt{9} = ' Note: If your signal is frequency-shifted, include baseline';
mb = msgbox(Prompt,'Zoom to SIGNAL region of interest');
pause
zoom off
delete(mb)
clear mb
% =====
% User input for frequency limits/region of interest
% =====
% Upper Bound of Actual Signal
PromptUP{1} = 'Draw an UPPER BOUND for your SIGNAL.';
PromptUP{2} = '';
PromptUP{3} = ' - Full screen viewing is recommended for higher accuracy.';
PromptUP{4} = '';
PromptUP{5} = ' - LEFT-CLICK to place points, RIGHT-CLICK when you are finished.';
mb = msgbox(PromptUP,'Draw UPPER BOUND for Calculation');
roiUP = images.roi.Polyline('color',UpClr);
draw(roiUP);
ROIUP = roiUP.Position;
ROIUP = sortrows(ROIUP); % must be in order for future interpolation
% Clean up the ROI and replot it
% Don't allow multiple y-values per x-value
ROIxstepUP = ROIUP(2:end,1)-ROIUP(1:(end-1),1);
if mean(ROIxstepUP)>0 % moving rightwards
ROIxstepUP = [1;ROIxstepUP];
ROIUP(ROIxstepUP<=0,:)=[];
else
ROIxstepUP = [-1;ROIxstepUP];
ROIUP(ROIxstepUP>=0,:)=[]; % moving leftwards
end
roiUP.Position=ROIUP; %update the drawing on the figure
ROIUPx = ROIUP(:,1).*1e3 -FirstTimePoint; % convert to nanoseconds
ROIUPy = ROIUP(:,2); % frequency in GHz
delete(mb)
clear mb
% Lower Bound of Actual Signal
PromptLOW = PromptUP;
PromptLOW{1} = 'Draw a LOWER BOUND for your SIGNAL.';
mb = msgbox(PromptLOW,'Draw LOWER BOUND for Calculation');
roiLOW = images.roi.Polyline('color',LowClr);
draw(roiLOW);
ROILOW = roiLOW.Position;
ROILOW = sortrows(ROILOW); % must be in order for future interpolation
% Clean up the ROI and replot it
% Don't allow multiple y-values per x-value
ROIxstepLOW = ROILOW(2:end,1)-ROILOW(1:(end-1),1);
if mean(ROIxstepLOW)>0 % moving rightwards
ROIxstepLOW = [1;ROIxstepLOW];
ROILOW(ROIxstepLOW<=0,:)=[];
else
ROIxstepLOW = [-1;ROIxstepLOW];
ROILOW(ROIxstepLOW>=0,:)=[]; % moving leftwards
end
roiLOW.Position=ROILOW; %update the drawing on the figure
ROILOWx = ROILOW(:,1).*1e3 - FirstTimePoint; % convert to nanoseconds, shift
ROILOWy = ROILOW(:,2); % frequency in GHz
delete(mb)
clear mb
% Check if user is okay with their frequency bounds before proceeding
UserPrompt2 = questdlg('Are you satisfied with the current signal frequency bounds?',...
'Use These Bounds?','Yes','No','Cancel','No');
switch UserPrompt2
case 'Yes'
disp('Freqency bounds entered, beginning extraction process.');
case 'No'
while strcmpi(UserPrompt2,'No')
delete(roiUP);
roiUP = images.roi.Polyline('color',UpClr);
delete(roiLOW);
roiLOW = images.roi.Polyline('color',LowClr);
mb = msgbox(PromptUP,'Draw UPPER BOUND for your SIGNAL');
% re-do the upper bound
draw(roiUP);
ROIUP = roiUP.Position;
ROIUP = sortrows(ROIUP); % must be in order for future interpolation
% Clean up the ROI and replot it
% Don't allow multiple y-values per x-value
ROIxstepUP = ROIUP(2:end,1)-ROIUP(1:(end-1),1);
if mean(ROIxstepUP)>0 % moving rightwards
ROIxstepUP = [1;ROIxstepUP];
ROIUP(ROIxstepUP<=0,:)=[];
else
ROIxstepUP = [-1;ROIxstepUP];
ROIUP(ROIxstepUP>=0,:)=[]; % moving leftwards
end
roiUP.Position=ROIUP; %update the drawing on the figure
ROIUPx = ROIUP(:,1).*1e3 - FirstTimePoint; % convert to nanoseconds, shift
ROIUPy = ROIUP(:,2); % frequency in GHz
delete(mb)
clear mb
% re-do the lower bound
mb = msgbox(PromptLOW,'Draw LOWER BOUND for your SIGNAL');
draw(roiLOW);
ROILOW = roiLOW.Position;
ROILOW = sortrows(ROILOW); % must be in order for future interpolation
% Clean up the ROI and replot it
% Don't allow multiple y-values per x-value
ROIxstepLOW = ROILOW(2:end,1)-ROILOW(1:(end-1),1);
if mean(ROIxstepLOW)>0 % moving rightwards
ROIxstepLOW = [1;ROIxstepLOW];
ROILOW(ROIxstepLOW<=0,:)=[];
else
ROIxstepLOW = [-1;ROIxstepLOW];
ROILOW(ROIxstepLOW>=0,:)=[]; % moving leftwards
end
roiLOW.Position=ROILOW; %update the drawing on the figure
ROILOWx = ROILOW(:,1).*1e3 - FirstTimePoint; % convert to nanoseconds, shift
ROILOWy = ROILOW(:,2); % frequency in GHz
delete(mb)
clear mb
% Ask if the new bounds are satisfactory
UserPrompt2 = questdlg('Are you satisfied with the new frequency bounds?',...
'Use These Bounds?','Yes','No','Cancel','No');
end
case 'Cancel'
close all; clear; clc;
error('User selected cancel. Terminating processes now.')
end
% Convert ROIUPy and ROILOWy to GHz
ROIUPy = ROIUPy./(Wavelength/2);
ROILOWy = ROILOWy./(Wavelength/2);
close all % close the figure
% clear the temporary variables to free up memory
% Don't clear f, it will be used in all future calcs
% Don't clear texample or zexample, they will be used in an output plot
clear TmpWin TmpOvr roi Prompt Prompt2 ROI
%% Crop the signal time according to the ROI drawn by the user
% Doing this will reduce memory costs in spectrogram calculations
% Determine minimum and maximum time in drawn ROI limits
ROIxMin = min([ROIUPx;ROILOWx]);
ROIxMax = max([ROIUPx;ROILOWx]);
% Crop signal to ROI left and right limits
[~,LeftInd] = min(abs(Signal(:,1)-ROIxMin));
[~,RightInd] = min(abs(Signal(:,1)-ROIxMax));
Signal=Signal(LeftInd:RightInd,:);
% Redefine DetectorVoltage with new time limits
DetectorVoltage = Signal(:,2);
% Determine indices of frequency that are closest to bounds
% Preallocate variables populated in for loop:
fIndUP = zeros(1,length(ROIUPx));
fIndLOW = zeros(1,length(ROILOWx));
% Upper ROI bound
parfor i = 1:length(ROIUPx)
[~,fIndUP(i)] = min(abs(f-ROIUPy(i)));
end
% Lower ROI bound
parfor i = 1:length(ROILOWx)
[~,fIndLOW(i)] = min(abs(f-ROILOWy(i)));
end
% NOTE: frequency indices can be used directly in spectrograms
BottomInd = min([fIndUP fIndLOW]);
UpperInd = max([fIndUP fIndLOW])+1; % add 1 to be conservative
% Vertically crop f and zexample, this will be useful later
zexample = zexample(1:UpperInd,:);
f = f(1:UpperInd,:); % Keep f in double precision
% Re-normalize zexample:
%zexample = zexample./max(max(zexample));
% Redefine ROI Upper and Lower bound frequencies using f
ROIUPy = f(fIndUP);
ROILOWy = f(fIndLOW); % remove (-) numbers, will be helpful in future steps
%% Prep for calculation loop: create progress bar and pre-allocate cells
% 60 varieties of spectrograms are to be calculated
% output variables will update after every calculation
% individual spectrograms are saved in cell format