Sign in to follow this  
TheMadScientist

Convolution not working

Recommended Posts

I thought I wrote this simple convolution function properly but it doesnt seem to work. When I convolve a sine wave with a step function or any signal I get alot of mess. void Convolve(float *InputSignal,long N,float *ImpulseResponse,long M,float *OutputSignal) { long i = 0; long j; ZeroFill(OutputSignal,N); while ( i < N+M) { j=0; while ( j < M) { if((i-j) >= 0) // I ignore the first negative indexes. { OutputSignal[i] = OutputSignal[i] + ImpulseResponse[j]*InputSignal[i-j]; } j++; } i++; } } The Reason I would like to use a convolution function is to impliment FIR filtering later on. If anyone sees what I am doing wrong I would really appreciate seeing the proper code. Thanks

Share this post


Link to post
Share on other sites
Put your code within [code][/code] tags to format it properly for this board.

Anyway, I don't see anything wrong with the algorithm. Convolving a sine wave with a step function *should* result in a lot of mess (depending on the precision of your numeric integration). If you were expecting another sinusoid out, you should have convolved the sine with an *impulse*.

That said, your code could stand some cleanup. First off, why not just use for loops? Also, for accumulating values, you can use the immediate operator '+=' which is a little more expressive IMHO. For "ignoring negative indices", you could just start 'j' off with the current 'i' value in the inner loop. Oh, and declare numeric values as 'int' by default unless you know you need something else - it's the shortest type identifier for a reason, to encourage you to use it as default.

Assuming plain C (you poor thing), you can still clean things up this far:


void Convolve (float *InputSignal, long N, float *ImpulseResponse, long M, float *OutputSignal) {
ZeroFill(OutputSignal,N);
for (int i = 0; i < N + M; ++i) {
for (int j = i; j < M; ++j) {
OutputSignal[i] += ImpulseResponse[j]*InputSignal[i-j];
}
}
}


In C++ I would make use of std::vector:


typedef std::vector<float> signal;
signal convolution(const signal& input, const signal& impulseResponse) {
int inputLength = input.size();
int impulseLength = impulseResponse.size();
int resultLength = inputLength + impulseLength;
signal result(resultLength);
for (int i = 0; i < resultLength; ++i) {
for (int j = i; j < impulseLength; ++j) {
result[i] += impulseResponse[j]*input[i-j];
}
}
return result;
}


That might not look much cleaner (it got rid of some parameters at least), but proper use of std::vector is likely to solve a lot of memory management issues in the context of the calling code. It also lets you safely return the result value instead of writing out data at some raw pointer address.

Share this post


Link to post
Share on other sites
You have no idea how much this helps me out.Zahlman_Rating++; Thanks alot. You spoke of negative indexes which brings up another question. In a recursive filter you have negative indexes as well like y[n-5] or nonexistend ones like the one just mentioned. What do you do at that point? Do you substitute y[n-5] with a zero? Excuse my lack of dsp knowledge but I havent had any formal classes and have been trying to learn everything from tutorials and the pdf version of the scientist and Engineer's guide to DSP.

Share this post


Link to post
Share on other sites
Yes, indices outside the range of the singal and/or impulse response are usually assumed to be zero. If you don't make this assumption and only calculate the parts with valid indices, you will miss M-1 samples at the beginning of the output signal as you simply can't calculate those samples.

I suggest you make a pre-processing step where you calculate the first M-1 samples using a special loop, and then convolve the rest as normal with the main loop. Then you don't have to extend the input signal with more values, nor do you have to insert, possibly costly, logics to determine how much of the impulse repsonse to use.

Share this post


Link to post
Share on other sites
Thanks All of this has solved my problems with convolution but I am still puzzled about negative indexes in recursive filters both the [n- and [y- indexes. If one of my coefficients calls for y[n-5] and n=0 I would assume this would be zero but wouldn't that mean the same end problems in convolution would arise in recursion? Well here is my code for my recursive filter and I believe it doesnt work because I am mishandling the indexes. The function runs
without error so if you try it in VC++ it will work just not with the desired
output. I use longs because some data I have processed is too large for ints.



void NotchBandPass(float *InputSignal,float FC,float BW,float *OutputSignal,long N)
{

long i = 0;

memset(OutputSignal,0,N);

float R = (float)(1 - (3*BW));

float K =(float) (1-2*R*(cos(2*PI*FC))+R*R);
K =(float) (K/(2-2*cos(2*PI*FC)));

float A0,A1,A2,B1,B2;
A0=(float)(1-K);
A1=(float)(2*(K-R)*cos(2*PI*FC));
A2=(float)(R*R - K);
B1=(float)(2*R*cos(2*PI*FC));
B2=(float)((-1)*(R*R));

i= 0;

while( i < N)
{

OutputSignal[i]=(float)( A0*InputSignal[i] + A1*InputSignal[i-1] + A2*InputSignal[i-2]
+ B1*OutputSignal[i] + B2*OutputSignal[i-1]);

i++;
}


}




[Edited by - TheMadScientist on June 12, 2006 12:57:42 PM]

Share this post


Link to post
Share on other sites
Quote:
Original post by TheMadScientist
If one of my coefficients calls for y[n-5] and n=0 I would assume this would be zero

I think I've found your problem: your array indexes are negative, so values they are pointing to are meaningless and have random values.
When i=0 you are using theese elements:
- InputSignal[0] ok
- InputSignal[-1] <<-- bad
- InputSignal[-2] <<-- bad
- OutputSignal[0] ok
- OutputSignal[-1] <<-- bad
Remember arrays first index is 0 and there is nothing before.

You should do something like this.
Suppose you need nNeg element before zero index.


float *InputSignal=new float[N+nNeg+1];
float *OutputSignal=new float[N+nNeg+1];

InputSignal+=nNeg+1;
OutputSignal+=nNeg+1;
// now InputSignal[-nNeg] exists because actually it is InputSignal[0]

...

NotchBandPass(InputSignal,FC,BW,OutputSignal,N);

...

// dispose mem (delete must use original array so I need to change the two pointers)
InputSignal-=nNeg+1;
OutputSignal-=nNeg+1;
delete [] InputSignal;
delete [] OutputSignal;




Now your NotchBandPass should work... if there are not other mistakes [grin]

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this