function HH = plot_2D_filt(h,M);
% Computes and plots the two-dimensional response of the filter
% with coefficients given in the matrix h
% M is the number of points in the frequency discretization

% defaults to 64 samples
if nargin < 2 ; 
    M = 64;
end

% Number of taps of the FIR filter
N = max(size(h));  % h should be square!

% Computes the frequency response
% (simple and straightforward way, use dct2 instead for efficiency)
H = zeros(M);
for p=1:M
    w1 = pi*((p-1)/M);
    cosw1k = cos(w1*[0:(N-1)]);
    for q=1:M
        w2=pi*((q-1)/M);
        cosw2l = cos(w2*[0:(N-1)])';
        H(p,q) = cosw1k*h*cosw2l;
    end
end

% Reflect boundaries
HH = [H(end:-1:2,end:-1:2) H(end:-1:2,1:end) ; H(1:end,end:-1:2) H];

% Axes and plot
t=linspace(-pi,pi,2*M-1);
mesh(t,t,HH);
xlabel('\omega_1');
ylabel('\omega_2');
zlabel('H(\omega_1,\omega_2)');