Forums

NaN issue with NOTCH filter

Started by jtp_1960 December 10, 2013
Hello!

I have an issue with a Notch filter which works OK for certain Hz/bw/Q but
after changing Hz or bw/Q sometimes a NaN value is given for the sample for
center frequency. 

http://i39.tinypic.com/32zrgk8.jpg

Source code:


double pi = 3.14159265358979323846f;
double Fs = 96000.0;
double ln = 0.69314718055994530942f;
double halfsq2 = 0.707106781186547524401;

for(int i = 1; i < 11; i++){
	double a0;
	double A;
	double w0;
	double f0;
	double bw = 0.0f;
	double G = 0.0f;
	int type = bandData[currentSpeakerMode][currentChannel][i].type;
	int BwQ = bandData[currentSpeakerMode][currentChannel][i].BwQNull;
	int swType;

	if(type < 3){ swType = 1;} // peaking
	else if(type == 3 || type == 4){ swType = 2;} //lp
	else if(type == 5 || type == 6){ swType = 3;} //hp
	else if(type == 7){ swType = 4;}//bp
	else if(type == 8 || type == 10 || type == 11){ swType = 5;}//ls
	else if(type == 9 || type == 12 || type == 13){ swType = 6;}//hs
	else if(type == 14){ swType = 7;}//notch
	else{ swType = 8;}//all pass

	if(bandActivity[i].Checked){
		f0 =
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].freqHz);
		if(bwOrQSliderNeedsToBeActive(i)){
			bw =
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].bwdB);
		}
		else if(type == 10 || type == 12){ bw = 6.0f;}
		else if(type == 11 || type == 13){ bw = 12.0f;}

		if(bw == 0.0f){
			if(swType > 1 && swType < 5){ bw = halfsq2;}
			else if(swType == 5 || swType == 6){ bw = 0.9f;}
			else if(swType == 7){ bw = 15.0f;}
		}
		else{
			bw /= 12.0f;
		}

		if(gainSliderNeedsToBeActive(i)){
			G =
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].gaindB);
		}

		if(swType == 1 || swType == 5 || swType == 6){
			A = pow(10, (G / 40.0f));
		}
		else{
			A = pow(10, (G / 20.0f));
		}

		w0 = 2*pi*f0/Fs;
		double alpha;
		double beta = -1.0f;

		if(swType == 5 || swType == 6){ // ls/hs
			alpha = sin(w0)/2 * sqrt((A + 1/A) * (1/bw - 1) + 2);
			beta = 2 * sqrt(A) * alpha; 
		}
		else if(BwQ == 2){ // bw
			alpha = sin(w0) * sinh( ln * bw * w0/sin(w0) );
		}
		else{ // q
			alpha = sin(w0) / (2 * bw);
		}

		...
		switch(swType){
		...
		case 7: // notch
			a0 = 1 + alpha;
			coefficient[i].b0 = 1;
			coefficient[i].b1 = -2 * cos(w0);
			coefficient[i].b2 = 1;
			coefficient[i].a1 = -2 * cos(w0);
			coefficient[i].a2 = 1 - alpha;
			break;
		...
		}

		// norm
		coefficient[i].b0 = coefficient[i].b0/a0;
		coefficient[i].b1 = coefficient[i].b1/a0;
		coefficient[i].b2 = coefficient[i].b2/a0;
		coefficient[i].a0 = 1.0f;
		coefficient[i].a1 = coefficient[i].a1/a0;
		coefficient[i].a2 = coefficient[i].a2/a0;

..

// calculate magnitude
numerator = coefficient[s].b0*coefficient[s].b0 +
coefficient[s].b1*coefficient[s].b1 + coefficient[s].b2*coefficient[s].b2 +

		2.0L*(coefficient[s].b0*coefficient[s].b1 +
coefficient[s].b1*coefficient[s].b2)*cosl(w) +
2.0L*coefficient[s].b0*coefficient[s].b2*cosl(2.0L*w);
denominator = 1.0L + coefficient[s].a1*coefficient[s].a1 +
coefficient[s].a2*coefficient[s].a2 + 
		2.0L*(coefficient[s].a1 + coefficient[s].a1*coefficient[s].a2)*cosl(w) +
2.0L*coefficient[s].a2*cosl(2.0L*w);
magnitude = 20*log10(sqrtl(numerator / denominator)); 

..


Is there something wrong in source code or is there other factors
involved?	 

_____________________________		
Posted through www.DSPRelated.com
On 12/10/13 8:20 AM, jtp_1960 wrote:
> Hello! > > I have an issue with a Notch filter which works OK for certain Hz/bw/Q but > after changing Hz or bw/Q sometimes a NaN value is given for the sample for > center frequency. > > http://i39.tinypic.com/32zrgk8.jpg > > Source code: > > > double pi = 3.14159265358979323846f; > double Fs = 96000.0; > double ln = 0.69314718055994530942f; > double halfsq2 = 0.707106781186547524401; > > for(int i = 1; i< 11; i++){ > double a0; > double A; > double w0; > double f0; > double bw = 0.0f; > double G = 0.0f; > int type = bandData[currentSpeakerMode][currentChannel][i].type; > int BwQ = bandData[currentSpeakerMode][currentChannel][i].BwQNull; > int swType; > > if(type< 3){ swType = 1;} // peaking > else if(type == 3 || type == 4){ swType = 2;} //lp > else if(type == 5 || type == 6){ swType = 3;} //hp > else if(type == 7){ swType = 4;}//bp > else if(type == 8 || type == 10 || type == 11){ swType = 5;}//ls > else if(type == 9 || type == 12 || type == 13){ swType = 6;}//hs > else if(type == 14){ swType = 7;}//notch > else{ swType = 8;}//all pass > > if(bandActivity[i].Checked){ > f0 = > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].freqHz); > if(bwOrQSliderNeedsToBeActive(i)){ > bw = > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].bwdB); > } > else if(type == 10 || type == 12){ bw = 6.0f;} > else if(type == 11 || type == 13){ bw = 12.0f;} > > if(bw == 0.0f){ > if(swType> 1&& swType< 5){ bw = halfsq2;} > else if(swType == 5 || swType == 6){ bw = 0.9f;} > else if(swType == 7){ bw = 15.0f;} > } > else{ > bw /= 12.0f; > } > > if(gainSliderNeedsToBeActive(i)){ > G = > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].gaindB); > } > > if(swType == 1 || swType == 5 || swType == 6){ > A = pow(10, (G / 40.0f)); > } > else{ > A = pow(10, (G / 20.0f)); > } > > w0 = 2*pi*f0/Fs; > double alpha; > double beta = -1.0f; > > if(swType == 5 || swType == 6){ // ls/hs > alpha = sin(w0)/2 * sqrt((A + 1/A) * (1/bw - 1) + 2); > beta = 2 * sqrt(A) * alpha; > } > else if(BwQ == 2){ // bw > alpha = sin(w0) * sinh( ln * bw * w0/sin(w0) ); > } > else{ // q > alpha = sin(w0) / (2 * bw); > } > > ... > switch(swType){ > ... > case 7: // notch > a0 = 1 + alpha; > coefficient[i].b0 = 1; > coefficient[i].b1 = -2 * cos(w0); > coefficient[i].b2 = 1; > coefficient[i].a1 = -2 * cos(w0); > coefficient[i].a2 = 1 - alpha; > break; > ... > } > > // norm > coefficient[i].b0 = coefficient[i].b0/a0; > coefficient[i].b1 = coefficient[i].b1/a0; > coefficient[i].b2 = coefficient[i].b2/a0; > coefficient[i].a0 = 1.0f; > coefficient[i].a1 = coefficient[i].a1/a0; > coefficient[i].a2 = coefficient[i].a2/a0; > > .. > > // calculate magnitude > numerator = coefficient[s].b0*coefficient[s].b0 + > coefficient[s].b1*coefficient[s].b1 + coefficient[s].b2*coefficient[s].b2 + > > 2.0L*(coefficient[s].b0*coefficient[s].b1 + > coefficient[s].b1*coefficient[s].b2)*cosl(w) + > 2.0L*coefficient[s].b0*coefficient[s].b2*cosl(2.0L*w); > denominator = 1.0L + coefficient[s].a1*coefficient[s].a1 + > coefficient[s].a2*coefficient[s].a2 + > 2.0L*(coefficient[s].a1 + coefficient[s].a1*coefficient[s].a2)*cosl(w) + > 2.0L*coefficient[s].a2*cosl(2.0L*w); > magnitude = 20*log10(sqrtl(numerator / denominator)); > > .. > > > Is there something wrong in source code or is there other factors > involved? >
exactly what variable becomes a NaN? -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On Tuesday, December 10, 2013 8:20:19 AM UTC-5, jtp_1960 wrote:
> Hello! > > > > I have an issue with a Notch filter which works OK for certain Hz/bw/Q but > > after changing Hz or bw/Q sometimes a NaN value is given for the sample for > > center frequency. > > > > http://i39.tinypic.com/32zrgk8.jpg > > > > Source code: > > > > > > double pi = 3.14159265358979323846f; > > double Fs = 96000.0; > > double ln = 0.69314718055994530942f; > > double halfsq2 = 0.707106781186547524401; > > > > for(int i = 1; i < 11; i++){ > > double a0; > > double A; > > double w0; > > double f0; > > double bw = 0.0f; > > double G = 0.0f; > > int type = bandData[currentSpeakerMode][currentChannel][i].type; > > int BwQ = bandData[currentSpeakerMode][currentChannel][i].BwQNull; > > int swType; > > > > if(type < 3){ swType = 1;} // peaking > > else if(type == 3 || type == 4){ swType = 2;} //lp > > else if(type == 5 || type == 6){ swType = 3;} //hp > > else if(type == 7){ swType = 4;}//bp > > else if(type == 8 || type == 10 || type == 11){ swType = 5;}//ls > > else if(type == 9 || type == 12 || type == 13){ swType = 6;}//hs > > else if(type == 14){ swType = 7;}//notch > > else{ swType = 8;}//all pass > > > > if(bandActivity[i].Checked){ > > f0 = > > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].freqHz); > > if(bwOrQSliderNeedsToBeActive(i)){ > > bw = > > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].bwdB); > > } > > else if(type == 10 || type == 12){ bw = 6.0f;} > > else if(type == 11 || type == 13){ bw = 12.0f;} > > > > if(bw == 0.0f){ > > if(swType > 1 && swType < 5){ bw = halfsq2;} > > else if(swType == 5 || swType == 6){ bw = 0.9f;} > > else if(swType == 7){ bw = 15.0f;} > > } > > else{ > > bw /= 12.0f; > > } > > > > if(gainSliderNeedsToBeActive(i)){ > > G = > > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].gaindB); > > } > > > > if(swType == 1 || swType == 5 || swType == 6){ > > A = pow(10, (G / 40.0f)); > > } > > else{ > > A = pow(10, (G / 20.0f)); > > } > > > > w0 = 2*pi*f0/Fs; > > double alpha; > > double beta = -1.0f; > > > > if(swType == 5 || swType == 6){ // ls/hs > > alpha = sin(w0)/2 * sqrt((A + 1/A) * (1/bw - 1) + 2); > > beta = 2 * sqrt(A) * alpha; > > } > > else if(BwQ == 2){ // bw > > alpha = sin(w0) * sinh( ln * bw * w0/sin(w0) ); > > } > > else{ // q > > alpha = sin(w0) / (2 * bw); > > } > > > > ... > > switch(swType){ > > ... > > case 7: // notch > > a0 = 1 + alpha; > > coefficient[i].b0 = 1; > > coefficient[i].b1 = -2 * cos(w0); > > coefficient[i].b2 = 1; > > coefficient[i].a1 = -2 * cos(w0); > > coefficient[i].a2 = 1 - alpha; > > break; > > ... > > } > > > > // norm > > coefficient[i].b0 = coefficient[i].b0/a0; > > coefficient[i].b1 = coefficient[i].b1/a0; > > coefficient[i].b2 = coefficient[i].b2/a0; > > coefficient[i].a0 = 1.0f; > > coefficient[i].a1 = coefficient[i].a1/a0; > > coefficient[i].a2 = coefficient[i].a2/a0; > > > > .. > > > > // calculate magnitude > > numerator = coefficient[s].b0*coefficient[s].b0 + > > coefficient[s].b1*coefficient[s].b1 + coefficient[s].b2*coefficient[s].b2 + > > > > 2.0L*(coefficient[s].b0*coefficient[s].b1 + > > coefficient[s].b1*coefficient[s].b2)*cosl(w) + > > 2.0L*coefficient[s].b0*coefficient[s].b2*cosl(2.0L*w); > > denominator = 1.0L + coefficient[s].a1*coefficient[s].a1 + > > coefficient[s].a2*coefficient[s].a2 + > > 2.0L*(coefficient[s].a1 + coefficient[s].a1*coefficient[s].a2)*cosl(w) + > > 2.0L*coefficient[s].a2*cosl(2.0L*w); > > magnitude = 20*log10(sqrtl(numerator / denominator)); > > > > .. > > > > > > Is there something wrong in source code or is there other factors > > involved? > > > > _____________________________ > > Posted through www.DSPRelated.com
I havn't "dug" through your code, but when a NAN ("Not A Number") shows up, you have likely divided by zero or attempted to perform a calculation with a number not in the domain of the function causing the issue. The "zero" doesn't have to be exactly zero but can be an extremely small number resulting from the differencing of two similarly valued numbers. Since you know the frequency that makes it blow up, just trace the code as it processes the data. This should point you towards the problem. IHTH, Clay
>On 12/10/13 8:20 AM, jtp_1960 wrote: >> Hello! >> >> I have an issue with a Notch filter which works OK for certain Hz/bw/Q
but
>> after changing Hz or bw/Q sometimes a NaN value is given for the sample
for
>> center frequency. >> >> http://i39.tinypic.com/32zrgk8.jpg >> >> Source code: >> >> >> double pi = 3.14159265358979323846f; >> double Fs = 96000.0; >> double ln = 0.69314718055994530942f; >> double halfsq2 = 0.707106781186547524401; >> >> for(int i = 1; i< 11; i++){ >> double a0; >> double A; >> double w0; >> double f0; >> double bw = 0.0f; >> double G = 0.0f; >> int type = bandData[currentSpeakerMode][currentChannel][i].type; >> int BwQ = bandData[currentSpeakerMode][currentChannel][i].BwQNull; >> int swType; >> >> if(type< 3){ swType = 1;} // peaking >> else if(type == 3 || type == 4){ swType = 2;} //lp >> else if(type == 5 || type == 6){ swType = 3;} //hp >> else if(type == 7){ swType = 4;}//bp >> else if(type == 8 || type == 10 || type == 11){ swType = 5;}//ls >> else if(type == 9 || type == 12 || type == 13){ swType = 6;}//hs >> else if(type == 14){ swType = 7;}//notch >> else{ swType = 8;}//all pass >> >> if(bandActivity[i].Checked){ >> f0 = >>
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].freqHz);
>> if(bwOrQSliderNeedsToBeActive(i)){ >> bw = >>
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].bwdB);
>> } >> else if(type == 10 || type == 12){ bw = 6.0f;} >> else if(type == 11 || type == 13){ bw = 12.0f;} >> >> if(bw == 0.0f){ >> if(swType> 1&& swType< 5){ bw = halfsq2;} >> else if(swType == 5 || swType == 6){ bw = 0.9f;} >> else if(swType == 7){ bw = 15.0f;} >> } >> else{ >> bw /= 12.0f; >> } >> >> if(gainSliderNeedsToBeActive(i)){ >> G = >>
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].gaindB);
>> } >> >> if(swType == 1 || swType == 5 || swType == 6){ >> A = pow(10, (G / 40.0f)); >> } >> else{ >> A = pow(10, (G / 20.0f)); >> } >> >> w0 = 2*pi*f0/Fs; >> double alpha; >> double beta = -1.0f; >> >> if(swType == 5 || swType == 6){ // ls/hs >> alpha = sin(w0)/2 * sqrt((A + 1/A) * (1/bw - 1) + 2); >> beta = 2 * sqrt(A) * alpha; >> } >> else if(BwQ == 2){ // bw >> alpha = sin(w0) * sinh( ln * bw * w0/sin(w0) ); >> } >> else{ // q >> alpha = sin(w0) / (2 * bw); >> } >> >> ... >> switch(swType){ >> ... >> case 7: // notch >> a0 = 1 + alpha; >> coefficient[i].b0 = 1; >> coefficient[i].b1 = -2 * cos(w0); >> coefficient[i].b2 = 1; >> coefficient[i].a1 = -2 * cos(w0); >> coefficient[i].a2 = 1 - alpha; >> break; >> ... >> } >> >> // norm >> coefficient[i].b0 = coefficient[i].b0/a0; >> coefficient[i].b1 = coefficient[i].b1/a0; >> coefficient[i].b2 = coefficient[i].b2/a0; >> coefficient[i].a0 = 1.0f; >> coefficient[i].a1 = coefficient[i].a1/a0; >> coefficient[i].a2 = coefficient[i].a2/a0; >> >> .. >> >> // calculate magnitude >> numerator = coefficient[s].b0*coefficient[s].b0 + >> coefficient[s].b1*coefficient[s].b1 +
coefficient[s].b2*coefficient[s].b2 +
>> >> 2.0L*(coefficient[s].b0*coefficient[s].b1 + >> coefficient[s].b1*coefficient[s].b2)*cosl(w) + >> 2.0L*coefficient[s].b0*coefficient[s].b2*cosl(2.0L*w); >> denominator = 1.0L + coefficient[s].a1*coefficient[s].a1 + >> coefficient[s].a2*coefficient[s].a2 + >> 2.0L*(coefficient[s].a1 + coefficient[s].a1*coefficient[s].a2)*cosl(w)
+
>> 2.0L*coefficient[s].a2*cosl(2.0L*w); >> magnitude = 20*log10(sqrtl(numerator / denominator)); >> >> .. >> >> >> Is there something wrong in source code or is there other factors >> involved? >> > >exactly what variable becomes a NaN? > > >-- > >r b-j rbj@audioimagination.com > >"Imagination is more important than knowledge." > > >
39Hz/0,502 -> sample 39 value = NaN numerator=-1,33226762955019E-15 enominnator=1,4255263636187E-13 39Hz/0,503 -> sample 39 value = -23,1945375055488 numerator=6,66133814775094E-16 denominnator=1,4388490399142E-13 40Hz/0,504 -> sample 40 value = NaN numerator=-6,66133814775094E-16 denominnator=1,58317803311547E-13 _____________________________ Posted through www.DSPRelated.com
>On 12/10/13 8:20 AM, jtp_1960 wrote: >> Hello! >> >> I have an issue with a Notch filter which works OK for certain Hz/bw/Q
but
>> after changing Hz or bw/Q sometimes a NaN value is given for the sample
for
>> center frequency. >> >> http://i39.tinypic.com/32zrgk8.jpg >> >> Source code: >> >> >> double pi = 3.14159265358979323846f; >> double Fs = 96000.0; >> double ln = 0.69314718055994530942f; >> double halfsq2 = 0.707106781186547524401; >> >> for(int i = 1; i< 11; i++){ >> double a0; >> double A; >> double w0; >> double f0; >> double bw = 0.0f; >> double G = 0.0f; >> int type = bandData[currentSpeakerMode][currentChannel][i].type; >> int BwQ = bandData[currentSpeakerMode][currentChannel][i].BwQNull; >> int swType; >> >> if(type< 3){ swType = 1;} // peaking >> else if(type == 3 || type == 4){ swType = 2;} //lp >> else if(type == 5 || type == 6){ swType = 3;} //hp >> else if(type == 7){ swType = 4;}//bp >> else if(type == 8 || type == 10 || type == 11){ swType = 5;}//ls >> else if(type == 9 || type == 12 || type == 13){ swType = 6;}//hs >> else if(type == 14){ swType = 7;}//notch >> else{ swType = 8;}//all pass >> >> if(bandActivity[i].Checked){ >> f0 = >>
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].freqHz);
>> if(bwOrQSliderNeedsToBeActive(i)){ >> bw = >>
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].bwdB);
>> } >> else if(type == 10 || type == 12){ bw = 6.0f;} >> else if(type == 11 || type == 13){ bw = 12.0f;} >> >> if(bw == 0.0f){ >> if(swType> 1&& swType< 5){ bw = halfsq2;} >> else if(swType == 5 || swType == 6){ bw = 0.9f;} >> else if(swType == 7){ bw = 15.0f;} >> } >> else{ >> bw /= 12.0f; >> } >> >> if(gainSliderNeedsToBeActive(i)){ >> G = >>
System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].gaindB);
>> } >> >> if(swType == 1 || swType == 5 || swType == 6){ >> A = pow(10, (G / 40.0f)); >> } >> else{ >> A = pow(10, (G / 20.0f)); >> } >> >> w0 = 2*pi*f0/Fs; >> double alpha; >> double beta = -1.0f; >> >> if(swType == 5 || swType == 6){ // ls/hs >> alpha = sin(w0)/2 * sqrt((A + 1/A) * (1/bw - 1) + 2); >> beta = 2 * sqrt(A) * alpha; >> } >> else if(BwQ == 2){ // bw >> alpha = sin(w0) * sinh( ln * bw * w0/sin(w0) ); >> } >> else{ // q >> alpha = sin(w0) / (2 * bw); >> } >> >> ... >> switch(swType){ >> ... >> case 7: // notch >> a0 = 1 + alpha; >> coefficient[i].b0 = 1; >> coefficient[i].b1 = -2 * cos(w0); >> coefficient[i].b2 = 1; >> coefficient[i].a1 = -2 * cos(w0); >> coefficient[i].a2 = 1 - alpha; >> break; >> ... >> } >> >> // norm >> coefficient[i].b0 = coefficient[i].b0/a0; >> coefficient[i].b1 = coefficient[i].b1/a0; >> coefficient[i].b2 = coefficient[i].b2/a0; >> coefficient[i].a0 = 1.0f; >> coefficient[i].a1 = coefficient[i].a1/a0; >> coefficient[i].a2 = coefficient[i].a2/a0; >> >> .. >> >> // calculate magnitude >> numerator = coefficient[s].b0*coefficient[s].b0 + >> coefficient[s].b1*coefficient[s].b1 +
coefficient[s].b2*coefficient[s].b2 +
>> >> 2.0L*(coefficient[s].b0*coefficient[s].b1 + >> coefficient[s].b1*coefficient[s].b2)*cosl(w) + >> 2.0L*coefficient[s].b0*coefficient[s].b2*cosl(2.0L*w); >> denominator = 1.0L + coefficient[s].a1*coefficient[s].a1 + >> coefficient[s].a2*coefficient[s].a2 + >> 2.0L*(coefficient[s].a1 + coefficient[s].a1*coefficient[s].a2)*cosl(w)
+
>> 2.0L*coefficient[s].a2*cosl(2.0L*w); >> magnitude = 20*log10(sqrtl(numerator / denominator)); >> >> .. >> >> >> Is there something wrong in source code or is there other factors >> involved? >> > >exactly what variable becomes a NaN? > > >-- > >r b-j rbj@audioimagination.com > >"Imagination is more important than knowledge." > > >
39Hz/0,502 -> sample 39 value = NaN numerator=-1,33226762955019E-15 denominnator=1,4255263636187E-13 39Hz/0,503 -> sample 39 value = -23,1945375055488 numerator=6,66133814775094E-16 denominnator=1,4388490399142E-13 40Hz/0,504 -> sample 40 value = NaN numerator=-6,66133814775094E-16 denominnator=1,58317803311547E-13 In some cases there is INF instead of NaN sample 113 value = -INF numerator=0 denominnator=1,03672626039497E-12 Juha _____________________________ Posted through www.DSPRelated.com
Carification:

>39Hz/0,502 -> sample 39 value = NaN numerator=-1,33226762955019E-15 >denominnator=1,4255263636187E-13 > >39Hz/0,503 -> sample 39 value = -23,1945375055488 >numerator=6,66133814775094E-16 denominnator=1,4388490399142E-13 > >40Hz/0,504 -> sample 40 value = NaN numerator=-6,66133814775094E-16 >denominnator=1,58317803311547E-13 > >In some cases there is INF instead of NaN > >sample 113 value = -INF numerator=0 denominnator=1,03672626039497E-12
"value" is the magnitude (dB) got from 20*log10(sqrtl(numerator / denominator)). Looks like the NaN / INF is result from sqrtl() function sample 39 (39Hz) value = NaN numerator=-6,66133814775094E-16 denominnator=3,574918139293E-14 sqrtl(numerator/denominnator) = NaN Juha _____________________________ Posted through www.DSPRelated.com
On Tuesday, December 10, 2013 12:41:49 PM UTC-5, jtp_1960 wrote:
> Carification: > > > > >39Hz/0,502 -> sample 39 value = NaN numerator=-1,33226762955019E-15 > > >denominnator=1,4255263636187E-13 > > > > > >39Hz/0,503 -> sample 39 value = -23,1945375055488 > > >numerator=6,66133814775094E-16 denominnator=1,4388490399142E-13 > > > > > >40Hz/0,504 -> sample 40 value = NaN numerator=-6,66133814775094E-16 > > >denominnator=1,58317803311547E-13 > > > > > >In some cases there is INF instead of NaN > > > > > >sample 113 value = -INF numerator=0 denominnator=1,03672626039497E-12 > > > > "value" is the magnitude (dB) got from 20*log10(sqrtl(numerator / > > denominator)). > > > > Looks like the NaN / INF is result from sqrtl() function > > sample 39 (39Hz) > > value = NaN > > numerator=-6,66133814775094E-16 > > denominnator=3,574918139293E-14 > > sqrtl(numerator/denominnator) = NaN > > > > Juha > > > > > > _____________________________ > > Posted through www.DSPRelated.com
Looks like you are trying to find the square root of a negative number. The proper domain for the square root function is the non-negative numbers. Clay
"jtp_1960" <20495@dsprelated> writes:

> Hello! > > I have an issue with a Notch filter which works OK for certain Hz/bw/Q but > after changing Hz or bw/Q sometimes a NaN value is given for the sample for > center frequency. > > http://i39.tinypic.com/32zrgk8.jpg > > Source code: > > > double pi = 3.14159265358979323846f; > double Fs = 96000.0; > double ln = 0.69314718055994530942f; > double halfsq2 = 0.707106781186547524401; > > for(int i = 1; i < 11; i++){ > double a0; > double A; > double w0; > double f0; > double bw = 0.0f; > double G = 0.0f; > int type = bandData[currentSpeakerMode][currentChannel][i].type; > int BwQ = bandData[currentSpeakerMode][currentChannel][i].BwQNull; > int swType; > > if(type < 3){ swType = 1;} // peaking > else if(type == 3 || type == 4){ swType = 2;} //lp > else if(type == 5 || type == 6){ swType = 3;} //hp > else if(type == 7){ swType = 4;}//bp > else if(type == 8 || type == 10 || type == 11){ swType = 5;}//ls > else if(type == 9 || type == 12 || type == 13){ swType = 6;}//hs > else if(type == 14){ swType = 7;}//notch > else{ swType = 8;}//all pass > > if(bandActivity[i].Checked){ > f0 = > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].freqHz); > if(bwOrQSliderNeedsToBeActive(i)){ > bw = > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].bwdB); > } > else if(type == 10 || type == 12){ bw = 6.0f;} > else if(type == 11 || type == 13){ bw = 12.0f;} > > if(bw == 0.0f){ > if(swType > 1 && swType < 5){ bw = halfsq2;} > else if(swType == 5 || swType == 6){ bw = 0.9f;} > else if(swType == 7){ bw = 15.0f;} > } > else{ > bw /= 12.0f; > } > > if(gainSliderNeedsToBeActive(i)){ > G = > System::Convert::ToDouble(bandData[currentSpeakerMode][currentChannel][i].gaindB); > } > > if(swType == 1 || swType == 5 || swType == 6){ > A = pow(10, (G / 40.0f)); > } > else{ > A = pow(10, (G / 20.0f)); > } > > w0 = 2*pi*f0/Fs; > double alpha; > double beta = -1.0f; > > if(swType == 5 || swType == 6){ // ls/hs > alpha = sin(w0)/2 * sqrt((A + 1/A) * (1/bw - 1) + 2); > beta = 2 * sqrt(A) * alpha; > } > else if(BwQ == 2){ // bw > alpha = sin(w0) * sinh( ln * bw * w0/sin(w0) ); > } > else{ // q > alpha = sin(w0) / (2 * bw); > } > > ... > switch(swType){ > ... > case 7: // notch > a0 = 1 + alpha; > coefficient[i].b0 = 1; > coefficient[i].b1 = -2 * cos(w0); > coefficient[i].b2 = 1; > coefficient[i].a1 = -2 * cos(w0); > coefficient[i].a2 = 1 - alpha; > break; > ... > } > > // norm > coefficient[i].b0 = coefficient[i].b0/a0; > coefficient[i].b1 = coefficient[i].b1/a0; > coefficient[i].b2 = coefficient[i].b2/a0; > coefficient[i].a0 = 1.0f; > coefficient[i].a1 = coefficient[i].a1/a0; > coefficient[i].a2 = coefficient[i].a2/a0; > > .. > > // calculate magnitude > numerator = coefficient[s].b0*coefficient[s].b0 + > coefficient[s].b1*coefficient[s].b1 + coefficient[s].b2*coefficient[s].b2 + > > 2.0L*(coefficient[s].b0*coefficient[s].b1 + > coefficient[s].b1*coefficient[s].b2)*cosl(w) + > 2.0L*coefficient[s].b0*coefficient[s].b2*cosl(2.0L*w); > denominator = 1.0L + coefficient[s].a1*coefficient[s].a1 + > coefficient[s].a2*coefficient[s].a2 + > 2.0L*(coefficient[s].a1 + coefficient[s].a1*coefficient[s].a2)*cosl(w) + > 2.0L*coefficient[s].a2*cosl(2.0L*w); > magnitude = 20*log10(sqrtl(numerator / denominator)); > > .. > > > Is there something wrong in source code or is there other factors > involved? > > _____________________________ > Posted through www.DSPRelated.com
Hi Juha, I recently had a similar problem (I was using Octave) and it turned out my filter was unstable (poles outside the unit circle). I haven't dug through your code to check them, but you might want to. -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
clay@claysturner.com writes:

> On Tuesday, December 10, 2013 12:41:49 PM UTC-5, jtp_1960 wrote: >> Carification: >> >> >> >> >39Hz/0,502 -> sample 39 value = NaN numerator=-1,33226762955019E-15 >> >> >denominnator=1,4255263636187E-13 >> >> > >> >> >39Hz/0,503 -> sample 39 value = -23,1945375055488 >> >> >numerator=6,66133814775094E-16 denominnator=1,4388490399142E-13 >> >> > >> >> >40Hz/0,504 -> sample 40 value = NaN numerator=-6,66133814775094E-16 >> >> >denominnator=1,58317803311547E-13 >> >> > >> >> >In some cases there is INF instead of NaN >> >> > >> >> >sample 113 value = -INF numerator=0 denominnator=1,03672626039497E-12 >> >> >> >> "value" is the magnitude (dB) got from 20*log10(sqrtl(numerator / >> >> denominator)). >> >> >> >> Looks like the NaN / INF is result from sqrtl() function >> >> sample 39 (39Hz) >> >> value = NaN >> >> numerator=-6,66133814775094E-16 >> >> denominnator=3,574918139293E-14 >> >> sqrtl(numerator/denominnator) = NaN >> >> >> >> Juha >> >> >> >> >> >> _____________________________ >> >> Posted through www.DSPRelated.com > > Looks like you are trying to find the square root of a negative number. The proper domain for the square root function is the non-negative numbers. > > Clay
Wups - didn't see this. Yeah, that would do it. -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
On Tue, 10 Dec 2013 17:09:28 -0500, Randy Yates wrote:

> clay@claysturner.com writes: > >> On Tuesday, December 10, 2013 12:41:49 PM UTC-5, jtp_1960 wrote: >>> Carification: >>> >>> >>> >>> >39Hz/0,502 -> sample 39 value = NaN >>> >numerator=-1,33226762955019E-15 >>> >>> >denominnator=1,4255263636187E-13 >>> >>> >>> > >>> >39Hz/0,503 -> sample 39 value = -23,1945375055488 >>> >>> >numerator=6,66133814775094E-16 denominnator=1,4388490399142E-13 >>> >>> >>> > >>> >40Hz/0,504 -> sample 40 value = NaN >>> >numerator=-6,66133814775094E-16 >>> >>> >denominnator=1,58317803311547E-13 >>> >>> >>> > >>> >In some cases there is INF instead of NaN >>> >>> >>> > >>> >sample 113 value = -INF numerator=0 >>> >denominnator=1,03672626039497E-12 >>> >>> >>> >>> "value" is the magnitude (dB) got from 20*log10(sqrtl(numerator / >>> >>> denominator)). >>> >>> >>> >>> Looks like the NaN / INF is result from sqrtl() function >>> >>> sample 39 (39Hz) >>> >>> value = NaN >>> >>> numerator=-6,66133814775094E-16 >>> >>> denominnator=3,574918139293E-14 >>> >>> sqrtl(numerator/denominnator) = NaN >>> >>> >>> >>> Juha >>> >>> >>> >>> >>> >>> _____________________________ >>> >>> Posted through www.DSPRelated.com >> >> Looks like you are trying to find the square root of a negative number. >> The proper domain for the square root function is the non-negative >> numbers. >> >> Clay > > Wups - didn't see this. Yeah, that would do it.
When you fix that, you may also want to note that 20 * log10(sqrtl(something)) is the same as 10 * log10(something) You'll save some processing time... -- Tim Wescott Wescott Design Services http://www.wescottdesign.com