function [df]=ddiff(f,dim,dn)
% [df]=pddiff[f,dim]
% Second order accurrate perodic diff(f,1,dim).  size(df)==size(f) 

perm = [dim:max(ndims(f),dim) 1:dim-1];
f = permute(f,perm);
        
siz = size(f);

df=ipermute(reshape([f([2:siz(1) 1],:)-f([siz(1) 1:siz(1)-1],:)]/2,siz),perm); 
