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
NaN issue with NOTCH filter
Started by ●December 10, 2013
Reply by ●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."
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.comI 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 ●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/Qbut>> after changing Hz or bw/Q sometimes a NaN value is given for the samplefor>> 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 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/Qbut>> after changing Hz or bw/Q sometimes a NaN value is given for the samplefor>> 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 ●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 ●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.comLooks 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 ●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.comHi 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
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. > > ClayWups - didn't see this. Yeah, that would do it. -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●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
