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	dimitri	1.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);