myconv2(z,h)

%MYCONV2 is a fft accelerated correspondence to
%   ZC = CONV2(Z,H)
%
%   Z is a Ny x Nx matrix and H is Ny x Nh matrix
%   ZC is a 2*Ny x (Nx+Nh-1) matix
%   
%   See also CONV2, MYCONVPER2, CONV, MYCONV, MYCONVPER
%
%   Here follows an examplification of usage:
%   
%   clear all;

%   % Create a "measurement" domain
%   Lx = 1e-3;
%   Ly = 1e-3;
%   Nx = 2^6;
%   Ny = 2^5;
%   x  = Lx*(0:Nx-1)/Nx;
%   dx = x(2)-x(1);
%   y  = Ly*(0:Ny-1)/Ny;
%   dy = y(2)-y(1);
% 
%   [X,Y] = meshgrid(x,y);
% 
%   % "Sample" the (periodic) signal
%   z  = cos(2*pi*2*X/Lx);
% 
%   % Filter weighting function
%   lcx = Lx/2;  % Cut-off acc to ISO 11562 Gaussian Profile Filter Lx/5
%   lcy = Ly/2;  
%   Xf = X-Lx/2; % Symmetric and ensure that Nx+Nh-1 = 2^k
%   Yf = Y-Ly/2;
%   variancex = 1/(pi*sqrt(2/log(2)))*lcx; 
%   variancey = 1/(pi*sqrt(2/log(2)))*lcy; 
%   hx = dx/(sqrt(2*pi)*variancex)*...
%     exp(-((Xf/variancex).^2)/2);
%   hy = dy/(sqrt(2*pi)*variancey)*...
%     exp(-((Yf/variancey).^2)/2);
%   Nh  = size(hx,2);
% 
%   % Total length of the padded arrays needed for the fft operations
%   N  = Nx+Nh-1;
% 
%   % Apply the filter in the spatial domain using conv2
%   tic;
%   zc = conv2(z,hx.*hy);
%   toc;
% 
%   % Apply the filter in the frequency domain
%   tic;
%   Z  = fft2([z     , zeros(Ny,Nh-1); zeros(Ny-1,N)].');
%   H  = fft2([hx.*hy, zeros(Ny,Nx-1); zeros(Ny-1,N)].');
%   zf = real(ifft2(Z.*H))';
%   toc;
% 
%   % Doing the same again but using the functional representation
%   tic;
%   zffunc = myconv2(z,hx.*hy);
%   toc;
% 
%   % Plotting the results
%   plot(...
%     1:N, zc(1,:)    , 'k.' ,...Filtered - spatial domain
%     1:N, zf(1,:)    , 'r--',...Filtered - freq. domain
%     1:N, zffunc(1,:), 'w:'  ...Filtered - freq. domain, func. repr.
%     );

%   Author(s): A. Almqvist
%   Copyright 2009- Andreas Almqvist
%   $Revision: 1$  $Date: 2009/11$
%   Homemade function
function zc = myconv2(z,h)

[Nzy,Nzx]   = size(z);
[Nhy,Nhx] = size(h);

% Total length of the padded arrays needed for the fft operations
N = Nzx+Nhx-1;    %M = Nzy+Nhy-1;

% Compute the convolution utilizing the fft
zc = real(...
       ifft2(...
         fft2([z, zeros(Nzy,Nhx-1); zeros(Nhy-1,N)].').*...
         fft2([h, zeros(Nhy,Nzx-1); zeros(Nzy-1,N)].')...
         )...
       )';