3
$\begingroup$

I'm trying to calculate convolution of two given vectors in MATLAB without using loop, and of course without the function conv itself, but I can't remove the last loop I've used in the code below. Any hint would be appreciated.

function [ convolve ] = myconv( x,y)
%This function calculate the convolution of the two inputs
lx=length(x);
ly=length(y);
lt=lx+ly;
x=[x zeros(1,ly)];
y=[y zeros(1,lx)];
convolve=zeros;
for n=2:lt
    convolve(n-1)=sum(x(1:n-1).*fliplr(y(1:n-1)));
end
$\endgroup$
  • $\begingroup$ general programming questions are offtopic here (might be better on StackOverflow, if improved), and asking for code written to a specification is also off-topic here; I don't really know what to salvage from this question, sorry. $\endgroup$ – Marcus Müller Mar 13 '17 at 11:15
  • $\begingroup$ Sorry,I'm new in here ,I will do better the next time. $\endgroup$ – MAh2014 Mar 13 '17 at 12:30
5
$\begingroup$

what about this:

function [convolve] = myconv(x, y)
    L = length(x) + length(y) - 1
    convolve = ifft(fft(x, L) .* fft(y, L));
end

Using the convolution theorem.

Regarding deconv function, check this out:

x = randn(10,1);
y = randn(10,1);

L = length(x) + length(y) - 1;

z = ifft(fft(x, L) .* fft(y, L));
z2 = conv(x, y);
sum(abs(z-z2))

x2 = ifft(fft(z, L) ./ fft(y, L));
x3 = deconv(z, y);
x0 = [x; zeros(L-length(x),1)];
sum(abs(x0-x2))
sum(abs(x3-x))

Program output:

ans =

   1.1616e-14

ans =

   5.4956e-15


ans =

   2.3924e-12

The first sum calculates the deviation between conv and my proposal. The second sum shows difference between the Fourier-Domain decovonvoled sequence and original sequence (note that the output of the FFT has L samples, where the last L-len(x) samples are 0). The third sum compares the output of deconv and the original sequence. You see that they all match up to floating point errors.

$\endgroup$
  • $\begingroup$ Yeah,I think that is the best,just one question,Is it also true for deconvolution ,I mean is$ deconv=ifft(fft(h)./fft(x))$,if $ h=conv(x,y)$.I tried but got different answers?Where am I doing wrong? $\endgroup$ – MAh2014 Mar 13 '17 at 12:23
  • $\begingroup$ it's better to explicitly constraint the output of the ifft to be real in case for for real inputs, as many implementations will treat it as complex valued due to very small roundoff errors introduced which practically being nonzero formally yields imaginary components... $\endgroup$ – Fat32 Mar 13 '17 at 12:28
  • $\begingroup$ @Fat32 that's correct. But I think, the OP should be able to do this on his own, in case it's a real signal (it's not stated its real). $\endgroup$ – Maximilian Matthé Mar 13 '17 at 12:36
  • $\begingroup$ @MAh2014 In principle this should work also. But, take care of the fft lengths, so it should be y=ifft(fft(h, L) ./ fft(x, L)); $\endgroup$ – Maximilian Matthé Mar 13 '17 at 12:41
  • $\begingroup$ @Fat32, Yes that's right I'm dealing with real signals,but how could I get the exact answer? $\endgroup$ – MAh2014 Mar 13 '17 at 12:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.