Direct Convolution in MATLAB

Discussion in 'Programmer's Corner' started by battlehands, Aug 18, 2011.

  1. battlehands

    Thread Starter New Member

    Jul 29, 2011
    13
    0
    I need to perform convolution between two a filter and 5 signals in MATLAB without using the conv or filter commands.

    My filter is a sampled cosine function multiplied by some constant.

    h(t) = K cos(t) , [-pi/2, pi/2]

    My signals are in columns of a matrix.

    Column 1 = subsignal1
    Column 2 = subsignal2

    and so forth.

    I would like to write a loop to perform each of these convolutions individually, or all at once if possible.

    Something like:

    /code

    x = the 1st column in my matrix;
    N = length(x);
    h1=fliplr(h);

    yc = x;
    for t = 3:N-1
    b = x(t-2:t+1)
    s = b*h;
    yc(t) = s;
    end

    /code
     
  2. steveb

    Senior Member

    Jul 3, 2008
    2,433
    469
    Write out the discrete time version of the convolution sum using the sigma notation. Then consider the time limits on your signals to change the sum to be a finite sum. Then implement the summation as a loop directly as written.

    The nice thing about the "sigma" notation is that it directly translates into a coded for-loop.
     
  3. battlehands

    Thread Starter New Member

    Jul 29, 2011
    13
    0
    I think the limits are what is giving me trouble. I have 882000 samples in my signal, so each subsignal has 176400 samples. Therefore, I need my filter to sweep from at least 1:176400.
     
  4. steveb

    Senior Member

    Jul 3, 2008
    2,433
    469
    OK, and you need to consider the limits of the impulse response function too.

    Try writing out the formula for the convolution sum step by step, starting with the general formula and then substituting in your actual x[n] and h[n], then it will be come clearer where the limits should be. For every n= 1,2,3 ... up to your maximum value, there is a relevant upper and lower limit for your sum. It's very hard to see it until you write it out formally, and think very carefully about where x[n] times h[n] becomes zero, hence allowing you to set the limits. In other words, there is no point (because it's impossible and doesn't change the answer) in summing zeros an infinite number of times.
     
  5. battlehands

    Thread Starter New Member

    Jul 29, 2011
    13
    0
    The definition Im following is:

    y(n) = Ʃ_k=-infty^infty x(k)h(n-k)

    by plugging in my actual x[n] and h[n] do you mean take the end points and start doing point-wise calculations for values of y?

    Do you have an example of what you are saying?
     
  6. steveb

    Senior Member

    Jul 3, 2008
    2,433
    469
    Well, your example is good to use.

    So, your first step is good. Now think about the limits on h and x. You already posted the limits on x above, so you need to also look at the limits on h.

    Before we do this, we need to think about one issue with Matlab. Matlab does not allow the specification of a vector at an index of 0. You must start at 1 which is not consistent with usual discrete time notation because you would like to have h[n] defined at h[0], and also the input signal at x[0] is relevant.

    So, I'm going to assume your signals start at x[0], but assume that the first value is zero (x[0]=0) and fortunately your impulse response function just happens to be zero, which is not always the case.

    So, your h[n] has 27 values for n=0, 1, 2 ... 26, with h[0]=h[26]=0. This means that n-k in your convolution sum must be in the range of 1 to 25, and you don't need to consider summing outside that range.

    Try to proceed from here.
     
  7. battlehands

    Thread Starter New Member

    Jul 29, 2011
    13
    0
    Something like this?

    /code

    [nxr, nxc] = size(sig);
    [nmr, nmc] = size(h);

    nyr = nxr - nmr + 1;
    nyc = nxc - nmc + 1;

    yc = zeros(nyr, nyc);

    for z = 1:nmr
    for a = 1:nmc
    subsig:),z) = sig(nmr-z+1:nxr-z+1, nmc-a+1:nxc-a+1);
    yc = y + subsig:),z) * h(z, a);
    end
    end

    /code
     
  8. steveb

    Senior Member

    Jul 3, 2008
    2,433
    469
    Before trying to code it, you should express it mathematically. Only when you have clearly specified what you are trying to do should you convert to code. It's still not clear to me that you have the limits correct, and I don't feel ready to try to interpret uncommented code until I see the mathematical expression clearly stated. Actually, once the clear mathematical statement is made, comments are not very critical if you use the same variable names.
     
  9. battlehands

    Thread Starter New Member

    Jul 29, 2011
    13
    0
    How can I do it by hand when my signal and filter both contain so many points?

    Im getting lost with all of these variables. h, j, k ,n
     
  10. steveb

    Senior Member

    Jul 3, 2008
    2,433
    469
    What I mean is to write it out mathematically using the sigma notation, just as you started to do. The key is to identify the correct limits as a function of n. The benefit of sigma notation is that it converts directly to a loop in code.

    The problem with your first formula

     y[n]=\sum_{k=-\infty}^{+\infty} x[k] \; h[n-k]

    is that you can't make a for loop that sums an infinite number of terms. However, there is no need to sum over infinity because the product x[k] h[n-k] is zero accept for particular values of k for any given n.

    There are two limits you are concerned with. The range of n and the range of k are both important. For example, it's clear that the output y[n] is zero for n less than zero. So, you don't care about calculating for negative n. Also, because your system impulse response function is finite with a length of 27 samples, there is no need to calculate for n > 176427.

    The range of k for any particular n is also important. So the limits on the sum could depend on the value of n. So the formula you should write out should look something like the following.


     y[n]=\sum_{k=L[n]}^{U[n]} x[k] \; h[n-k]

    where L[n] and U[n] are the functions for the lower and upper limits respectively. Once you know L and U, you can directly write a for loop in code that does the job. I know it's really tricky to think about the first time you do this, but this is the process that leads to the answer. You really have to think carefully on these types of problems.
     
Loading...