I read chapter 16 from dspguide.com, and from that I did a
windowed-sinc function, calculated by this :
for (i=0; i<*M/2; i++)
{
h[i]=sin(2*pi*fc*(i-*M/2))/(i-*M/2);
h[*M-i-1]=h[i];
}
h[*M/2]=2*pi*fc;
the problem is, the central bin (h[*M/2]) seems to be a little less
than a billion times too big compared to the rest of the kernel, and
then, the rest of the kernel doesn't seem to have the right shape for a
low-pass FIR filter kernel. see anything that ain't right about it?
btw, here, pi=atan(1.0)*4.0, fc is the cutoff frequency expressed as a
fraction of the sampling frequency, and *M is the length of the kernel
to be generated, an always is an odd number.
Windowed-sinc function
Started by ●August 31, 2005
Reply by ●August 31, 20052005-08-31
Michel Rouzic wrote:> I read chapter 16 from dspguide.com, and from that I did a > windowed-sinc function, calculated by this : > > for (i=0; i<*M/2; i++) > { > h[i]=sin(2*pi*fc*(i-*M/2))/(i-*M/2); > h[*M-i-1]=h[i]; > } > h[*M/2]=2*pi*fc; > > the problem is, the central bin (h[*M/2]) seems to be a little less > than a billion times too big compared to the rest of the kernel, and > then, the rest of the kernel doesn't seem to have the right shape for a > low-pass FIR filter kernel. see anything that ain't right about it? > > btw, here, pi=atan(1.0)*4.0, fc is the cutoff frequency expressed as a > fraction of the sampling frequency, and *M is the length of the kernel > to be generated, an always is an odd number.Look to your indexing. Is there a divide-by-almost-zero problem? Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●August 31, 20052005-08-31
Michel Rouzic wrote:> I read chapter 16 from dspguide.com, and from that I did a > windowed-sinc function, calculated by this : > > for (i=0; i<*M/2; i++) > { > h[i]=sin(2*pi*fc*(i-*M/2))/(i-*M/2); > h[*M-i-1]=h[i]; > } > h[*M/2]=2*pi*fc; > > the problem is, the central bin (h[*M/2]) seems to be a little less > than a billion times too big compared to the rest of the kernel, and > then, the rest of the kernel doesn't seem to have the right shape for a > low-pass FIR filter kernel. see anything that ain't right about it?This is a variation of sinc(x)=sin(x)/x for x==0. You need to apply L'Hopital's rule on the expression, and see what value the expression should have for i== *M/2, and set it explicitly. By the way, you are programming in C, right? This *M/2 construct seems a bit odd. You *could* risk that the compiler divides the value in the pointer M by 2 before de-referencing it... It will look better (and probably run faster) if you make some local variable Mm = (*M)/2 before this loop. Rune
Reply by ●August 31, 20052005-08-31
"Michel Rouzic" <Michel0528@yahoo.fr> wrote in message news:1125479097.536311.73010@g47g2000cwa.googlegroups.com...>I read chapter 16 from dspguide.com, and from that I did a > windowed-sinc function, calculated by this : > > for (i=0; i<*M/2; i++) > { > h[i]=sin(2*pi*fc*(i-*M/2))/(i-*M/2); > h[*M-i-1]=h[i]; > } > h[*M/2]=2*pi*fc; > > the problem is, the central bin (h[*M/2]) seems to be a little less > than a billion times too big compared to the rest of the kernel, and > then, the rest of the kernel doesn't seem to have the right shape for a > low-pass FIR filter kernel. see anything that ain't right about it? > > btw, here, pi=atan(1.0)*4.0, fc is the cutoff frequency expressed as a > fraction of the sampling frequency, and *M is the length of the kernel > to be generated, an always is an odd number. >Is not 2*pi*fc included in the denominator "x" for sinx/x ???
Reply by ●August 31, 20052005-08-31
Michel Rouzic wrote:> I read chapter 16 from dspguide.com, and from that I did a > windowed-sinc function, calculated by this : > > for (i=0; i<*M/2; i++) > { > h[i]=sin(2*pi*fc*(i-*M/2))/(i-*M/2); > h[*M-i-1]=h[i]; > } > h[*M/2]=2*pi*fc; > > the problem is, the central bin (h[*M/2]) seems to be a little less > than a billion times too big compared to the rest of the kernel, and > then, the rest of the kernel doesn't seem to have the right shape for a > low-pass FIR filter kernel. see anything that ain't right about it?I don't think you have the right normalization (the right multiplicative constant) for the argument of sin. Imagine that fc happens to be 44100 Hz (or 8000, or 48000, or any other *integer* value). (i-*M/2) is always an integer, so you end up computing sin (2*pi*k), with k an integer value, which gives you all zeros as the result. Ok, so perhaps fc is not exactly an integer -- or it is, but the internal floating-point representation does not allow you to *exactly* represent that integer number. Then, you're getting basically "random" numbers if fc is sufficiently large -- 2*pi times a large number, you're basically getting whatever the "remainder" is -- i.e., 2*pi*K modulo 2*pi Seems like you'd have to double check your "normalization" multiplic. coefficients. HTH, Carlos --
Reply by ●August 31, 20052005-08-31
"Michel Rouzic" schrieb> I read chapter 16 from dspguide.com, > [...] > btw, here, pi=atan(1.0)*4.0, >the standard math library offers a constant for that, M_PI, IIRC. no need to calculate.> fc is the cutoff frequency expressed as a fraction of > the sampling frequency, >dspguide has it as a fraction of Nyquist, so 0<=fc<0.5, yours will be 0<=fc<1.0> and *M is the length of the kernel to be generated, > an always is an odd number. >dspguide's examples are unfortunately in BASIC (no quite suited to production code, they are also rather meant for illustration), and there the index goes (translated into C) for (i=0; i<=M; i++) instead of the C-like for (i=0; i<M; i++) Maybe this helps. Martin
Reply by ●August 31, 20052005-08-31
Martin Blume wrote:> "Michel Rouzic" schrieb > >>I read chapter 16 from dspguide.com, >>[...] >>btw, here, pi=atan(1.0)*4.0, >> > > the standard math library offers a constant for that, > M_PI, IIRC. no need to calculate.No it doesn't (unless they added it in C99, which I somewhat doubt it). M_PI is a non-standard extension in Borland compilers (not sure if others as well), NOT a stdandard C library feature. Carlos --
Reply by ●August 31, 20052005-08-31
Martin Blume wrote:> "Michel Rouzic" schrieb > >>I read chapter 16 from dspguide.com, >>[...] >>btw, here, pi=atan(1.0)*4.0, >> > > the standard math library offers a constant for that, > M_PI, IIRC. no need to calculate.That depends on the compiler since it is not part of the C standard. GCC includes it, but others do not. I include this in my headers when needed: #ifndef M_PI #define M_PI (3.1415926535897932384626433832795) #endif -- Phil Frisbie, Jr. Hawk Software http://www.hawksoft.com
Reply by ●August 31, 20052005-08-31
Michel Rouzic wrote:> I read chapter 16 from dspguide.com, and from that I did a > windowed-sinc function, calculated by this : > > for (i=0; i<*M/2; i++) > { > h[i]=sin(2*pi*fc*(i-*M/2))/(i-*M/2); > h[*M-i-1]=h[i]; > } > h[*M/2]=2*pi*fc; > > the problem is, the central bin (h[*M/2]) seems to be a little less > than a billion times too big compared to the rest of the kernel, and > then, the rest of the kernel doesn't seem to have the right shape for a > low-pass FIR filter kernel. see anything that ain't right about it?I ran into the same problem when trying to create a windowed sinc FIR generator and referencing chapter 16 of dspguide.com. Try: h[i] = sin(2*pi*fc*(i-M/2))/(pi*(i-M/2)); and h[M/2] = 2*fc I posted source code a few months back on the musicdsp.org source code archive for generating a windowed sinc FIR filter in response to a conversation there. Have a look there for more details..
Reply by ●September 1, 20052005-09-01
Jerry Avins wrote:> Michel Rouzic wrote: > > I read chapter 16 from dspguide.com, and from that I did a > > windowed-sinc function, calculated by this : > > > > for (i=3D0; i<*M/2; i++) > > { > > h[i]=3Dsin(2*pi*fc*(i-*M/2))/(i-*M/2); > > h[*M-i-1]=3Dh[i]; > > } > > h[*M/2]=3D2*pi*fc; > > > > the problem is, the central bin (h[*M/2]) seems to be a little less > > than a billion times too big compared to the rest of the kernel, and > > then, the rest of the kernel doesn't seem to have the right shape for a > > low-pass FIR filter kernel. see anything that ain't right about it? > > > > btw, here, pi=3Datan(1.0)*4.0, fc is the cutoff frequency expressed as a > > fraction of the sampling frequency, and *M is the length of the kernel > > to be generated, an always is an odd number. > > Look to your indexing. Is there a divide-by-almost-zero problem? > > Jerry > -- > Engineering is the art of making what you want from things you can get. > =AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF==AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF= =AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF I don't think so (idk). What are you thinking about?
