diff --git a/change_log.txt b/change_log.txt
index 98cbe095..177f0098 100644
--- a/change_log.txt
+++ b/change_log.txt
@@ -2,6 +2,34 @@
 J. James Jun
 
 --------------------------------------------------------------------
+[2019/1/5: v3.2.5]
+Centered cluster average waveforms are properly saved. 
+  `trWav_spk_clu` for filtered waveforms and `trWav_raw_clu` for raw waveforms
+  `tmrWav_spk_clu` and `tmrWav_raw_clu` have worked properly, which saves
+  waveforms at all sites.
+Selecting 'vpp' (peak-to-peak voltage) calcultes different values
+  depending on whether you are displaying filtered or raw waveforms.
+  When 'pca' or 'ppca' (private pca) features are selected, PCA values are computed for the filtered waveforms whether or not the raw waveforms are displayed.
+Manual GUI: The x-axis labels on the waveform view are displayed without overlapping for smaller sized monitors.
+  To avoid overlapping, the text size is decresased and text angle is rotated (new feature since Matlab 2016b).
+
+# New features
+Manual GUI: 'ppca' (private PCA) option is added under "Projection" menu.
+  When one unit is selected, filtered spike waveforms are projected to the first and second principal vectors of the selected unit for each site.
+  When two units are selected, filtered spike waveforms are projected to the first principal axes of the two selected units for each site.
+"jrc export-chan myparam.prm chan-list" exports channel(s) on "chan-list" to the Matlab workspace and to a binary file (_ch#.jrc).
+  To export more than one channel, use "," or ":" without space. 
+  For example, run "jrc export-chan 33,34,35" or "jrc export-chan 33:35" to export channels 33-35.
+"jrc export-lfp" to export the LFP waveforms (ordered by the site numbers) to the Matlab workspace
+  "mnLfp" has a nSamples x nSites dimension, and stores the LFP waveforms in the original data type (16-bit integer by default).
+  "mrSiteXY" has a nSites x 2 dimension, and stores x and y coordinate for the sites (in micrometers)
+  To extract a column of sites located at x=0, run "mnLfp(:, mrSiteXY(:,1)==0)"
+
+[2019/1/4: v3.2.4]
+# Bug fixes
+Manual GUI: projection display is correctly updated.
+  Previously the second unit selection (shown in red) required two right clicks to update.
+
 [2019/1/3: v3.2.3]
 Bugfix: incorrect "gather" function was called for R2017b when Bigdata toolbox is installed.
 Faster and more robust fft_clean operation (run if `fft_thresh` > 0)
diff --git a/default.prm b/default.prm
index 698c4741..6fe4e70e 100644
--- a/default.prm
+++ b/default.prm
@@ -16,7 +16,7 @@ fInverse_file = 0;          % Set to 1 to flip the polarity of the signal
 header_offset = 0;          % File header offset (set to 0 for files containing no header info (e.g. WHISPER format)
 
 % #Execution parameters
-version = 'v3.2.3';         % JRCLUST version. Updated on Jan 3, 2018
+version = 'v3.2.5';         % JRCLUST version. Updated on Jan 8, 2018
 fVerbose = 0;               % Verbose flag, set to 0 to suppress displaying extra information 
 fParfor = 1;                % Use Multiple CPU cores (if parallel processing toolbox is installed)
 fGpu = 1;                   % Use GPU if parallel processing toolbox is installed
@@ -50,7 +50,7 @@ fText = 1;                  % Show spike count per unit (*, press 'n' to toggle
 
 % #Pre-processing (run "jrc spikesort" after change)
 tlim_load = [];             % Time range to load in sec (uses all time range if set to empty)
-vcFilter = 'ndiff';         % {'ndiff', 'sgdiff', 'bandpass', 'fir1', 'user', 'none'}
+vcFilter = 'ndiff';         % {'ndiff', 'sgdiff', 'bandpass', 'fir1', 'user', 'fftdiff', 'none'}
 nDiff_filt = 2;			    % Differentiation filter for vcFilter='sgdiff', ignored otherwise. Set to [] to disable. 2n+1 samples used for centered differentiation
 freqLim = [300, 3000];      % Frequency cut-off limit for vcFilter='bandpass', ignored otherwise. 
 filtOrder = 3;              % Bandpass filter order
diff --git a/jrc3.m b/jrc3.m
index d5e2d738..b9c3193e 100644
--- a/jrc3.m
+++ b/jrc3.m
@@ -159,6 +159,8 @@
         trWav_raw = load_bin_(strrep(P.vcFile_prm, '.prm', '_spkraw.jrc'), 'int16', S0.dimm_spk);
         assignWorkspace_(trWav_raw);  
     case {'export-spkwav', 'spkwav'}, export_spkwav_(P, vcArg2); % export spike waveforms
+    case {'export-chan'}, export_chan_(P, vcArg2); % export channels
+    case {'export-car'}, export_car_(P, vcArg2); % export common average reference
     case {'export-spkwav-diff', 'spkwav-diff'}, export_spkwav_(P, vcArg2, 1); % export spike waveforms        
     case 'export-spkamp', export_spkamp_(P, vcArg2); %export microvolt unit
     case {'export-csv', 'exportcsv'}, export_csv_(P);
@@ -168,6 +170,7 @@
     case {'export-fet', 'export-features', 'export-feature'}, export_fet_(P);
     case 'export-diff', export_diff_(P); %spatial differentiation for two column probe
     case 'import-lfp', import_lfp_(P); 
+    case 'export-lfp', export_lfp_(P); 
     case 'drift', plot_drift_(P);  
     case 'plot-rd', plot_rd_(P);
     otherwise, fError = 1;
@@ -1510,7 +1513,7 @@ function detect_(P)
 P.vcFile_prm = vcFile_prm;
 assert_(isfield(P, 'vcFile'), sprintf('Check "%s" file syntax', vcFile_prm));
 
-if ~exist_file_(P.vcFile)
+if ~exist_file_(P.vcFile) && isempty(get_(P, 'csFile_merge'))
     P.vcFile = replacePath_(P.vcFile, vcFile_prm);
     if ~exist_file_(P.vcFile)
         fprintf('vcFile not specified. Assuming multi-file format ''csFiles_merge''.\n');
@@ -3297,8 +3300,9 @@ function close_hFig_traces_(hFig, event)
         for i=1:size(mnWav2,2)
             mnWav2(:,i) = conv(mnWav2(:,i), vrFilter, 'same'); 
         end        
-    case 'ndiff'
-        mnWav2 = ndiff_(mnWav2, P.nDiff_filt);
+    case 'ndiff', mnWav2 = ndiff_(mnWav2, P.nDiff_filt);
+    case 'fftdiff', mnWav2 = fftdiff_(mnWav2, P);        
+%     case 'fftdiff', mnWav2 = fftdiff__(gather_(mnWav2), P.freqLim(2)/P.sRateHz/2);        
     case {'sgdiff', 'sgfilt'}
         mnWav2 = sgfilt_(mnWav2, P.nDiff_filt);
     case 'bandpass'
@@ -3327,6 +3331,60 @@ function close_hFig_traces_(hFig, event)
 end %func
 
 
+%--------------------------------------------------------------------------
+function mnWav1 = fftdiff_(mnWav, P)
+    
+fGpu = isGpu_(mnWav);
+nLoads_gpu = get_set_(P, 'nLoads_gpu', 8);  % GPU load limit    
+
+% [fGpu, nLoads_gpu] = deal(0, 1); %debug
+
+nSamples = size(mnWav,1);
+[nLoad1, nSamples_load1, nSamples_last1] = partition_load_(nSamples, round(nSamples/nLoads_gpu));
+mnWav1 = zeros(size(mnWav), 'like', mnWav);    
+freqLim_ = P.freqLim / (P.sRateHz / 2);
+for iLoad = 1:nLoad1
+    iOffset = (iLoad-1) * nSamples_load1;
+    if iLoad<nLoad1
+        vi1 = (1:nSamples_load1) + iOffset;
+    else
+        vi1 = (1:nSamples_last1) + iOffset;
+    end
+    mnWav1_ = mnWav(vi1,:);
+    if fGpu % use GPU
+        try 
+            mnWav1(vi1,:) = fftdiff__(mnWav1_, freqLim_);
+        catch
+            fGpu = 0;
+        end
+    end
+    if ~fGpu % use CPU 
+        mnWav1(vi1,:) = fftdiff__(gather_(mnWav1_), freqLim_);
+    end
+end %for
+end %func
+
+
+%--------------------------------------------------------------------------
+function mnWav1 = fftdiff__(mnWav, freqLim_)
+% apply fft to diffrentiate
+% mnWav = gather_(mnWav);
+
+n = size(mnWav,1);
+% n1 = round(n/2*freqLim_(1));
+% n2 = round(n/2*diff(freqLim_));
+
+n1 = round(n/2 * freqLim_(2));
+npow2 = 2^nextpow2(n);
+% w = single([linspace(0, 1, n2), linspace(1, 0, n2)])';
+% w = [zeros(n1, 1, 'single'); w; zeros(npow2-2*n1-4*n2, 1, 'single'); -w; zeros(n1, 1, 'single')];
+% w = single(pi*1i) * w;
+w = single(pi*1i) * single([linspace(0, 1, n1), linspace(1, -1, npow2-2*n1), linspace(-1, 0, n1)]');
+mnWav1 = real(ifft(bsxfun(@times, fft(single(mnWav), npow2), w), 'symmetric'));
+mnWav1 = cast(mnWav1(1:n,:), class_(mnWav));
+end %func
+
+
 %--------------------------------------------------------------------------
 function mnWav1 = fread_(fid_bin, dimm_wav, vcDataType)
 % Get around fread bug (matlab) where built-in fread resize doesn't work
@@ -3979,6 +4037,8 @@ function add_menu_(hFig, P)
     'Checked', if_on_off_(P.vcFet_show, {'vpp', 'vmin'}));
 uimenu(mh_proj, 'Label', 'pca', 'Callback', @(h,e)proj_view_(h), ...
     'Checked', if_on_off_(P.vcFet_show, {'pca'}));
+uimenu(mh_proj, 'Label', 'ppca', 'Callback', @(h,e)proj_view_(h), ...
+    'Checked', if_on_off_(P.vcFet_show, {'ppca', 'private pca'}));
 % uimenu(mh_proj, 'Label', 'cov', 'Callback', @(h,e)proj_view_(h), ...
 %     'Checked', if_on_off_(P.vcFet_show, {'cov', 'spacetime'}));
 
@@ -4601,26 +4661,8 @@ function plot_FigTime_(S0)
 % S_clu = S0.S_clu;
 P = S0.P;
 if ~isfield(P, 'vcFet_show'), P.vcFet_show = 'vpp'; end
-% vrFet1 = []; vrTime1=[];
-% iFet = str2double(P.vcFet_show(end));
 [vrFet1, viSpk1] = getFet_clu_(iClu, iSite, S0);
 vrTime1 = double(S0.viTime_spk(viSpk1)) / P.sRateHz;
-        
-% if nargin < 2
-%    % get background feature      
-%    switch lower(P.vcFet_show)
-%        case {'vmin', 'vpp'}
-%             vrFet1 = S0.cvrVpp_site{iSite};   
-%        otherwise
-%            error('not implemented yet'); % list of background spikes
-% %                 vrFet1 = S0.cmrFet_site{iSite}(:,iFet);
-%     end
-%    vrTime1 = S0.cvrTime_site{iSite};
-%    viSpk1 = []; %@TODO: Sclu.cviSpk_site{iSite};
-% else        
-%     [vrFet1, viSpk1] = getFet_clu_(iClu, iSite);
-%     vrTime1 = double(S0.viTime_spk(viSpk1)) / P.sRateHz;
-% end
 
 % label
 switch lower(P.vcFet_show)
@@ -5026,6 +5068,15 @@ function update_plot2_proj_(vrX, vrY)
         if ~isempty(iClu2)  
             [mrMin2, mrMax2] = pca_pc_spk_(viSpk02, viSites0, mrPv1, mrPv2);
         end 
+                
+    case {'ppca', 'private pca'} %channel by channel pca. do it by channel
+        % determine pca vector from cluster 1
+        [mrPv1, mrPv2] = pca_pv_clu_(viSites0, iClu1, iClu2);            
+        [mrMin0, mrMax0] = pca_pc_spk_(viSpk00, viSites0, mrPv1, mrPv2); %getall spikes whose center lies in certain range
+        [mrMin1, mrMax1] = pca_pc_spk_(viSpk01, viSites0, mrPv1, mrPv2); %getall spikes whose center lies in certain range
+        if ~isempty(iClu2)              
+            [mrMin2, mrMax2] = pca_pc_spk_(viSpk02, viSites0, mrPv1, mrPv2);
+        end
         
     otherwise % generic
         [mrMin0, mrMax0] = getFet_spk_(viSpk00, viSites0, S0); %getall spikes whose center lies in certain range
@@ -5112,7 +5163,7 @@ function update_plot2_proj_(vrX, vrY)
 % if nargin<3, fWav_raw_show = P.fWav_raw_show; end
 if nargin<3, S0 = []; end
 if isempty(S0), S0 = get(0, 'UserData'); end
-if nargin<4, fWav_raw_show = 0; end
+if nargin<4, fWav_raw_show = get_set_(S0.P, 'fWav_raw_show', 0); end
 
 % unique exception handling %171201 JJJ
 [viSites1_uniq, ~, viiSites1_uniq] = unique(viSites1);
@@ -5533,11 +5584,26 @@ function rescale_FigWav_(event, S0, P)
 function rescale_FigTime_(event, S0, P)
 % rescale_FigTime_(event, S0, P)
 % rescale_FigTime_(maxAmp, S0, P)
-if nargin<3, P = S0.P; end
+if nargin<2, S0 = []; end
+if nargin<3, P = []; end
+
+if isempty(S0), S0 = get0_(); end
+if isempty(P), P = S0.P; end
+S_clu = S0.S_clu;
 
 [S_fig, maxAmp_prev] = set_fig_maxAmp_('FigTime', event);
 ylim_(S_fig.hAx, [0, 1] * S_fig.maxAmp);
 imrect_set_(S_fig.hRect, [], [0, S_fig.maxAmp]);
+iSite = S_clu.viSite_clu(S0.iCluCopy);
+
+% switch lower(P.vcFet_show)
+%     case {'vpp', 'vmin'} %voltage feature
+%         vcYlabel = sprintf('Site %d (\\mu%s)', iSite, P.vcFet_show);
+%     otherwise %other feature options
+%         vcYlabel = sprintf('Site %d (%s)', iSite, P.vcFet_show);
+% end
+% ylabel(S_fig.hAx, vcYlabel);
+
 end %func
 
 
@@ -5757,7 +5823,7 @@ function keyPressFcn_FigProj_(hFig, event)
 figure_wait_(1);
 switch lower(event.Key)
     case {'uparrow', 'downarrow'}
-        rescale_FigProj_(event, hFig, S_fig);
+        rescale_FigProj_(event, hFig, S_fig, S0);
 
     case {'leftarrow', 'rightarrow'} % change channels
         fPlot = 0;
@@ -5830,9 +5896,10 @@ function keyPressFcn_FigProj_(hFig, event)
 if nargin<2, hFig = []; end
 if nargin<3, S_fig = []; end
 if nargin<4, S0 = []; end
-
 if isempty(hFig) || isempty(S_fig), [hFig, S_fig] = get_fig_cache_('FigProj'); end
 if isempty(S0), S0 = get(0, 'UserData'); end
+P = S0.P;
+
 if isnumeric(event)
     S_fig.maxAmp = event;
 else
@@ -5841,6 +5908,14 @@ function keyPressFcn_FigProj_(hFig, event)
 vhPlot = [S_fig.hPlot0, S_fig.hPlot1, S_fig.hPlot2];
 if isempty(S0.iCluPaste), vhPlot(end) = []; end
 rescaleProj_(vhPlot, S_fig.maxAmp, S0.P);
+switch lower(P.vcFet_show)
+    case {'vpp', 'vmin', 'vmax'}
+        S_fig.vcXLabel = 'Site # (%0.0f \\muV; upper: V_{min}; lower: V_{max})';
+        S_fig.vcYLabel = 'Site # (%0.0f \\muV_{min})';
+    otherwise
+        S_fig.vcXLabel = sprintf('Site # (%%0.0f %s; upper: %s1; lower: %s2)', P.vcFet_show, P.vcFet_show, P.vcFet_show);
+        S_fig.vcYLabel = sprintf('Site # (%%0.0f %s)', P.vcFet_show);    
+end
 xlabel(S_fig.hAx, sprintf(S_fig.vcXLabel, S_fig.maxAmp));   
 ylabel(S_fig.hAx, sprintf(S_fig.vcYLabel, S_fig.maxAmp));  
 if nargout==0
@@ -5913,6 +5988,9 @@ function button_FigWavCor_(xyPos, vcButton)
         mrFet1 = calc_cov_spk_(viSpk1, iSite);
     case {'pca', 'gpca'}
         mrFet1 = pca_pc_spk_(viSpk1, iSite);
+    case {'ppca', 'private pca'}
+        [mrPv1, mrPv2] = pca_pv_clu_(iSite, S0.iCluCopy);
+        mrFet1 = pca_pc_spk_(viSpk1, iSite, mrPv1, mrPv2);
     otherwise
         error('not implemented yet');
 end
@@ -7017,7 +7095,12 @@ function update_FigTime_()
 S_clu = S0.S_clu;
 P = S0.P;
 iClu1 = S0.iCluCopy;
-trWav12 = tnWav2uV_(tnWav_spk_sites_(S_clu.cviSpk_clu{iClu1}, site12, S0), P);
+if ismember(P.vcFet_show, {'pca', 'ppca', 'gpca'})
+    fWav_raw_show = 0;
+else
+    fWav_raw_show = get_set_(P, 'fWav_raw_show', 0);
+end
+trWav12 = tnWav2uV_(tnWav_spk_sites_(S_clu.cviSpk_clu{iClu1}, site12, S0, fWav_raw_show), P);
 if diff(site12) == 0, trWav12(:,2,:) = trWav12(:,1,:); end
 vxPoly = (mrPolyPos([1:end,1],1) - site12_show(1)) * S1.maxAmp;
 vyPoly = (mrPolyPos([1:end,1],2) - site12_show(2)) * S1.maxAmp;
@@ -7049,8 +7132,13 @@ function update_FigTime_()
             vcXlabel = sprintf('Site %d (cov2)', site12(1));
         end  
         
-    case 'pca'
-        [mrAmin12, mrAmax12] = pca_pc_spk_(S_clu.cviSpk_clu{iClu1}, site12);
+    case {'pca', 'ppca', 'gpca'}
+        if strcmpi(P.vcFet_show, 'ppca')
+            [mrPv1, mrPv2] = pca_pv_clu_(site12, S0.iCluCopy, S0.iCluPaste);
+            [mrAmin12, mrAmax12] = pca_pc_spk_(S_clu.cviSpk_clu{iClu1}, site12, mrPv1, mrPv2);
+        else
+            [mrAmin12, mrAmax12] = pca_pc_spk_(S_clu.cviSpk_clu{iClu1}, site12);
+        end
         [mrAmin12, mrAmax12] = multifun_(@(x)abs(x'), mrAmin12, mrAmax12);
         vyPlot = mrAmin12(:,2);
         vcYlabel = sprintf('Site %d (PC1)', site12(2));
@@ -7555,6 +7643,7 @@ function update_spikes_(varargin)
     nClu_pre = size(S_clu.trWav_spk_clu, 3);
     vlClu_update((1:nClu) > nClu_pre) = 1;
     [tmrWav_spk_clu, tmrWav_raw_clu] = deal(S_clu.tmrWav_spk_clu, S_clu.tmrWav_raw_clu);
+    [trWav_spk_clu, trWav_raw_clu] = deal(S_clu.trWav_spk_clu, S_clu.trWav_raw_clu);
     [tmrWav_raw_lo_clu, tmrWav_raw_hi_clu] = deal(S_clu.tmrWav_raw_lo_clu, S_clu.tmrWav_raw_hi_clu);    
 else
     vlClu_update = true(nClu, 1);
@@ -11380,18 +11469,18 @@ function raw_waveform_(hMenu)
 %--------------------------------------------------------------------------
 function [mrPv1, mrPv2] = pca_pv_spk_(viSpk1, viSites1, tnWav_spk1)
 % if viSite not found just set to zero
-
-if nargin<3
-    tnWav_spk1 = permute(tnWav_spk_sites_(viSpk1, viSites1), [1,3,2]);
-end
-nT = size(tnWav_spk1, 1);
-nSites = numel(viSites1);
-mrPv_global = get0_('mrPv_global');
+nSites = numel(viSites1);  
+S0 = get0_();
+mrPv_global = get_(S0, 'mrPv_global'); 
 if ~isempty(mrPv_global)
     % show global pca
     mrPv1 = repmat(mrPv_global(:,1), [1, nSites]);
     mrPv2 = repmat(mrPv_global(:,2), [1, nSites]);
 else
+    if nargin<3
+        tnWav_spk1 = permute(tnWav_spk_sites_(viSpk1, viSites1, S0, 0), [1,3,2]);
+    end
+    nT = size(tnWav_spk1, 1);      
     % show site pca
     [mrPv1, mrPv2] = deal(zeros(nT, nSites, 'single'));
     for iSite1=1:nSites
@@ -11401,6 +11490,36 @@ function raw_waveform_(hMenu)
 end %func
 
 
+%--------------------------------------------------------------------------
+function [mrPv1, mrPv2] = pca_pv_clu_(viSites, iClu1, iClu2)
+% [mrPv1, mrPv2] = pca_pv_clu_(viSites, iClu1): return two prinvec per site from same clu
+% [mrPv1, mrPv2] = pca_pv_clu_(viSites, iClu1, iClu2): return one prinvec per clu per site
+
+if nargin<3, iClu2 = []; end
+S0 = get0_();
+S_clu = S0.S_clu;
+% show site pca
+MAX_SAMPLE = 10000; %for pca
+viSpk1 = subsample_vr_(S_clu.cviSpk_clu{iClu1}, MAX_SAMPLE);
+tnWav_spk1 = permute(tnWav_spk_sites_(viSpk1, viSites, S0, 0), [1,3,2]);
+[nT, nSites] = deal(size(tnWav_spk1, 1), numel(viSites));
+[mrPv1, mrPv2] = deal(zeros(nT, nSites, 'single'));
+
+if isempty(iClu2)
+    for iSite1=1:nSites
+        [mrPv1(:,iSite1), mrPv2(:,iSite1)] = pca_pv_(tnWav_spk1(:,:,iSite1));
+    end %for
+else
+    viSpk2 = subsample_vr_(S_clu.cviSpk_clu{iClu2}, MAX_SAMPLE);    
+    tnWav_spk2 = permute(tnWav_spk_sites_(S_clu.cviSpk_clu{iClu2}, viSites, S0, 0), [1,3,2]);
+    for iSite1=1:nSites
+        mrPv1(:,iSite1) = pca_pv_(tnWav_spk1(:,:,iSite1));
+        mrPv2(:,iSite1) = pca_pv_(tnWav_spk2(:,:,iSite1));
+    end %for
+end
+end %func
+
+
 %--------------------------------------------------------------------------
 function [vrPv1, vrPv2] = pca_pv_(mr1)
 MAX_SAMPLE = 10000; %for pca
@@ -11456,17 +11575,14 @@ function raw_waveform_(hMenu)
 function proj_view_(hMenu)
 P = get0_('P');
 vcFet_show = lower(get(hMenu, 'Label'));
-switch vcFet_show
-    case {'vpp', 'pca', 'cov'}, P.vcFet_show = vcFet_show;
-    otherwise, error('proj_view_:Invalid vcFet_show');
-end
+P.vcFet_show = vcFet_show;
 vhMenu = hMenu.Parent.Children;
 for iMenu=1:numel(vhMenu)
     vhMenu(iMenu).Checked = if_on_off_(vhMenu(iMenu).Label, vcFet_show);
 end
 % auto-scale the view
-set0_(P);
-auto_scale_proj_time_();
+S0 = set0_(P);
+button_CluWav_simulate_(S0.iCluCopy, S0.iCluPaste, S0);
 end %func
 
 
@@ -14062,39 +14178,40 @@ function auto_scale_proj_time_(S0, fPlot)
 if nargin<2, fPlot = 0; end
 
 autoscale_pct = get_set_(S0.P, 'autoscale_pct', 99.5);
-[hFig, S_fig] = get_fig_cache_('FigProj');
-% vrY = [S_fig.hPlot0.YData(:); S_fig.hPlot1.YData(:); S_fig.hPlot2.YData(:)];
-[mrMin0, mrMax0, mrMin1, mrMax1, mrMin2, mrMax2] = fet2proj_(S0, S_fig.viSites_show);
-% mrAmp = [mrMin0, mrMax0, mrMin1, mrMax1, mrMin2, mrMax2];
-% mrAmp = [mrMin1, mrMax1, mrMin2, mrMax2];
-% S_fig.maxAmp = quantile(mrAmp(:), autoscale_pct/100);
+[hFig_proj, S_fig_proj] = get_fig_cache_('FigProj');
+[mrMin0, mrMax0, mrMin1, mrMax1, mrMin2, mrMax2] = fet2proj_(S0, S_fig_proj.viSites_show);
 if isempty(mrMin2) || isempty(mrMax2)
     cmrAmp = {mrMin1, mrMax1};
 else
     cmrAmp = {mrMin1, mrMax1, mrMin2, mrMax2};
 end
-S_fig.maxAmp = max(cellfun(@(x)quantile(x(:), autoscale_pct/100), cmrAmp));
-% S_fig.maxAmp %debug
-set(hFig, 'UserData', S_fig);
+S_fig_proj.maxAmp = max(cellfun(@(x)quantile(x(:), autoscale_pct/100), cmrAmp));
+set(hFig_proj, 'UserData', S_fig_proj);
 
-[hFig, S_fig] = get_fig_cache_('FigTime');
+
+% Update time
+[hFig_time, S_fig_time] = get_fig_cache_('FigTime');
 iSite = S0.S_clu.viSite_clu(S0.iCluCopy);
-[vrFet0, vrTime0] = getFet_site_(iSite, [], S0);    % plot background    
+% [vrFet0, vrTime0] = getFet_site_(iSite, [], S0);    % plot background    
 [vrFet1, vrTime1, vcYlabel, viSpk1] = getFet_site_(iSite, S0.iCluCopy, S0); % plot iCluCopy
 if isempty(S0.iCluPaste)
-    vrFet = [vrFet0(:); vrFet1(:)];
+%     vrFet = [vrFet0(:); vrFet1(:)];
+    cvrFet = {vrFet1};
 else
     [vrFet2, vrTime2, vcYlabel, viSpk2] = getFet_site_(iSite, S0.iCluPaste, S0); % plot iCluCopy
-    vrFet = [vrFet0(:); vrFet1(:); vrFet2(:)];
+%     vrFet = [vrFet0(:); vrFet1(:); vrFet2(:)];
+    cvrFet = {vrFet1, vrFet2};
 end
-S_fig.maxAmp = quantile(vrFet, autoscale_pct/100);
-set(hFig, 'UserData', S_fig);
+% S_fig_time.maxAmp = quantile(vrFet, autoscale_pct/100);
+S_fig_time.maxAmp = max(cellfun(@(x)quantile(x(:), autoscale_pct/100), cvrFet));
+set(hFig_time, 'UserData', S_fig_time);
 
+% plot
 if fPlot
-    keyPressFcn_cell_(get_fig_cache_('FigWav'), {'j', 't'}); 
+    keyPressFcn_cell_(get_fig_cache_('FigWav'), {'j', 't'}, S0); 
 else
-    rescale_FigProj_(S_fig.maxAmp);    
-    rescale_FigTime_(S_fig.maxAmp, S0, S0.P);
+    rescale_FigProj_(S_fig_proj.maxAmp, hFig_proj, S_fig_proj, S0);    
+    rescale_FigTime_(S_fig_time.maxAmp, S0, S0.P);
 end
 end %func
 
@@ -14965,11 +15082,20 @@ function uistack_(h, vc)
 if isempty(S_clu), S_clu = get0_('S_clu'); end
 
 if fText
-    csText_clu = arrayfun(@(i)sprintf('%d (%d)', i, S_clu.vnSpk_clu(i)), 1:S_clu.nClu, 'UniformOutput', 0);
+    csText_clu = arrayfun(@(i)sprintf('%d(%d)', i, S_clu.vnSpk_clu(i)), 1:S_clu.nClu, 'UniformOutput', 0);
 else
     csText_clu = arrayfun(@(i)sprintf('%d', i), 1:S_clu.nClu, 'UniformOutput', 0);
 end
-set(S_fig.hAx, 'Xtick', 1:S_clu.nClu, 'XTickLabel', csText_clu)
+set(S_fig.hAx, 'Xtick', 1:S_clu.nClu, 'XTickLabel', csText_clu, 'FontSize', 8);
+try
+    if fText
+        xtickangle(S_fig.hAx, -20); 
+    else
+        xtickangle(S_fig.hAx, 0); 
+    end
+catch; 
+end
+
 S_fig.fText = fText;
 if nargout==0
     hFig = get_fig_cache_('FigWav');
@@ -17865,8 +17991,8 @@ function plot_aux_rate_(fSelectedUnit)
 % 9/29/17 JJJ: Displaying the version number of the program and what's used. #Tested
 function [vcVer, vcDate, vcVer_used] = jrc_version_(vcFile_prm)
 if nargin<1, vcFile_prm = ''; end
-vcVer = 'v3.2.3';
-vcDate = '1/3/2018';
+vcVer = 'v3.2.5';
+vcDate = '1/8/2018';
 vcVer_used = '';
 if nargout==0
     fprintf('%s (%s) installed\n', vcVer, vcDate);
@@ -18869,6 +18995,9 @@ function make_trial_(vcFile_prm, fImec)
 % 12/20/17 JJJ: support combined recordings and header containing files
 % 12/15/17 JJJ: Load a single channel in memory efficient way
 function vrWav = load_bin_chan_(P, iChan)
+% vrWav = load_bin_chan_(P, 0): return average of all valid channels
+% vrWav = load_bin_chan_(P, iChan): return specific channel
+if nargin<2, iChan = 0; end
 vrWav = [];
 
 % Join multiple recordings if P.csFile_merge is set
@@ -18879,29 +19008,88 @@ function make_trial_(vcFile_prm, fImec)
     for iFile=1:numel(csFile_bin)
         P_ = setfield(P, 'vcFile', csFile_bin{iFile});
         cvrWav{iFile} = load_bin_chan_(P_, iChan);
+        fprintf('.');
     end
     vrWav = cell2mat_(cvrWav);
     return;
 end
 
-if ~exist_file_(P.vcFile), return; end    
-try        
-    nBytes = getBytes_(P.vcFile);
-    if isempty(nBytes), return; end
-    header_offset = get_(P, 'header_offset', 0);
-    bytesPerSample = bytesPerSample_(P.vcDataType);
-    nSamples = floor((nBytes-header_offset) / bytesPerSample / P.nChans);
-
-    fid = fopen(P.vcFile, 'r');     
-    fseek(fid, bytesPerSample * (iChan-1) + header_offset, 'bof');
-    vrWav = fread(fid, [nSamples, 1], ['*', lower(P.vcDataType)], (P.nChans-1) * bytesPerSample);     
-    fclose(fid);
+try       
+    vrWav = load_bin_paged_(P, iChan);    
 catch
+    disperr_();
     return;
 end
 end %func
 
 
+%--------------------------------------------------------------------------
+% 01/08/18: read and reduce
+function mn = load_bin_paged_(P, viChan, nBytes_page)
+% ~35x slower than RAM indexing
+% mn = load_bin_reduce_(P, 0)
+% mn = load_bin_reduce_(P, viChans)
+% mn = load_bin_reduce_(P, viChans)
+LOAD_FACTOR = 5;
+if nargin<3, nBytes_page = []; end
+if isempty(nBytes_page)
+    S = memory();
+    nBytes_page = floor(S.MaxPossibleArrayBytes() / LOAD_FACTOR);
+end
+bytesPerSample = bytesPerSample_(P.vcDataType);
+nSamples_page = floor(nBytes_page / P.nChans / bytesPerSample);
+mn = []; 
+
+% Determine number of samples
+if ~exist_file_(P.vcFile) || isempty(viChan), return; end    
+nBytes = getBytes_(P.vcFile);
+if isempty(nBytes), return; end
+header_offset = get_(P, 'header_offset', 0);
+nSamples = floor((nBytes-header_offset) / bytesPerSample / P.nChans);
+
+% Loading loop
+fid = fopen(P.vcFile, 'r');
+try
+    if header_offset>0, fseek(fid, header_offset, 'bof'); end
+    nPages = ceil(nSamples / nSamples_page);
+    if viChan(1) == 0
+        fMean = 1;
+        viChan = P.viSite2Chan;
+        viChan(P.viSiteZero) = [];    
+    else
+        fMean = 0;
+    end
+    if nPages == 1
+        mn = fread(fid, [P.nChans, nSamples], ['*', lower(P.vcDataType)]);
+        mn = mn(viChan,:);
+        if fMean, mn = cast(mean(mn), P.vcDataType); end
+        mn = mn';
+    else
+        if fMean
+            mn = zeros([nSamples, 1], P.vcDataType);
+        else
+            mn = zeros([nSamples, numel(viChan)], P.vcDataType);
+        end
+        for iPage = 1:nPages
+            if iPage < nPages
+                nSamples_ = nSamples_page;        
+            else
+                nSamples_ = nSamples - nSamples_page * (iPage-1);
+            end
+            vi_ = (1:nSamples_) + (iPage-1) * nSamples_page;    
+            mn_ = fread(fid, [P.nChans, nSamples_], ['*', lower(P.vcDataType)]);
+            mn_ = mn_(viChan,:);
+            if fMean, mn_ = cast(mean(mn_), P.vcDataType); end
+            mn(vi_,:) = mn_';
+        end
+    end
+catch
+    disperr_();
+end
+fclose_(fid);
+end %func
+
+
 %--------------------------------------------------------------------------
 % 12/20/17 JJJ: Export to LFP file
 function import_lfp_(P)
@@ -18948,6 +19136,26 @@ function import_lfp_(P)
 end %func
 
 
+%--------------------------------------------------------------------------
+function export_lfp_(P)
+% export LFP waveform to workspace (ordered by the site numbers)
+
+P.vcFile_lfp = strrep(P.vcFile_prm, '.prm', '.lfp.jrc');
+if ~exist_file_(P.vcFile_lfp)
+    import_lfp_(P)
+end
+mnLfp = load_bin_(P.vcFile_lfp, P.vcDataType);
+nSamples = floor(size(mnLfp,1) / P.nChans);
+mnLfp = reshape(mnLfp(1:P.nChans*nSamples), P.nChans, nSamples)';
+
+mnLfp = mnLfp(:, P.viSite2Chan);
+mrSiteXY = P.mrSiteXY;
+assignWorkspace_(mnLfp, mrSiteXY);
+fprintf('\tmnLfp has nSamples x nSites dimension, sites are ordered from the bottom to top, left to right\n');
+fprintf('\tmrSiteXY (site positions) has nSites x 2 dimension; col.1: x-coord, col.2: y-coord (um)\n');
+end %func
+
+
 %--------------------------------------------------------------------------
 function flag = is_new_imec_(vcFile)
 flag = ~isempty(regexp(lower(vcFile), '.imec.ap.bin$'));
@@ -19129,6 +19337,48 @@ function export_quality_(varargin)
 end
 end %func
 
+
+%--------------------------------------------------------------------------
+function export_chan_(P, vcArg1)
+% export list of channels to a bin file, use fskip?
+if isempty(vcArg1)
+    vcArg1 = inputdlg('Which channel(s) to export (separate by commas or space)', 'Channel', 1, {num2str(P.nChans)});
+    if isempty(vcArg1), return; end
+    vcArg1 = vcArg1{1};
+end
+viChan = str2num(vcArg1);
+if isnan(viChan), fprintf(2, 'Must provide a channel number\n'); return; end
+if any(viChan > P.nChans | viChan < 0)
+    fprintf(2, 'Exceeding nChans (=%d).\n', P.nChans);
+    return; 
+end
+try    
+    vcChan_ = sprintf('%d-', viChan);
+    vcFile_out = strrep(P.vcFile_prm, '.prm', sprintf('_ch%s.jrc', vcChan_(1:end-1)));
+    mn = load_bin_chan_(P, viChan);
+    write_bin_(vcFile_out, mn);    
+    if numel(viChan) == 1
+        eval(sprintf('ch%d = mn;', viChan));
+        eval(sprintf('assignWorkspace_(ch%d);', viChan));
+    else
+        assignWorkspace_(mn);        
+    end
+catch
+    fprintf('Out of memory, exporting individual channels\n'); 
+    for iChan1 = 1:numel(viChan)
+        iChan = viChan(iChan1);
+        vcFile_out = strrep(P.vcFile_prm, '.prm', sprintf('_ch%d.jrc', iChan));
+        fprintf('Loading chan %d(%d/%d) from %s\n\t', iChan, iChan1, numel(viChan), P.vcFile_prm);
+        t1 = tic;
+        vn_ = load_bin_chan_(P, iChan);
+        fprintf('\n\ttook %0.1fs\n', toc(t1));
+        write_bin_(vcFile_out, vn_);
+    end %for
+end
+end %func
+
+
+
 %     case 'o' %overlap waveforms across sites
 %         hFig_temp = figure; hold on;
 %         if isempty(S0.iCluPaste)
diff --git a/prb/imec2aux.prb b/prb/imec2aux.prb
new file mode 100644
index 00000000..384d4c97
--- /dev/null
+++ b/prb/imec2aux.prb
@@ -0,0 +1,26 @@
+% Order of the probe sites in the recording file
+channels = 1 + [40  103 39	104	41	102	38	105	42	101	37	106	43	100	36	107, 44	99	35	108	45	98	34	109	46	97	33	110	47	96	32	111, 48	127	63	112	49	126	62	113	50	125	61	114	51	124	60	115, 52	116	59	123	53	117	58	122	54	118	57	121	55	119	56	120, 8	71	7	72	9	70	6	73	10	74	5	69	11	75	4	68, 12	76	3	67	13	77	2	66	14	78	1	65	15	79	0	64, 16	80	31	95	17	81	30	94	18	82	29	93	19	83	28	92,     20	84	27	91	21	85	26	90	22	86	25	89	23	87	24	88];
+channels(end+1:end+4) = 129:132;
+
+% Site location in micrometers (x and y)
+geometry = zeros(numel(channels), 2);
+geometry(1:2:128,1) = 0;
+geometry(2:2:128,1) = 28;
+geometry(1:2:128,2) = 20*(0:63);
+geometry(2:2:128,2) = geometry(1:2:128,2);
+
+% Reference sites are being excluded
+ref_sites = [1 18 33 50 65 82 97 114];
+channels(ref_sites) = []; 
+geometry(ref_sites,:) = []; 
+
+% Recording contact pad size in micrometers. Height x width
+pad = [12 12];
+
+% Default prm
+maxSite = 6.5;
+um_per_pix = 20;
+nSites_ref = 8;
+
+% Shanks
+% shank = ones(size(channels)); shank(geometry(:,1)>0)=2;
\ No newline at end of file