nt_dsample

[y,yy,yyy]=nt_dsample(x,dsr,method) - downsample by averaging neighboring samples

  y: downsampled data
  yy: instantaneous power of residual
  yyy: original minus downsampled
 
  x: data to downsample (2 or 3D)
  dsr: downsampling ratio (must be integer)
  method: 'smooth' or 'resample' [default: 'smooth']

 Downsampling is performed along columns.  If size(x,1) is not multiple of
 factor, it is truncated.
 
 If method=smooth, data are lowpass filtered by convolution with a square window.
 This ensures minimal temporal distortion of the waveform. However it does not
 strongly attenuate frequency components beyond the Nyquist frequency, so
 it is not optimal from a frequency-domain point of view. 
 
 If method=resample, downsampling is performed with a version of
 resample() modified so that the lowpass corner is 0.8 times Nyquist.
 
 If output is not assigned, plot downsampled signal and residual power.

 NoiseTools

0001 function [y,yy,yyy]=nt_dsample(x,dsr,method)
0002 %[y,yy,yyy]=nt_dsample(x,dsr,method) - downsample by averaging neighboring samples
0003 %
0004 %  y: downsampled data
0005 %  yy: instantaneous power of residual
0006 %  yyy: original minus downsampled
0007 %
0008 %  x: data to downsample (2 or 3D)
0009 %  dsr: downsampling ratio (must be integer)
0010 %  method: 'smooth' or 'resample' [default: 'smooth']
0011 %
0012 % Downsampling is performed along columns.  If size(x,1) is not multiple of
0013 % factor, it is truncated.
0014 %
0015 % If method=smooth, data are lowpass filtered by convolution with a square window.
0016 % This ensures minimal temporal distortion of the waveform. However it does not
0017 % strongly attenuate frequency components beyond the Nyquist frequency, so
0018 % it is not optimal from a frequency-domain point of view.
0019 %
0020 % If method=resample, downsampling is performed with a version of
0021 % resample() modified so that the lowpass corner is 0.8 times Nyquist.
0022 %
0023 % If output is not assigned, plot downsampled signal and residual power.
0024 %
0025 % NoiseTools
0026 
0027 if nargin<3||isempty(method); method='smooth'; end
0028 if nargin<2; error('!'); end
0029 if dsr==1; 
0030     y=x;
0031     yy=zeros(size(x));
0032     yyy=x;
0033     return; 
0034 end
0035 if dsr ~= round(dsr); error('factor must be integer'); end
0036 
0037 if iscell(x)
0038     switch nargout
0039         case 0
0040             for iTrial=1:numel(x)
0041                 nt_dsample(x{iTrial},dsr,method);
0042                 pause
0043             end
0044         case 1
0045            y={};
0046            for iTrial=1:size(x)
0047                y{iTrial}=nt_dsample(x{iTrial}.factor,method);
0048            end
0049         case 2
0050            y={};yy={};
0051            for iTrial=1:size(x)
0052                [y{iTrial},yy{iTrial}]=nt_dsample(x{iTrial}.factor,method);
0053            end
0054         case 3
0055            y={};yy={};
0056            for iTrial=1:size(x)
0057                [y{iTrial},yy{iTrial},yyy{iTrial}]=nt_dsample(x{iTrial}.factor,method);
0058            end
0059         otherwise
0060             error('!');
0061     end
0062     return
0063 end
0064 
0065 if ndims(x)>2
0066     d=size(x);
0067     switch nargout
0068         case 0
0069             for iTrial=1:size(x(:,:,:),3)
0070                 nt_dsample(x(:,:,iTrial),dsr,method);
0071                 pause
0072             end
0073         case 1
0074             x=x(:,:);
0075             y=nt_dsample(x,dsr);
0076             y=reshape(y,[size(y,1),d(2:end)]);
0077         case 2
0078             x=x(:,:);
0079             [y,yy]=nt_dsample(x,dsr);
0080             y=reshape(y,[size(y,1),d(2:end)]);
0081             yy=reshape(yy,[size(yy,1),d(2:end)]);
0082         case 3
0083             x=x(:,:);
0084             [y,yy,yyy]=nt_dsample(x,dsr);
0085             y=reshape(y,[size(y,1),d(2:end)]);
0086             yy=reshape(yy,[size(yy,1),d(2:end)]);
0087             yyy=reshape(yyy,[size(yyy,1),d(2:end)]);
0088         otherwise
0089             error('!'); 
0090     end
0091     return
0092 end
0093 
0094 if nargout==0
0095     [y,yy]=nt_dsample(x,dsr,method);
0096     figure(100); clf
0097     subplot 211; 
0098     plot(y);  xlabel('sample'); title('downsampled');
0099     subplot 212;
0100     plot(yy ./ (yy+y.^2)); xlabel('sample'); title(['proportion power lost (overall: ',num2str(mean(yy(:)/(mean(yy(:))+mean(y(:).^2)))), ')']); 
0101     ylim([0 1.1]);
0102     if size(x,2)>1;
0103         hold on; plot(mean(yy ./ (yy+y.^2),2), ':k','linewidth', 2);
0104         figure(101); clf;
0105         plot(mean(yy(:,:)./(mean(yy(:,:))+mean(y(:,:).^2))), '.-k')
0106         xlabel('channel'); title('proportion power lost per channel');
0107     end
0108     return
0109 end
0110 
0111 
0112 [m,n]=size(x);
0113 a=floor(m/dsr);
0114 x=x(1:a*dsr,:);
0115 x=reshape(x,[dsr,a,n]);
0116 y=mean(x,1);
0117 
0118 if nargout>=2
0119     yy=mean(bsxfun(@minus,x,y).^2,1);
0120     yy=shiftdim(yy,1);
0121 end
0122 if nargout==3
0123     yyy=reshape(bsxfun(@minus,x,y),[m,n]);
0124 end
0125     
0126 y=shiftdim(y,1);
0127 
0128