4th Order Chebyshev VST C++

[SharkeyO]

Hi friends,

I tried to implement a 4th order chebyshev-IIR-filter (IIR Form I), but all i get from it is digital white noise (or artefacts). Can anyone help me out?

Please have a look at the code:

Code: Select all

for (int32 i = 0; i < numChannels; i++)
					{
						int32 sampleFrames = data.numSamples;
						float* ptrIn  = in[i];
						float* ptrOut = out[i];
						float tmp;

						// el coefficientes
						float a0 = 0.0042412407f;
						float a1 = 0.016964963f;
						float a2 = 0.025447443f;
						float a3 = 0.016964963f;
						float a4 = 0.0042412407f;
						
						float b0 = 1.f;
						float b1 = -2.7280328f;
						float b2 = 3.2549775f;
						float b3 = -1.9259477f;
						float b4 = 0.47514287f;

						while (--sampleFrames >= 0)
						{
							// apply gain
							tmp=(*ptrIn++);
							(*ptrOut++) = (b0*tmp+b1*hist1+b2*hist2+b3*hist3+b4*hist4)/(a0*tmp+a1*hist1+a2*hist2+a3*hist3+a4*hist4);

							hist3=hist2;
							hist2=hist1;
							hist1=tmp;

						}
					}

Do i misunderstand something.

Thanks for your help.

André

[asomers]

Perhaps I'm misunderstanding your use of direct-form 1, but it looks like there are a number of flaws here.

First of all, what you have posted is not an IIR filter. You are only storing input samples. Also it seems that you are only updating 3 of your 4 history samples.

Next, it seems there is a misinterpretation about moving from transfer function to difference equation. The IIR (denominator) part of the transfer function is supposed to be subtracted from the FIR part (numerator), not divided. see https://ccrma.stanford.edu/~jos/fp/Z_Tr ... tions.html (focus on the relationship between the top eqn y(n) and the bottom eqn H(z).

I don't know why you have an a0 coefficient at all. A fourth order filter will have 5 feedforward coefficients (b0 - b4) and 4 feedback coefficients (a1 - b4). Only past output samples are involved in the feedback.

A review of this might be useful. https://ccrma.stanford.edu/~jos/fp/Direct_Form_I.html

[antto]

also, if NumChannels is greater than 1, you'll get another problem..

[duncanparsons]

First thing I noted was the channels issue! André, you'll need to keep separate histories for each channel, and do it across process calls. An easy way to do it could be to have arrays of histories (say InHist[chans][order] & OutHist[chans][order]) as part of the class, and reset everything to zero in resume. You might want to do denormal checking as well.

HTH

[SharkeyO]

Hey,

Thank you guys. I just tried to do like you told me and this is what i come up with. But i now don't get anything out of the speakers. The following code shows what i have done:

Code: Select all

for (int32 i = 0; i < numChannels; i++)
					 {
						
						int32 sampleFrames = data.numSamples;
						float* ptrIn1  = in[0];
						float* ptrOut1 = out[0];
						float* ptrIn2  = in[1];
						float* ptrOut2 = out[1];
                        float tmpl;
						float tmpr;

						// el coefficientes
						float a0 = 0.0042412407f;
						float a1 = 0.016964963f;
						float a2 = 0.025447443f;
						float a3 = 0.016964963f;
						float a4 = 0.0042412407f;
                  
						float b0 = 1.f;
						float b1 = -2.7280328f;
						float b2 = 3.2549775f;
						float b3 = -1.9259477f;
						float b4 = 0.47514287f;

						while (--sampleFrames >= 0)
						{
						  // apply gain
				          if(i==0)
						   {
						    
							tmpl = (b0*(*ptrIn1++)+b1*histl1+b2*histl2+b3*histl3+b4*histl4)-(a1*histl1+a2*histl2+a3*histl3+a4*histl4);
                            
							(*ptrOut1++)=tmpl;

						    histl4=histl3;
							histl3=histl2;
						    histl2=histl1;
						    histl1=tmpl;
						   } 

						  if(i==1)
						   {
						    							tmpr = (b0*(*ptrIn2++)+b1*histr1+b2*histr2+b3*histr3+b4*histr4)-(a1*histr1+a2*histr2+a3*histr3+a4*histr4);
                            
							(*ptrOut2++)=tmpr;

						    histr4=histr3;
							histr3=histr2;
						    histr2=histr1;
						    histr1=tmpr;
						   } 
					    } 

Where is the problem?!

Shouldn't this be the right code per channel:

Code: Select all

//tmpl = (a0*(*ptrIn1++)+a1*histl1+a2*histl2+a3*histl3+a4*histl4)/(b0+b1*histl1+b2*histl2+b3*histl3+b4*histl4);

Thank you for your help.

[DaveHoskins]

The history of the output and input needs to be separate:

Code: Select all

output = a0 * input + a1 * x1 + a2 * x2 + b1 * y1 + b2 * y2;
x2 = x1; x1 = input;
y2 = y1; y1 = output;

More info here:
http://www.dspguide.com/ch19/1.htm

Those 'if' statements look cumbersome also, instead of the 'if' statements why not have something like:
float hist1[2];
and index it by 'i'
You can also have
float* ptrIn = in;

Cheers,
Dave.

[antto]

SharkeyO: you didn't get anything from the speakers, because of this:

Code: Select all

for (int i = 0; i < NumChannels; i++)
{
	// the first time this loop is processed, "i" is 0, and it'll be incremented after the loop finishes (or before the next repetition, doesn't matter)
	...
	if (i == 1)
	{
		// this is where you've placed your code, but
		// "i" is still 0 ..
	}
}

[SharkeyO]

Thanks for your help. i can hear something now, that sounds like it's filtered. But i don't think it sound's correct. As i have been measuring the signal i t looks like i do have a problem with the coeffs, because the filter should be a LP filtering cutoff at 600Hz - I measured a cutoff at 4000Hz, which is very strange. I calculated me a chebychef filter of 6th order at http://www.dsptutor.freeuk.com/IIRFilte ... es102.html

the coeefs i got out look like

a[0] = 8.684399E-6
a[1] = 5.2106396E-5
a[2] = 1.30266E-4
a[3] = 1.73688E-4
a[4] = 1.30266E-4
a[5] = 5.2106396E-5
a[6] = 8.684399E-6

b[0] = 1.0
b[1] = -5.248061
b[2] = 11.779576
b[3] = -14.449645
b[4] = 10.205772
b[5] = -3.9329484
b[6] = 0.6459298

i now do post the left channel of my processing routine

Code: Select all

// el coefficientes
float a0 = 0.000008684399f;
float a1 = 0.000052106396f;
float a2 = 0.000130266f;
float a3 = 0.000173688f;
float a4 = 0.000130266f;
float a5 = 0.000052106396f;
float a6 = 0.000008684399f;
                 
float b0 = 1.f;
float b1 = -5.248061f;
float b2 = 11.779576f;
float b3 = -14.449645f;
float b4 = 10.205772f;
float b5 = -3.9329484f;
float b6 = 0.6459298f;

Code: Select all

inputl=(*ptrIn1++);
outputl = a0 * inputl + a1 * x1l + a2 * x2l + a3 * x3l + a4 * x4l + a5 * x5l + a6 * x6l 
- b1 * y1l - b2 * y2l - b3 * y3l - b4 * y4l - b5 * y5l - b6 * y6l;
					  
					  x6l = x5l;
					  x5l = x4l;
					  x4l = x3l;
					  x3l = x2l;
					  x2l = x1l; 
					  x1l = inputl;
                      
					  y6l = y5l; 
					  y5l = y4l; 
					  y4l = y3l; 
					  y3l = y2l; 
					  y2l = y1l; 
					  y1l = outputl;

					  (*ptrOut1++)=outputl;

What have i done wrong again?
Thanks for your patience.

Cheers
André

[antto]

SharkeyO, nothing wrong with the filter
what's wrong is that the coefficients don't match with the sampling rate

the filter design aplet you linked, uses a fixed Sampling Rate of 8KHz
so, when you entered 600Hz into that, it made you a filter with cutoff frequency of (600/8000)=0.075
yes, 0.075 is the actual frequency (that's not in Hz)
so if you run that filter at 44100Hz sampling rate, the cutoff will correspond to (0.075*44100)=3307.5Hz (yup, that's kinda close to what you measured)

you can design a filter with this aplet, which will work at 44100 and filter at 600Hz but you gotta do some math and enter a different cutoff frequency there
that would be ((600/44100) * 8000)=108.843Hz

but then again, the result will be correct only when you run the filter at 44100Hz sampling rate

the solution is to find the algorithms for calculating the filter coefficients, and you'll be able to run the filter at any freq/samplerate

[SharkeyO]

Yeah, i just didn't see it. Thank you for your help.

[asomers]

That's because that website uses a 8000Hz sampling rate! You could oversample that filter (I think? someone correct me). Or, just use a filter designer that lets you specify the sampling rate. I just tried this one: http://www-users.cs.york.ac.uk/~fisher/ ... /mkfscript

Note the different format of difference equation (signs flipped on feedback terms), and the gain term that scales the feedforward terms.

EDIT Antto beat me to it