v4_3deg/matlab/xpolate.m

function y=xpolate(x,m,n)

% xpolate(x,m,n)
%             extra/interpolate using tony's 2-D fft method
%             replaces nans in a 2-D array with sensible values
%
%       x     is the 2-D array to be filled
%       m,n   optional arguments to specify size of 2-D fft
%             if n is not specified, n=m;
%             if m and n are not specified, [m n]=size(x);
%
%             the algorithm truncates progressively fewer
%             high frequency components in order to achieve
%             smooth transitions to areas with no data

[M N]=size(x);

if     nargin==2, n=m;
elseif nargin==1, m=M; n=N;
end

y=nan*ones(m,n); y(1:M,1:N)=x;

[X,Y]=meshgrid(1:n,1:m); X=X-1; Y=Y-1; % distance from top left corner
ix=find(isnan(y));                   % location of nans
y(ix)=mmean(x)*ones(size(ix));       % replace nans by mean value
mc=fix(m/2)+1; nc=fix(n/2)+1;        % center value
w=zeros(m,n);
maxx=mmax(x); minx=mmin(x);

for i=1:(max(mc,nc)-1)
  
  % take 2-d fft of y
  
  f=fftshift(fft2(y));
  
  % truncate high frequency

  mx=max((mc-i),1):min((mc+i),m);
  nx=max((nc-i),1):min((nc+i),n);
  w(mx,nx)=ones(length(mx),length(nx));
  f=f.*w;
  
  % map back onto y
  
  f=real(ifft2(ifftshift(f)));       % revert to space domain
  y(ix)=f(ix);                       % replace missing values
  y(find(y>maxx))=maxx;
  y(find(y<minx))=minx;
  
end

y=y(1:M,1:N);
1	function y=xpolate(x,m,n)
2
3	% xpolate(x,m,n)
4	% extra/interpolate using tony's 2-D fft method
5	% replaces nans in a 2-D array with sensible values
6	%
7	% x is the 2-D array to be filled
8	% m,n optional arguments to specify size of 2-D fft
9	% if n is not specified, n=m;
10	% if m and n are not specified, [m n]=size(x);
11	%
12	% the algorithm truncates progressively fewer
13	% high frequency components in order to achieve
14	% smooth transitions to areas with no data
15
16	[M N]=size(x);
17
18	if nargin==2, n=m;
19	elseif nargin==1, m=M; n=N;
20	end
21
22	y=nan*ones(m,n); y(1:M,1:N)=x;
23
24	[X,Y]=meshgrid(1:n,1:m); X=X-1; Y=Y-1; % distance from top left corner
25	ix=find(isnan(y)); % location of nans
26	y(ix)=mmean(x)*ones(size(ix)); % replace nans by mean value
27	mc=fix(m/2)+1; nc=fix(n/2)+1; % center value
28	w=zeros(m,n);
29	maxx=mmax(x); minx=mmin(x);
30
31	for i=1:(max(mc,nc)-1)
32
33	% take 2-d fft of y
34
35	f=fftshift(fft2(y));
36
37	% truncate high frequency
38
39	mx=max((mc-i),1):min((mc+i),m);
40	nx=max((nc-i),1):min((nc+i),n);
41	w(mx,nx)=ones(length(mx),length(nx));
42	f=f.*w;
43
44	% map back onto y
45
46	f=real(ifft2(ifftshift(f))); % revert to space domain
47	y(ix)=f(ix); % replace missing values
48	y(find(y>maxx))=maxx;
49	y(find(y<minx))=minx;
50
51	end
52
53	y=y(1:M,1:N);