To my understanding simply extend the input arrays such that you have a couple of points from the first few samples added onto the end adding 2*pi to those angles and a few points from the end of your arrays at the beginning subtracting 2*pi from those:
theta = linspace(0,2*pi-0.001,31);
f = 12*cos(theta) + 5*sin(3*theta).^2;
Theta = [theta(end-3:end)-2*pi,theta,theta(1:3)+2*pi];
F1 = [f(end-3:end),f,f(1:3)];
THETA = linspace(0,2*pi,301);
F2 = 12*cos(THETA) + 5*sin(3*THETA).^2;
Fi = interp1(Theta,F1,THETA,'pchip');
That way you force the interpolation to be periodic, you only need to braket your theta and f array with three (possibly 2 but I'm too lazy to check) points to get the cubic and spline interpolation-methods to give you a periodic interpolation.
HTH
