MATLAB Examples

## Contents

```function [z,s,exitflag] = smoothn(varargin) ```
```%SMOOTHN Robust spline smoothing for 1-D to N-D data. % SMOOTHN provides a fast, automatized and robust discretized spline % smoothing for data of arbitrary dimension. % % Z = SMOOTHN(Y) automatically smoothes the uniformly-sampled array Y. Y % can be any N-D noisy array (time series, images, 3D data,...). Non % finite data (NaN or Inf) are treated as missing values. % % Z = SMOOTHN(Y,S) smoothes the array Y using the smoothing parameter S. % S must be a real positive scalar. The larger S is, the smoother the % output will be. If the smoothing parameter S is omitted (see previous % option) or empty (i.e. S = []), it is automatically determined by % minimizing the generalized cross-validation (GCV) score. % % Z = SMOOTHN(Y,W) or Z = SMOOTHN(Y,W,S) smoothes Y using a weighting % array W of positive values, that must have the same size as Y. Note % that a nil weight corresponds to a missing value. % % If you want to smooth a vector field or multicomponent data, Y must be % a cell array. For example, if you need to smooth a 3-D vectorial flow % (Vx,Vy,Vz), use Y = {Vx,Vy,Vz}. The output Z is also a cell array which % contains the smoothed components. See examples 5 to 8 below. % % [Z,S] = SMOOTHN(...) also returns the calculated value for the % smoothness parameter S so that you can fine-tune the smoothing % subsequently if needed. % % 1) Robust smoothing % ------------------- % Z = SMOOTHN(...,'robust') carries out a robust smoothing that minimizes % the influence of outlying data. % % An iteration process is used with the 'ROBUST' option, or in the % presence of weighted and/or missing values. Z = SMOOTHN(...,OPTIONS) % smoothes with the termination parameters specified in the structure % OPTIONS. OPTIONS is a structure of optional parameters that change the % default smoothing properties. It must be last input argument. % --- % The structure OPTIONS can contain the following fields: % ----------------- % OPTIONS.TolZ: Termination tolerance on Z (default = 1e-3), % OPTIONS.TolZ must be in ]0,1[ % OPTIONS.MaxIter: Maximum number of iterations allowed % (default = 100) % OPTIONS.Initial: Initial value for the iterative process % (default = original data, Y) % OPTIONS.Weight: Weight function for robust smoothing: % 'bisquare' (default), 'talworth' or 'cauchy' % ----------------- % % [Z,S,EXITFLAG] = SMOOTHN(...) returns a boolean value EXITFLAG that % describes the exit condition of SMOOTHN: % 1 SMOOTHN converged. % 0 Maximum number of iterations was reached. % % 2) Different spacing increments % ------------------------------- % SMOOTHN, by default, assumes that the spacing increments are constant % and equal in all the directions (i.e. dx = dy = dz = ...). If the % increments differ from one direction to the other, use the following % field in OPTIONS: % OPTIONS.Spacing' = [d1 d2 d3...], % where dI represents the spacing between points in the Ith dimension. % % Important note: d1 is the spacing increment for the first % non-singleton dimension (i.e. the vertical direction for matrices). % % 3) References (please refer to the two following papers) % ------------- % 1) Garcia D, Robust smoothing of gridded data in one and higher % dimensions with missing values. Computational Statistics & Data % Analysis, 2010;54:1167-1178. % <a % href="matlab:web('http://www.biomecardio.com/pageshtm/publi/csda10.pdf')">PDF download</a> % 2) Garcia D, A fast all-in-one method for automated post-processing of % PIV data. Exp Fluids, 2011;50:1247-1259. % <a % href="matlab:web('http://www.biomecardio.com/pageshtm/publi/media11.pdf')">PDF download</a> % % Examples: % -------- % %--- Example #1: smooth a curve --- % x = linspace(0,100,2^8); % y = cos(x/10)+(x/50).^2 + randn(size(x))/10; % y([70 75 80]) = [5.5 5 6]; % z = smoothn(y); % Regular smoothing % zr = smoothn(y,'robust'); % Robust smoothing % subplot(121), plot(x,y,'r.',x,z,'k','LineWidth',2) % axis square, title('Regular smoothing') % subplot(122), plot(x,y,'r.',x,zr,'k','LineWidth',2) % axis square, title('Robust smoothing') % % %--- Example #2: smooth a surface --- % xp = 0:.02:1; % [x,y] = meshgrid(xp); % f = exp(x+y) + sin((x-2*y)*3); % fn = f + randn(size(f))*0.5; % fs = smoothn(fn); % subplot(121), surf(xp,xp,fn), zlim([0 8]), axis square % subplot(122), surf(xp,xp,fs), zlim([0 8]), axis square % % %--- Example #3: smooth an image with missing data --- % n = 256; % y0 = peaks(n); % y = y0 + randn(size(y0))*2; % I = randperm(n^2); % y(I(1:n^2*0.5)) = NaN; % lose 1/2 of data % y(40:90,140:190) = NaN; % create a hole % z = smoothn(y); % smooth data % subplot(2,2,1:2), imagesc(y), axis equal off % title('Noisy corrupt data') % subplot(223), imagesc(z), axis equal off % title('Recovered data ...') % subplot(224), imagesc(y0), axis equal off % title('... compared with the actual data') % % %--- Example #4: smooth volumetric data --- % [x,y,z] = meshgrid(-2:.2:2); % xslice = [-0.8,1]; yslice = 2; zslice = [-2,0]; % vn = x.*exp(-x.^2-y.^2-z.^2) + randn(size(x))*0.06; % subplot(121), slice(x,y,z,vn,xslice,yslice,zslice,'cubic') % title('Noisy data') % v = smoothn(vn); % subplot(122), slice(x,y,z,v,xslice,yslice,zslice,'cubic') % title('Smoothed data') % % %--- Example #5: smooth a cardioid --- % t = linspace(0,2*pi,1000); % x = 2*cos(t).*(1-cos(t)) + randn(size(t))*0.1; % y = 2*sin(t).*(1-cos(t)) + randn(size(t))*0.1; % z = smoothn({x,y}); % plot(x,y,'r.',z{1},z{2},'k','linewidth',2) % axis equal tight % % %--- Example #6: smooth a 3-D parametric curve --- % t = linspace(0,6*pi,1000); % x = sin(t) + 0.1*randn(size(t)); % y = cos(t) + 0.1*randn(size(t)); % z = t + 0.1*randn(size(t)); % u = smoothn({x,y,z}); % plot3(x,y,z,'r.',u{1},u{2},u{3},'k','linewidth',2) % axis tight square % % %--- Example #7: smooth a 2-D vector field --- % [x,y] = meshgrid(linspace(0,1,24)); % Vx = cos(2*pi*x+pi/2).*cos(2*pi*y); % Vy = sin(2*pi*x+pi/2).*sin(2*pi*y); % Vx = Vx + sqrt(0.05)*randn(24,24); % adding Gaussian noise % Vy = Vy + sqrt(0.05)*randn(24,24); % adding Gaussian noise % I = randperm(numel(Vx)); % Vx(I(1:30)) = (rand(30,1)-0.5)*5; % adding outliers % Vy(I(1:30)) = (rand(30,1)-0.5)*5; % adding outliers % Vx(I(31:60)) = NaN; % missing values % Vy(I(31:60)) = NaN; % missing values % Vs = smoothn({Vx,Vy},'robust'); % automatic smoothing % subplot(121), quiver(x,y,Vx,Vy,2.5), axis square % title('Noisy velocity field') % subplot(122), quiver(x,y,Vs{1},Vs{2}), axis square % title('Smoothed velocity field') % % %--- Example #8: smooth a 3-D vector field --- % load wind % original 3-D flow % % -- Create noisy data % % Gaussian noise % un = u + randn(size(u))*8; % vn = v + randn(size(v))*8; % wn = w + randn(size(w))*8; % % Add some outliers % I = randperm(numel(u)) < numel(u)*0.025; % un(I) = (rand(1,nnz(I))-0.5)*200; % vn(I) = (rand(1,nnz(I))-0.5)*200; % wn(I) = (rand(1,nnz(I))-0.5)*200; % % -- Visualize the noisy flow (see CONEPLOT documentation) % figure, title('Noisy 3-D vectorial flow') % xmin = min(x(:)); xmax = max(x(:)); % ymin = min(y(:)); ymax = max(y(:)); % zmin = min(z(:)); % daspect([2,2,1]) % xrange = linspace(xmin,xmax,8); % yrange = linspace(ymin,ymax,8); % zrange = 3:4:15; % [cx cy cz] = meshgrid(xrange,yrange,zrange); % k = 0.1; % hcones = coneplot(x,y,z,un*k,vn*k,wn*k,cx,cy,cz,0); % set(hcones,'FaceColor','red','EdgeColor','none') % hold on % wind_speed = sqrt(un.^2 + vn.^2 + wn.^2); % hsurfaces = slice(x,y,z,wind_speed,[xmin,xmax],ymax,zmin); % set(hsurfaces,'FaceColor','interp','EdgeColor','none') % hold off % axis tight; view(30,40); axis off % camproj perspective; camzoom(1.5) % camlight right; lighting phong % set(hsurfaces,'AmbientStrength',.6) % set(hcones,'DiffuseStrength',.8) % % -- Smooth the noisy flow % Vs = smoothn({un,vn,wn},'robust'); % % -- Display the result % figure, title('3-D flow smoothed by SMOOTHN') % daspect([2,2,1]) % hcones = coneplot(x,y,z,Vs{1}*k,Vs{2}*k,Vs{3}*k,cx,cy,cz,0); % set(hcones,'FaceColor','red','EdgeColor','none') % hold on % wind_speed = sqrt(Vs{1}.^2 + Vs{2}.^2 + Vs{3}.^2); % hsurfaces = slice(x,y,z,wind_speed,[xmin,xmax],ymax,zmin); % set(hsurfaces,'FaceColor','interp','EdgeColor','none') % hold off % axis tight; view(30,40); axis off % camproj perspective; camzoom(1.5) % camlight right; lighting phong % set(hsurfaces,'AmbientStrength',.6) % set(hcones,'DiffuseStrength',.8) % % See also SMOOTH1Q, DCTN, IDCTN. % % -- Damien Garcia -- 2009/03, revised 2014/02 % website: <a % href="matlab:web('http://www.biomecardio.com')">www.BiomeCardio.com</a> ```

## Check input arguments

```error(nargchk(1,5,nargin)); OPTIONS = struct; NArgIn = nargin; if isstruct(varargin{end}) % SMOOTHN(...,OPTIONS) OPTIONS = varargin{end}; NArgIn = NArgIn-1; end idx = cellfun(@ischar,varargin); if any(idx) % SMOOTHN(...,'robust',...) assert(sum(idx)==1 && strcmpi(varargin{idx},'robust'),... 'Wrong syntax. Read the help text for instructions.') isrobust = true; assert(find(idx)==NArgIn,... 'Wrong syntax. Read the help text for instructions.') NArgIn = NArgIn-1; else isrobust = false; end ```
```Error using smoothn (line 223) Not enough input arguments. ```

## Test & prepare the variables

```%--- % y = array to be smoothed y = varargin{1}; if ~iscell(y), y = {y}; end sizy = size(y{1}); ny = numel(y); % number of y components for i = 1:ny assert(isequal(sizy,size(y{i})),... 'Data arrays in Y must have the same size.') y{i} = double(y{i}); end noe = prod(sizy); % number of elements if noe==1, z = y{1}; s = []; exitflag = true; return, end %--- % Smoothness parameter and weights W = ones(sizy); s = []; if NArgIn==2 if isempty(varargin{2}) || isscalar(varargin{2}) % smoothn(y,s) s = varargin{2}; % smoothness parameter else % smoothn(y,W) W = varargin{2}; % weight array end elseif NArgIn==3 % smoothn(y,W,s) W = varargin{2}; % weight array s = varargin{3}; % smoothness parameter end assert(isnumeric(W),'W must be a numeric array') assert(isnumeric(s),'S must be a numeric scalar') assert(isequal(size(W),sizy),... 'Arrays for data and weights (Y and W) must have same size.') assert(isempty(s) || (isscalar(s) && s>0),... 'The smoothing parameter S must be a scalar >0') %--- % Field names in the structure OPTIONS OptionNames = fieldnames(OPTIONS); idx = ismember(OptionNames,... {'TolZ','MaxIter','Initial','Weight','Spacing','Order'}); assert(all(idx),... ['''' OptionNames{find(~idx,1)} ''' is not a valid field name for OPTIONS.']) %--- % "Maximal number of iterations" criterion if ~ismember('MaxIter',OptionNames) OPTIONS.MaxIter = 100; % default value for MaxIter end assert(isnumeric(OPTIONS.MaxIter) && isscalar(OPTIONS.MaxIter) &&... OPTIONS.MaxIter>=1 && OPTIONS.MaxIter==round(OPTIONS.MaxIter),... 'OPTIONS.MaxIter must be an integer >=1') %--- % "Tolerance on smoothed output" criterion if ~ismember('TolZ',OptionNames) OPTIONS.TolZ = 1e-3; % default value for TolZ end assert(isnumeric(OPTIONS.TolZ) && isscalar(OPTIONS.TolZ) &&... OPTIONS.TolZ>0 && OPTIONS.TolZ<1,'OPTIONS.TolZ must be in ]0,1[') %--- % "Initial Guess" criterion if ~ismember('Initial',OptionNames) isinitial = false; else isinitial = true; z0 = OPTIONS.Initial; if ~iscell(z0), z0 = {z0}; end nz0 = numel(z0); % number of y components assert(nz0==ny,... 'OPTIONS.Initial must contain a valid initial guess for Z') for i = 1:nz0 assert(isequal(sizy,size(z0{i})),... 'OPTIONS.Initial must contain a valid initial guess for Z') z0{i} = double(z0{i}); end end %--- % "Weight function" criterion (for robust smoothing) if ~ismember('Weight',OptionNames) OPTIONS.Weight = 'bisquare'; % default weighting function else assert(ischar(OPTIONS.Weight),... 'A valid weight function (OPTIONS.Weight) must be chosen') assert(any(ismember(lower(OPTIONS.Weight),... {'bisquare','talworth','cauchy'})),... 'The weight function must be ''bisquare'', ''cauchy'' or '' talworth''.') end %--- % "Order" criterion (by default m = 2) % Note: m = 0 is of course not recommended! if ~ismember('Order',OptionNames) m = 2; % order else m = OPTIONS.Order; if ~ismember(m,0:2); error('MATLAB:smoothn:IncorrectOrder',... 'The order (OPTIONS.order) must be 0, 1 or 2.') end end %--- % "Spacing" criterion d = ndims(y{1}); if ~ismember('Spacing',OptionNames) dI = ones(1,d); % spacing increment else dI = OPTIONS.Spacing; assert(isnumeric(dI) && length(dI)==d,... 'A valid spacing (OPTIONS.Spacing) must be chosen') end dI = dI/max(dI); %--- % Weights. Zero weights are assigned to not finite values (Inf or NaN), % (Inf/NaN values = missing data). IsFinite = isfinite(y{1}); for i = 2:ny, IsFinite = IsFinite & isfinite(y{i}); end nof = nnz(IsFinite); % number of finite elements W = W.*IsFinite; assert(all(W(:)>=0),'Weights must all be >=0') W = W/max(W(:)); %--- % Weighted or missing data? isweighted = any(W(:)<1); %--- % Automatic smoothing? isauto = isempty(s); ```

## Create the Lambda tensor

```%--- % Lambda contains the eingenvalues of the difference matrix used in this % penalized least squares process (see CSDA paper for details) d = ndims(y{1}); Lambda = zeros(sizy); for i = 1:d siz0 = ones(1,d); siz0(i) = sizy(i); Lambda = bsxfun(@plus,Lambda,... (2-2*cos(pi*(reshape(1:sizy(i),siz0)-1)/sizy(i)))/dI(i)^2); end if ~isauto, Gamma = 1./(1+s*Lambda.^m); end ```

## Upper and lower bound for the smoothness parameter

The average leverage (h) is by definition in [0 1]. Weak smoothing occurs if h is close to 1, while over-smoothing appears when h is near 0. Upper and lower bounds for h are given to avoid under- or over-smoothing. See equation relating h to the smoothness parameter for m = 2 (Equation #12 in the referenced CSDA paper).

```N = sum(sizy~=1); % tensor rank of the y-array hMin = 1e-6; hMax = 0.99; if m==0 % Not recommended. For mathematical purpose only. sMinBnd = 1/hMax^(1/N)-1; sMaxBnd = 1/hMin^(1/N)-1; elseif m==1 sMinBnd = (1/hMax^(2/N)-1)/4; sMaxBnd = (1/hMin^(2/N)-1)/4; elseif m==2 sMinBnd = (((1+sqrt(1+8*hMax^(2/N)))/4/hMax^(2/N))^2-1)/16; sMaxBnd = (((1+sqrt(1+8*hMin^(2/N)))/4/hMin^(2/N))^2-1)/16; end ```

## Initialize before iterating

```%--- Wtot = W; %--- Initial conditions for z if isweighted %--- With weighted/missing data % An initial guess is provided to ensure faster convergence. For that % purpose, a nearest neighbor interpolation followed by a coarse % smoothing are performed. %--- if isinitial % an initial guess (z0) has been already given z = z0; else z = InitialGuess(y,IsFinite); end else z = cell(size(y)); for i = 1:ny, z{i} = zeros(sizy); end end %--- z0 = z; for i = 1:ny y{i}(~IsFinite) = 0; % arbitrary values for missing y-data end %--- tol = 1; RobustIterativeProcess = true; RobustStep = 1; nit = 0; DCTy = cell(1,ny); vec = @(x) x(:); %--- Error on p. Smoothness parameter s = 10^p errp = 0.1; opt = optimset('TolX',errp); %--- Relaxation factor RF: to speedup convergence RF = 1 + 0.75*isweighted; ```

## Main iterative process

```%--- while RobustIterativeProcess %--- "amount" of weights (see the function GCVscore) aow = sum(Wtot(:))/noe; % 0 < aow <= 1 %--- while tol>OPTIONS.TolZ && nit<OPTIONS.MaxIter nit = nit+1; for i = 1:ny DCTy{i} = dctn(Wtot.*(y{i}-z{i})+z{i}); end if isauto && ~rem(log2(nit),1) %--- % The generalized cross-validation (GCV) method is used. % We seek the smoothing parameter S that minimizes the GCV % score i.e. S = Argmin(GCVscore). % Because this process is time-consuming, it is performed from % time to time (when the step number - nit - is a power of 2) %--- fminbnd(@gcv,log10(sMinBnd),log10(sMaxBnd),opt); end for i = 1:ny z{i} = RF*idctn(Gamma.*DCTy{i}) + (1-RF)*z{i}; end % if no weighted/missing data => tol=0 (no iteration) tol = isweighted*norm(vec([z0{:}]-[z{:}]))/norm(vec([z{:}])); z0 = z; % re-initialization end exitflag = nit<OPTIONS.MaxIter; if isrobust %-- Robust Smoothing: iteratively re-weighted process %--- average leverage h = 1; for k = 1:N if m==0 % not recommended - only for numerical purpose h0 = 1/(1+s/dI(k)^(2^m)); elseif m==1 h0 = 1/sqrt(1+4*s/dI(k)^(2^m)); elseif m==2 h0 = sqrt(1+16*s/dI(k)^(2^m)); h0 = sqrt(1+h0)/sqrt(2)/h0; else error('m must be 0, 1 or 2.') end h = h*h0; end %--- take robust weights into account Wtot = W.*RobustWeights(y,z,IsFinite,h,OPTIONS.Weight); %--- re-initialize for another iterative weighted process isweighted = true; tol = 1; nit = 0; %--- RobustStep = RobustStep+1; RobustIterativeProcess = RobustStep<4; % 3 robust steps are enough. else RobustIterativeProcess = false; % stop the whole process end end ```

## Warning messages

```%--- if isauto if abs(log10(s)-log10(sMinBnd))<errp warning('MATLAB:smoothn:SLowerBound',... ['S = ' num2str(s,'%.3e') ': the lower bound for S ',... 'has been reached. Put S as an input variable if required.']) elseif abs(log10(s)-log10(sMaxBnd))<errp warning('MATLAB:smoothn:SUpperBound',... ['S = ' num2str(s,'%.3e') ': the upper bound for S ',... 'has been reached. Put S as an input variable if required.']) end end if nargout<3 && ~exitflag warning('MATLAB:smoothn:MaxIter',... ['Maximum number of iterations (' int2str(OPTIONS.MaxIter) ') has been exceeded. ',... 'Increase MaxIter option (OPTIONS.MaxIter) or decrease TolZ (OPTIONS.TolZ) value.']) end if numel(z)==1, z = z{:}; end ```

## GCV score

```%--- function GCVscore = gcv(p) % Search the smoothing parameter s that minimizes the GCV score %--- s = 10^p; Gamma = 1./(1+s*Lambda.^m); %--- RSS = Residual sum-of-squares RSS = 0; if aow>0.9 % aow = 1 means that all of the data are equally weighted % very much faster: does not require any inverse DCT for kk = 1:ny RSS = RSS + norm(DCTy{kk}(:).*(Gamma(:)-1))^2; end else % take account of the weights to calculate RSS: for kk = 1:ny yhat = idctn(Gamma.*DCTy{kk}); RSS = RSS + norm(sqrt(Wtot(IsFinite)).*... (y{kk}(IsFinite)-yhat(IsFinite)))^2; end end %--- TrH = sum(Gamma(:)); GCVscore = RSS/nof/(1-TrH/noe)^2; end ```
```end function W = RobustWeights(y,z,I,h,wstr) % One seeks the weights for robust smoothing... ABS = @(x) sqrt(sum(abs(x).^2,2)); r = cellfun(@minus,y,z,'UniformOutput',0); % residuals r = cellfun(@(x) x(:),r,'UniformOutput',0); rI = cell2mat(cellfun(@(x) x(I),r,'UniformOutput',0)); MMED = median(rI); % marginal median AD = ABS(bsxfun(@minus,rI,MMED)); % absolute deviation MAD = median(AD); % median absolute deviation %-- Studentized residuals u = ABS(cell2mat(r))/(1.4826*MAD)/sqrt(1-h); u = reshape(u,size(I)); switch lower(wstr) case 'cauchy' c = 2.385; W = 1./(1+(u/c).^2); % Cauchy weights case 'talworth' c = 2.795; W = u<c; % Talworth weights case 'bisquare' c = 4.685; W = (1-(u/c).^2).^2.*((u/c)<1); % bisquare weights otherwise error('MATLAB:smoothn:IncorrectWeights',... 'A valid weighting function must be chosen') end W(isnan(W)) = 0; % NOTE: % ---- % The RobustWeights subfunction looks complicated since we work with cell % arrays. For better clarity, here is how it would look like without the % use of cells. Much more readable, isn't it? % % function W = RobustWeights(y,z,I,h) % % weights for robust smoothing. % r = y-z; % residuals % MAD = median(abs(r(I)-median(r(I)))); % median absolute deviation % u = abs(r/(1.4826*MAD)/sqrt(1-h)); % studentized residuals % c = 4.685; W = (1-(u/c).^2).^2.*((u/c)<1); % bisquare weights % W(isnan(W)) = 0; % end end ```

## Initial Guess with weighted/missing data

```function z = InitialGuess(y,I) ny = numel(y); %-- nearest neighbor interpolation (in case of missing values) if any(~I(:)) z = cell(size(y)); if license('test','image_toolbox') for i = 1:ny [z{i},L] = bwdist(I); z{i} = y{i}; z{i}(~I) = y{i}(L(~I)); end else % If BWDIST does not exist, NaN values are all replaced with the % same scalar. The initial guess is not optimal and a warning % message thus appears. z = y; for i = 1:ny z{i}(~I) = mean(y{i}(I)); end warning('MATLAB:smoothn:InitialGuess',... ['BWDIST (Image Processing Toolbox) does not exist. ',... 'The initial guess may not be optimal; additional',... ' iterations can thus be required to ensure complete',... ' convergence. Increase ''MaxIter'' criterion if necessary.']) end else z = y; end %-- coarse fast smoothing using one-tenth of the DCT coefficients siz = size(z{1}); z = cellfun(@(x) dctn(x),z,'UniformOutput',0); for k = 1:ndims(z{1}) for i = 1:ny z{i}(ceil(siz(k)/10)+1:end,:) = 0; z{i} = reshape(z{i},circshift(siz,[0 1-k])); z{i} = shiftdim(z{i},1); end end z = cellfun(@(x) idctn(x),z,'UniformOutput',0); end ```

## DCTN

```function y = dctn(y) %DCTN N-D discrete cosine transform. % Y = DCTN(X) returns the discrete cosine transform of X. The array Y is % the same size as X and contains the discrete cosine transform % coefficients. This transform can be inverted using IDCTN. % % Reference % --------- % Narasimha M. et al, On the computation of the discrete cosine % transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936. % % Example % ------- % RGB = imread('autumn.tif'); % I = rgb2gray(RGB); % J = dctn(I); % imshow(log(abs(J)),[]), colormap(jet), colorbar % % The commands below set values less than magnitude 10 in the DCT matrix % to zero, then reconstruct the image using the inverse DCT. % % J(abs(J)<10) = 0; % K = idctn(J); % figure, imshow(I) % figure, imshow(K,[0 255]) % % -- Damien Garcia -- 2008/06, revised 2011/11 % -- www.BiomeCardio.com -- y = double(y); sizy = size(y); y = squeeze(y); dimy = ndims(y); % Some modifications are required if Y is a vector if isvector(y) dimy = 1; if size(y,1)==1, y = y.'; end end % Weighting vectors w = cell(1,dimy); for dim = 1:dimy n = (dimy==1)*numel(y) + (dimy>1)*sizy(dim); w{dim} = exp(1i*(0:n-1)'*pi/2/n); end % --- DCT algorithm --- if ~isreal(y) y = complex(dctn(real(y)),dctn(imag(y))); else for dim = 1:dimy siz = size(y); n = siz(1); y = y([1:2:n 2*floor(n/2):-2:2],:); y = reshape(y,n,[]); y = y*sqrt(2*n); y = ifft(y,[],1); y = bsxfun(@times,y,w{dim}); y = real(y); y(1,:) = y(1,:)/sqrt(2); y = reshape(y,siz); y = shiftdim(y,1); end end y = reshape(y,sizy); end ```

## IDCTN

```function y = idctn(y) %IDCTN N-D inverse discrete cosine transform. % X = IDCTN(Y) inverts the N-D DCT transform, returning the original % array if Y was obtained using Y = DCTN(X). % % Reference % --------- % Narasimha M. et al, On the computation of the discrete cosine % transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936. % % Example % ------- % RGB = imread('autumn.tif'); % I = rgb2gray(RGB); % J = dctn(I); % imshow(log(abs(J)),[]), colormap(jet), colorbar % % The commands below set values less than magnitude 10 in the DCT matrix % to zero, then reconstruct the image using the inverse DCT. % % J(abs(J)<10) = 0; % K = idctn(J); % figure, imshow(I) % figure, imshow(K,[0 255]) % % See also DCTN, IDSTN, IDCT, IDCT2, IDCT3. % % -- Damien Garcia -- 2009/04, revised 2011/11 % -- www.BiomeCardio.com -- y = double(y); sizy = size(y); y = squeeze(y); dimy = ndims(y); % Some modifications are required if Y is a vector if isvector(y) dimy = 1; if size(y,1)==1 y = y.'; end end % Weighing vectors w = cell(1,dimy); for dim = 1:dimy n = (dimy==1)*numel(y) + (dimy>1)*sizy(dim); w{dim} = exp(1i*(0:n-1)'*pi/2/n); end % --- IDCT algorithm --- if ~isreal(y) y = complex(idctn(real(y)),idctn(imag(y))); else for dim = 1:dimy siz = size(y); n = siz(1); y = reshape(y,n,[]); y = bsxfun(@times,y,w{dim}); y(1,:) = y(1,:)/sqrt(2); y = ifft(y,[],1); y = real(y*sqrt(2*n)); I = (1:n)*0.5+0.5; I(2:2:end) = n-I(1:2:end-1)+1; y = y(I,:); y = reshape(y,siz); y = shiftdim(y,1); end end y = reshape(y,sizy); end ```