myconvper(z,h)

%MYCONVPER is a periodical fft accelerated correspondence to
%   ZC = CONV(Z,H)
%   
%   Here follows an examplification of usage:
%   
%   clear all;
%   
%   %  Create a "measurement" domain
%   Lx = 1e-3;
%   Nx = 2^15;
%   x  = Lx*(0:Nx-1)/Nx;
%   dx = x(2)-x(1);
%   
%   % "Sample" the (periodic) signal
%   zo = cos(2*pi*2*x/Lx)+...
%        0.5*cos(2*pi*5*(x-0.14*Lx)/Lx)+...
%        0.25*cos(2*pi*20*(x-0.3*Lx)/Lx)+...
%        0.25*cos(2*pi*30*(x-0.7*Lx)/Lx);
% 
%   z  = [zo(end-Nx/2+1:end),zo,zo(1:Nx/2)]; % Periodic ext. removing edge-effects
%   Nz = length(z);
%   
%   % Filter weighting function
%   lcx = Lx/5; % Cut-off acc to ISO 11562 Gaussian Profile Filter Lx/5
%   xf  = [x x(end)+dx]-Lx/2; % Symmetric and ensure that Nx+Nh-1 = 2^k
%   variancex = 1/(pi*sqrt(2/log(2)))*lcx; 
%   h = dx/(sqrt(2*pi)*variancex)*...
%        exp(-(((xf)/variancex).^2)/2);
%   Nh  = length(h);
% 
%   % Total length of the padded arrays needed for the fft operations
%   N  = Nz+Nh-1;
%   
%   % Apply the filter in the spatial domain using conv
%   tic;
%   zc = conv(z,h);
%   zc = zc(Nx+1:N-Nx); 
%   toc;
% 
%   % Apply the filter in the frequency domain
%   tic;
%   Z  = fft([z zeros(1,Nh-1)]);
%   H  = fft([h zeros(1,Nz-1)]);
%   zf = real(ifft(Z.*H));
%   zf = zf(Nx+1:N-Nx); 
%   toc;
% 
%   % Doing the same again but using the functional representation
%   tic;
%   zffunc = myconvper(zo,h);
%   toc;
% 
%   % Plotting the results
%   plot(...
%        x, zo    , 'b-' ,...Original data
%        x, zc    , 'k.' ,...Filtered - spatial domain
%        x, zf    , 'r--',...Filtered - freq. domain
%        x, zffunc, 'w:'  ...Filtered - freq. domain, func. repr.
%        );

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

[Ny,Nx] = size(z);
z  = [z(:,end-Nx/2+1:end),z,z(:,1:Nx/2)]; % Periodic extension to remove
                                          % the edge effects
zc = myconv(z,h);

Nz = size(z,2);
Nh = size(h,2);
N  = Nz+Nh-1;

% zc = real(ifft(fft([z zeros(Ny,Nh-1)]').*fft([h zeros(Ny,Nz-1)]')))';

zc = zc(:,Nx+1:N-Nx);