Reply by jtp_1960 December 10, 20132013-12-10
Thanks for the good hints!

After debugging the code as suggested, so far I found one bug in my
coding:

if(bw == 0.0f){
    if(swType > 1 && swType < 5){ bw = halfsq2;} // lp/hp
    else if(swType == 5 || swType == 6){ bw = 0.9f;}  // ls/hs
    else if(swType == 7){ bw = 30.0f;}  // notch
}
else{ <-- BUG (should be done for ls/hs types only)
    bw /= 12.0f;
}

This change helped for issues above 70Hz but still issues below it. I just
need to debug some more.

Juha
	 

_____________________________		
Posted through www.DSPRelated.com
Reply by Tim Wescott December 10, 20132013-12-10
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
Reply by Randy Yates December 10, 20132013-12-10
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
Reply by Randy Yates December 10, 20132013-12-10
"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
Reply by December 10, 20132013-12-10
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
Reply by jtp_1960 December 10, 20132013-12-10
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
Reply by jtp_1960 December 10, 20132013-12-10
>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
Reply by jtp_1960 December 10, 20132013-12-10
>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
Reply by December 10, 20132013-12-10
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
Reply by robert bristow-johnson December 10, 20132013-12-10
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."