function [xFine, fFine] = m20120528_03_FourierInterpolation(f, xMin, xMax, NFine, boundary) % Interpolate a periodic function from one uniform grid to a finer grid, % using Fourier transforms to achieve spectral accuracy. % Possible values for boundary: % 1 = first element is xMin, last element is xMax-dx. % 2 = first element is xMin+dx, last element is xMax. % 3 = first element is xMin, last element is xMax, so f(1)=f(end) if NFine==numel(f) fFine=f; switch boundary case 1 xFine = linspace(xMin,xMax,NFine+1); xFine(end)=[]; case 2 xFine = linspace(xMin,xMax,NFine+1); xFine(1)=[]; case 3 xFine = linspace(xMin,xMax,NFine); otherwise error('Invalid value of boundary') end return end if NFine < numel(f) error('Not set up yet for reducing resolution') end column=1; row=2; if iscolumn(f) orientation=column; f=f'; else orientation=row; end % Internally, everything will be done with rows. We'll transpose at the % end if needed; switch boundary case 1 xFine = linspace(xMin,xMax,NFine+1); xFine(end)=[]; case 2 xFine = linspace(xMin,xMax,NFine+1); xFine(1)=[]; f = [f(end), f(1:(end-1))]; case 3 if abs(f(1)- f(end))> 1e-13 error('For boundary=3, f(1) must equal f(end)') end xFine = linspace(xMin,xMax,NFine); f(end)=[]; NFine=NFine-1; otherwise error('Invalid value of boundary') end NCoarse = numel(f); ff=fft(f); if mod(NCoarse,2)==0 halfN = ceil(NCoarse/2); ffFine = [ff(1:(halfN)), 0.5*ff(halfN+1), zeros(1,NFine-NCoarse-1), 0.5*ff(halfN+1), ff((halfN+2):end)] * NFine/NCoarse; else halfN = ceil(NCoarse/2); ffFine = [ff(1:halfN), zeros(1,NFine-NCoarse), ff((halfN+1):end)] * NFine/NCoarse; end fFine = ifft(ffFine); switch boundary case 2 fFine = [fFine(2:end), fFine(1)]; case 3 fFine = [fFine, fFine(1)]; end if orientation==column xFine=xFine'; fFine=fFine'; end end