function xx = frot(im,theta)
% Fourier domain rotations
[N M] = size(im);
[nn mm] = ndgrid([-N/2:N/2-1],[-M/2:M/2-1]);
e1 = exp(-i*2*pi*tan(theta/2)*nn.*mm./M);
e2 = exp(i*2*pi*sin(theta)*nn.*mm./N);
r1 = zeros(size(im));
r2 = zeros(size(im));
xx = zeros(size(im));
for lp=1:N
  r1(lp,:) = ifft(fftshift(fftshift(fft(im(lp,:))).*e1(lp,:)));
end
for lp=1:M
  r2(:,lp) = ifft(fftshift(fftshift(fft(r1(:,lp))).*e2(:,lp)));
end
for lp=1:N
  xx(lp,:) = ifft(fftshift(fftshift(fft(r2(lp,:))).*e1(lp,:)));
end