0

I want to take the discrete convolution of two 1-D vectors. The vectors correspond to intensity data as a function of frequency. My goal is to take the convolution of one intensity vector B with itself and then take the convolution of the result with the original vector B, and so on, each time taking the convolution of the result with the original vector B. I want the final result to have the same length as the original vector B.

I am starting with a code in IDL that I am trying to modify to MATLAB. The relevant part of the code reads:

for i=1,10 do begin
if i lt 2 then A=B else A=Rt
n=i+1
Rt=convol(A,B,s,center=0,/edge_zero)
...
endfor

which I have rewritten in MATLAB

for i = 2:11
    if i < 3
        A = B; % indices start at 0, not 1
    else
        A = Rt;
    end
    n = i + 1;

    % Scale by 1/s
    Rt = (1/s).*conv(A,B);
    ...
end

But I am not sure how to incorporate the zero-padding that uses the edge_zero option. In IDL, the convolution calculates the values of elements at the edge of the vector as if the vector were padded with zeroes. The optional third option for the convolution function in MATLAB includes the option 'same', which returns the central part of the convolution of the same size as u for conv(u,v), but that doesn't seem to be the correct way about this problem. How do I do an analogous zero-padding in MATLAB?

idladle
  • 43
  • 4
  • 1
    MATLAB always uses zero padding for convolution. Why do you say that 'same' is not the correct way? Are you getting different results? – Cris Luengo Mar 21 '19 at 03:38

1 Answers1

0

Here is a code I needed for my doctoral research that I know does the zero padding correctly. I hope it helps.

function conv_out=F_convolve_FFT(f,g,dw,flipFlag),
 if(nargin<4), flipFlag==0; end;

 % length of function f to be convolved, initialization of conv_out, and padding p
 N=length(f); conv_out=zeros(1,N); p=zeros(1,N);

 % uncomment if working with tensor f,g's of rank 3 or greater
 % f=F_reduce_rank(f); g=F_reduce_rank(g);

 % padding. also, this was commented out: EN=length(fp);
 fp = [f p]; gp = [g p];

 % if performing convolution of the form c(w) = int(f(wp)g(w+wp),wp) = int(f(w-wp)gbar(wp),wp) due to reverse-order domain on substitution
 if(flipFlag==1), gp=fliplr(gp); end;

 % perform the convolution. You do NOT need to multiply the invocation of "F_convolve_FFT(f,g,dw,flipFlag)" in your program by "dx", the finite-element.
 c1 = ifft(fft(fp).*fft(gp))*dw;

 % if performing "reverse" convolution, an additional circshift is necessary to offset the padding
 if(flipFlag==1), c1=circshift(c1',N)'; end;

 % offset the padding by dNm
 if(mod(N,2)==0), dNm=N/2; elseif(mod(N,2)==1), dNm=(N-1)/2; end;

 % padding. also, this was commented out: EN=length(fp);
 conv_out(:)=c1(dNm+1:dNm+N);

return
Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
ignoramus
  • 13
  • 3