function H = BP_fn(bw,gamma,psi,lambda,theta)
sigma = lambda/pi*sqrt(log(2)/2)*(2^bw+1)/(2^bw-1);
sigma_x = sigma;
sigma_y = sigma/gamma;
sz = fix(8*max(sigma_y,sigma_x));
if mod(sz,2)==0, sz = sz+1;end
[x y] = meshgrid(-fix(sz/2):fix(sz/2),fix(sz/2):-1:fix(-sz/2));
% Rotation
x_theta = x*cos(theta)+y*sin(theta);
y_theta = -x*sin(theta)+y*cos(theta);
H = exp(-0.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2))
.*cos(2*pi/lambda*x_theta+psi);