What am I missing normalizing a FIR filter?

I analytically derive a simple raised cosine filter with a center frequency
and width by taking the inverse transform of the idealized profile in
spectral space.  The transformed analytical function is normalized so the
maximum amplitude is one.  It is then digitized and truncated at one of the
zeros.  The mean of the truncated filter is then subtracted from each
element.  I apply this filter to a triangular wave of amplitude 1 centered
about 0.5 (not zero).  The convolution (multiplying the truncated filter
length successivly along each point of the triangular wave) gives the
correct shape (it pulls out the 1, 3, 5, etc. harmonics just fine) but the
resulting amplitude is absurdly magnified!  I expected an amplitude less
than unity, but for the fundamental frequency its up near 200!

What kind of normalization am I missing?

Brian Reinhold wrote:

> I analytically derive a simple raised cosine filter with a center frequency
> and width by taking the inverse transform of the idealized profile in
> spectral space.  The transformed analytical function is normalized so the
> maximum amplitude is one.  It is then digitized and truncated at one of the
> zeros.  The mean of the truncated filter is then subtracted from each
> element.  I apply this filter to a triangular wave of amplitude 1 centered
> about 0.5 (not zero).  The convolution (multiplying the truncated filter
> length successivly along each point of the triangular wave) gives the
> correct shape (it pulls out the 1, 3, 5, etc. harmonics just fine) but the
> resulting amplitude is absurdly magnified!  I expected an amplitude less
> than unity, but for the fundamental frequency its up near 200!
> 
> What kind of normalization am I missing?

"... So that the maximum amplitude is one." Amplitude of what? The 
signal, or each coefficient? The filter's output is the sum of all the 
accumulation (sum) of all the products. (That's what "MAC" stands for.)

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

"Brian Reinhold" <breinhold@comcast.net> wrote in message
news:_J3Lb.762808$HS4.6034109@attbi_s01...
> I analytically derive a simple raised cosine filter with a center
frequency
> and width by taking the inverse transform of the idealized profile in
> spectral space.  The transformed analytical function is normalized so the
> maximum amplitude is one.  It is then digitized and truncated at one of
the
> zeros.  The mean of the truncated filter is then subtracted from each
> element.  I apply this filter to a triangular wave of amplitude 1 centered
> about 0.5 (not zero).  The convolution (multiplying the truncated filter
> length successivly along each point of the triangular wave) gives the
> correct shape (it pulls out the 1, 3, 5, etc. harmonics just fine) but the
> resulting amplitude is absurdly magnified!  I expected an amplitude less
> than unity, but for the fundamental frequency its up near 200!
>
> What kind of normalization am I missing?

Brian,

Not a very clear description.  Hard to duplicate from this.
"simple raised cosine filter" "with a center frequency and width" ...  ???
"the idealized profile in spectral space" is.. ?
"the transformed analytical function" is.. ?
"digitized" do you mean sampled?  not the same thing....
"one of the zeros" ... must refer to the vague analytical function above...

Fred

Sorry about the confusion.

The raised cosine filter is

H(w) = .5 * [1 - cos((w - w0)/wf)]  for  w0 - pi*wf < w < w0 + pi*wf,
H(w) = 0 otherwise.

where 'w0' is the center frequency
    'wf' is the "width" of the filter and
    'w' is, well, the frequency.
This function should rise continuously from 0 at 'w = w0-pi*wf', reach its
maximum at 'w = w0' and return back to zero at 'w = w0 + pi*wf'.  At least
that's the idea.  So what I want it to do is retain the 'w0' frequency and
get rid of everything else.  Of course the smaller 'wf' is the more taps I
will need on this.

The inverse transform of this function (the filter to be convolved with the
input signal) is

h(t) = cos(w0 * t) * sin(pi*wf*t) /[pi * (1 - wf*wf*t*t) * t]

The maximum amplitude of the function is 1 if I divide the whole thing by
'wf'.

So that's what I mean by the maximum amplitude of the filter is 1.

Since 'wf << w0' for my applications, a good way to truncate this filter
that gets the bulk of the meaty part is to take the zeros of 'sin(pi*wf*t).
They are given by
't = n/wf' where 'n > 1' and 'n' is an integer.  Picking 'n=2' seems to do
quite well. Then I have the range of 't' or the convolution range over which
the filter will be used. After making discrete, this gives me the number of
taps.

By 'digitizing' I should say making 'h(t)' discrete. 't' is only valid for
discrete times.

So I convolve h(nt) * x(nt) = y(nt).  Bad notation, perhaps, but you know,
for each y(nt) I have to sum up the product of x(nt) * h(nt - T) where each
sum is over the truncated range of the filter.  Lots of computations.

So what I meant by all these amplitudes is that
h(t) has a maximum value of 1 and a mean of zero over the truncated range.
x(t) is a triangle wave with a maximum value of 1 and a minimum of 0, so
it's mean over the truncated range of the filter would be .5 * (the number
of taps).  But I figured the 0 mean of the filter should take care of it.

However, even though this effectively filters out the proper harmonics of
the trianglular wave, the resulting amplitude of the components is much
greater than 1. Near 200 for the fundamental!  Is that due to the mean of
x(nt)?

In the end I wish I could find a narrow band pass filter that had very few
taps that would do the same job as this filter.  Perhaps some IIR filter?
But I have no idea how to derive one.  It is conceptually pretty easy to see
how inverse transforming a spectral window over the desired frequencies
works (even if it is a pain trying to get all my 2*pi factors, angular
versus standard frequency notations, and other normalizations correct). But
I have no idea what to base an IIR filter on.

All I want to do is grab an FSK signal out of a bunch of noise.  I thought
my integrate and dump would do the job by itself (since you are essentially
applying a matched filter to the signal) but when I simply sang into a mike
on top of the signal, the noises I could make readily defeated the FSK
detector.  Granted, my voice noise was a few orders of magnitude more
intense than the signal but it was also at much lower frequencies.  So began
the search for an input filter BEFORE the integrate and dump.

Brian

"Brian Reinhold" <breinhold@comcast.net> wrote in message
news:3pcLb.314270$_M.1815441@attbi_s54...
> Sorry about the confusion.
>
> The raised cosine filter is
>
> H(w) = .5 * [1 - cos((w - w0)/wf)]  for  w0 - pi*wf < w < w0 + pi*wf,
> H(w) = 0 otherwise.
>
> where 'w0' is the center frequency
>     'wf' is the "width" of the filter and
>     'w' is, well, the frequency.
> This function should rise continuously from 0 at 'w = w0-pi*wf', reach its
> maximum at 'w = w0' and return back to zero at 'w = w0 + pi*wf'.  At least
> that's the idea.  So what I want it to do is retain the 'w0' frequency and
> get rid of everything else.  Of course the smaller 'wf' is the more taps I
> will need on this.
>
> The inverse transform of this function (the filter to be convolved with
the
> input signal) is
>
> h(t) = cos(w0 * t) * sin(pi*wf*t) /[pi * (1 - wf*wf*t*t) * t]
>
> The maximum amplitude of the function is 1 if I divide the whole thing by
> 'wf'.
>
> So that's what I mean by the maximum amplitude of the filter is 1.
>
> Since 'wf << w0' for my applications, a good way to truncate this filter
> that gets the bulk of the meaty part is to take the zeros of
'sin(pi*wf*t).
> They are given by
> 't = n/wf' where 'n > 1' and 'n' is an integer.  Picking 'n=2' seems to do
> quite well. Then I have the range of 't' or the convolution range over
which
> the filter will be used. After making discrete, this gives me the number
of
> taps.
>
> By 'digitizing' I should say making 'h(t)' discrete. 't' is only valid for
> discrete times.
>
> So I convolve h(nt) * x(nt) = y(nt).  Bad notation, perhaps, but you know,
> for each y(nt) I have to sum up the product of x(nt) * h(nt - T) where
each
> sum is over the truncated range of the filter.  Lots of computations.
>
> So what I meant by all these amplitudes is that
> h(t) has a maximum value of 1 and a mean of zero over the truncated range.
> x(t) is a triangle wave with a maximum value of 1 and a minimum of 0, so
> it's mean over the truncated range of the filter would be .5 * (the number
> of taps).  But I figured the 0 mean of the filter should take care of it.
>
> However, even though this effectively filters out the proper harmonics of
> the trianglular wave, the resulting amplitude of the components is much
> greater than 1. Near 200 for the fundamental!  Is that due to the mean of
> x(nt)?
>
> In the end I wish I could find a narrow band pass filter that had very few
> taps that would do the same job as this filter.  Perhaps some IIR filter?
> But I have no idea how to derive one.  It is conceptually pretty easy to
see
> how inverse transforming a spectral window over the desired frequencies
> works (even if it is a pain trying to get all my 2*pi factors, angular
> versus standard frequency notations, and other normalizations correct).
But
> I have no idea what to base an IIR filter on.
>
> All I want to do is grab an FSK signal out of a bunch of noise.  I thought
> my integrate and dump would do the job by itself (since you are
essentially
> applying a matched filter to the signal) but when I simply sang into a
mike
> on top of the signal, the noises I could make readily defeated the FSK
> detector.  Granted, my voice noise was a few orders of magnitude more
> intense than the signal but it was also at much lower frequencies.  So
began
> the search for an input filter BEFORE the integrate and dump.
>
 Brian,

OK - I'm looking at it.  First off, it appears you aren't generating the
raised cosine you had in mind.

It seems to me that you want something more like:
H[w0]=1 but it looks like you have H[w0]=0
How about:
H(i) = .5 * [1 + cos(2*pi*(w(i)-w0)/wf)];
and, then you have to include the flipped values at H[-w] as well
H(N-i) = .5 * [1 + cos(2*pi*(w(i)-w0)/wf)];

Next, scaling between transforms is often a question.  You started with a
filter with peak value of 1.0 in frequency (which is how I prefer to scale
things usually).  Then, when you divided by wf to get a peak value of 1.0 in
time, you scaled up the result by wf!  So, that accounts for the large
scaling of the result.

Then, you find that you have to truncate the result in the time domain
because you started with a finite-length spectral character in the frequency
domain.  Doing this changes the raised cosine in frequency.

You'd do better to do this: Design a FIR bandpass filter with the length you
want and the bandpass character you want.  There are quite a few programs on
the web that will do this for you.  It's a whole lot easier and gets rid of
a lot of the steps you're needing to take with the process you've defined.
Here's one of many:

http://www.systolix.co.uk/feintro.htm

Using this or some other design program you'll find that a decent FIR filter
will be of length 100 or more depending on your spec's.  On the other hand
an IIR filter might be done in 3 biquads or 15 coefficients.

Fred

Thanks Fred.

I made  a typo.  I do exactly as you suggested (including the mirror profile
to make the final form real) except that I don't have the 2*pi factor. I put
that in later.

But I shouldn't normalize the final real expression such that it's maximum
value is one (by dividing by 'wf')?  What I did find was that I got a much
closer answer if I normalized the convolution by the integral of  the
truncated version of h(t)*h(t). But the  answer was still off the
theoretical.  I will try it without first normalizing by 'wf' and see if
that works.

In the end a filter with that many taps will kill me.  I guess I'd
betterlook for some IIR stuff.

Brian
>
> OK - I'm looking at it.  First off, it appears you aren't generating the
> raised cosine you had in mind.
>
> It seems to me that you want something more like:
> H[w0]=1 but it looks like you have H[w0]=0
> How about:
> H(i) = .5 * [1 + cos(2*pi*(w(i)-w0)/wf)];
> and, then you have to include the flipped values at H[-w] as well
> H(N-i) = .5 * [1 + cos(2*pi*(w(i)-w0)/wf)];
>
> Next, scaling between transforms is often a question.  You started with a
> filter with peak value of 1.0 in frequency (which is how I prefer to scale
> things usually).  Then, when you divided by wf to get a peak value of 1.0
in
> time, you scaled up the result by wf!  So, that accounts for the large
> scaling of the result.
>
> Then, you find that you have to truncate the result in the time domain
> because you started with a finite-length spectral character in the
frequency
> domain.  Doing this changes the raised cosine in frequency.
>
> You'd do better to do this: Design a FIR bandpass filter with the length
you
> want and the bandpass character you want.  There are quite a few programs
on
> the web that will do this for you.  It's a whole lot easier and gets rid
of
> a lot of the steps you're needing to take with the process you've defined.
> Here's one of many:
>
> http://www.systolix.co.uk/feintro.htm
>
> Using this or some other design program you'll find that a decent FIR
filter
> will be of length 100 or more depending on your spec's.  On the other hand
> an IIR filter might be done in 3 biquads or 15 coefficients.
>
> Fred
>
>

"Brian Reinhold" <breinhold@comcast.net> wrote in message
news:ITwLb.5510$5V2.10651@attbi_s53...
> Thanks Fred.
>
> I made  a typo.  I do exactly as you suggested (including the mirror
profile
> to make the final form real) except that I don't have the 2*pi factor. I
put
> that in later.
>
> But I shouldn't normalize the final real expression such that it's maximum
> value is one (by dividing by 'wf')?  What I did find was that I got a much
> closer answer if I normalized the convolution by the integral of  the
> truncated version of h(t)*h(t). But the  answer was still off the
> theoretical.  I will try it without first normalizing by 'wf' and see if
> that works.
>
> In the end a filter with that many taps will kill me.  I guess I'd
> betterlook for some IIR stuff.


Brian,

h(t)*h(t) ?????  what's that about?

> But I shouldn't normalize the final real expression such that it's maximum
> value is one (by dividing by 'wf')?

No, you shouldn't.  Consider this:

Start with a temporal cosine of peak amplitude 1.0 with an integral number
of cycles in 256 seconds- period of some integer multiple of 1/256Hz like
10/256.
Window it (multiply it) with a rectangular window of length 256 and
amplitude 1.0 - from t=0 to t=255.
This does nothing to the original samples within that time span.

Now, let's do the corresponding dual operation in the frequency domain.
FFT the window.  You get a single sample at index 1 (f=0) with amplitude
256 - not 1.0
FFT the cosine.  You get a pair of samples at indices 11 and 247 with
amplitudes of 128 each - not 1.0
Convolve the filter with the cosine in frequency.  You get a pair of samples
of amplitude 128*256.
You can either look at this and decide to scale by 1/256 right away or
inverse transform the result with scaling of 1/256 in *addition* to the
built in scaling typical in an IFFT.
You get a cosine in time of amplitude 1.0 - which is what you better get!

The bottom line:
If you start with a signal or function that has peak amplitude of 1.0 in one
domain then it is unlikely to have peak amplitude of 1.0 in the dual domain.
You generally can't have both.  You defined the raised cosine to have a peak
of 1.0 in frequency - so don't expect the peak in time to also be 1.0.
Since the dual transform pair scales together, dividing by wf divides in
both domains so the 1.0 cosine you started with is thus perturbed by the
same factor.

Fred

THanks for all your comments and help Fred.  It has been very helpful for an
unseasoned DSP person like me.

I realize there is an amplitude normalization problem, just as there is with
fourier integrals.  So I tried to normalize h(t) that I had derived
analytically such that if I convolved it with itself, the maximum amplitude
of the convolution (which would occur when the two 'functions' were exactly
in phase) would be 1.  The scaling factor would be computed by
Integral[h(t)*h(t)dt] from plus or minus infinity (where '*' means product
and not convolution.  Of course I truncate the range to that of the
truncated function. I removed the 'wf' normalization.  In this way if the
signal was identical to h(t), then it should be 'passed through' perfectly.
At least that was the idea.  It did not work unless I normaized by 'wf'
first, however, but I don't rule out a program bug yet.

It's probably something stupid that I am doing in the code, because every
other aspect of the filter behaves as anticipated!

Brian


"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:6NydnXlCNfUE0Z3dRVn-jg@centurytel.net...
>
> "Brian Reinhold" <breinhold@comcast.net> wrote in message
> news:ITwLb.5510$5V2.10651@attbi_s53...
> > Thanks Fred.
> >
> > I made  a typo.  I do exactly as you suggested (including the mirror
> profile
> > to make the final form real) except that I don't have the 2*pi factor. I
> put
> > that in later.
> >
> > But I shouldn't normalize the final real expression such that it's
maximum
> > value is one (by dividing by 'wf')?  What I did find was that I got a
much
> > closer answer if I normalized the convolution by the integral of  the
> > truncated version of h(t)*h(t). But the  answer was still off the
> > theoretical.  I will try it without first normalizing by 'wf' and see if
> > that works.
> >
> > In the end a filter with that many taps will kill me.  I guess I'd
> > betterlook for some IIR stuff.
>
>
> Brian,
>
> h(t)*h(t) ?????  what's that about?
>
> > But I shouldn't normalize the final real expression such that it's
maximum
> > value is one (by dividing by 'wf')?
>
> No, you shouldn't.  Consider this:
>
> Start with a temporal cosine of peak amplitude 1.0 with an integral number
> of cycles in 256 seconds- period of some integer multiple of 1/256Hz like
> 10/256.
> Window it (multiply it) with a rectangular window of length 256 and
> amplitude 1.0 - from t=0 to t=255.
> This does nothing to the original samples within that time span.
>
> Now, let's do the corresponding dual operation in the frequency domain.
> FFT the window.  You get a single sample at index 1 (f=0) with amplitude
> 256 - not 1.0
> FFT the cosine.  You get a pair of samples at indices 11 and 247 with
> amplitudes of 128 each - not 1.0
> Convolve the filter with the cosine in frequency.  You get a pair of
samples
> of amplitude 128*256.
> You can either look at this and decide to scale by 1/256 right away or
> inverse transform the result with scaling of 1/256 in *addition* to the
> built in scaling typical in an IFFT.
> You get a cosine in time of amplitude 1.0 - which is what you better get!
>
> The bottom line:
> If you start with a signal or function that has peak amplitude of 1.0 in
one
> domain then it is unlikely to have peak amplitude of 1.0 in the dual
domain.
> You generally can't have both.  You defined the raised cosine to have a
peak
> of 1.0 in frequency - so don't expect the peak in time to also be 1.0.
> Since the dual transform pair scales together, dividing by wf divides in
> both domains so the 1.0 cosine you started with is thus perturbed by the
> same factor.
>
> Fred
>
>
>
>
>
>

"Brian Reinhold" <breinhold@comcast.net> wrote in message
news:wb2Mb.19459$nt4.41120@attbi_s51...
> THanks for all your comments and help Fred.  It has been very helpful for
an
> unseasoned DSP person like me.
>
> I realize there is an amplitude normalization problem, just as there is
with
> fourier integrals.  So I tried to normalize h(t) that I had derived
> analytically such that if I convolved it with itself, the maximum
amplitude
> of the convolution (which would occur when the two 'functions' were
exactly
> in phase) would be 1.

I don't see why normalizing that way makes sense.

For example, if you convolve 1 2 3 with itself, the peak result is 10.  It
should be.  If you convolve 1 2 3 with 4 4 4, the result has a peak of 24.
Now, the convolution of 1 2 3 and 4 4 4 should not be 24/10 = 2.4  - there's
no reason I can think of for this normalization.  Maybe it's just that I
don't know of one.

You may need to differentiate between continuous, infinite analytical
expressions and discrete, finite or possibly periodic expressions.

Here is your original post heavily edited with changes made so the process
is correct:

I analytically derive a simple raised cosine filter with peak 1.0 and a
center frequency
and width according to my specifications.  This is a bandpass filter with
band center gain of 1.0

I get the time domain version by taking the inverse transform of this
idealized spectral function.  The IFFT that I use scales by 1/N internally
as is the usual thing (but not what FFTW does).

The transformed temporal function or impulse response has been obviously
pre-normalized so the maximum amplitude if the filter at band center is 1.0,
so no further normalization is necessary.

The continuous, infinite impulse response is then sampled and truncated.
Just to be sure the stopband at f=0 has zero response, the very small mean
of the truncated filter is then subtracted from each filter coefficient
divided by the length of the filter -> i.e. subtracting mean/N.

I apply this filter to a cosine wave of amplitude 1 centered at amplitude
1.0 with frequency 0.25 (the center of the bandpass) - so it is sampled
exactly at 4 samples per period and is indistinguishable from a triangular
wave which would be undersampled.

The convolution (multiplying the truncated filter length successivly along
each point of the triangular wave) gives the correct temporal shape of a
cosine sampled 4 times per cycle.  Obviously no harmonics can be seen
because the sample rate doesn't allow.  The resulting amplitude is +1 to -1
at the peaks as expected.  The average value of 1.0 of the input is
effectively filtered out.

Now, if the bandpass is narrow enough to attenuate the harmonics, even if we
move the bandpass to a lower frequency so the relative sample rate is
higher, no harmonics will emerge the bandpass.

Fred

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:_L6dnSpEcJxclpzdRVn-hw@centurytel.net...
>

Brian,

I want to modify this one thing that I said:

> The continuous, infinite impulse response is then sampled and truncated.
> Just to be sure the stopband at f=0 has zero response, the very small mean
> of the truncated filter is then subtracted from each filter coefficient
> divided by the length of the filter -> i.e. subtracting mean/N.

It would have been correct to say "subtracting the mean" if the objective
was to set the response at f=0 identically to zero.

However, I think this is an unusual thing to to and isn't guaranteed to work
in the general case.
What if the filter unit sample response is 1 1 1 ? Then the mean is 1 and
subtracting the mean yields all zero coefficients.  I rather think that the
response at f=0 is the same sort of thing as the response at any other
stopband frequency.  So, you just accept it along with all the other
non-zero responses in the stop band and don't subtract the mean.  In the
case you have, it probably doesn't hurt either - the "correction" should be
pretty small.

Fred