function [out] = hp_filter(in)

% Get image size
[Y,X] = size(in);

% Build filter
[Xm,Ym] = meshgrid(1:X,1:Y);
Xm = (Xm - X/2) / X;
Ym = (Ym - Y/2) / Y;
filt = cos(pi * Xm) .* cos(pi * Ym);
filt = (1 - filt) .* (2 - filt);

% Apply filter
out = in .* filt;