Interpolated frequency sampling linear-phase FIR design

Friends,

I face the following problem. Because several vectors of differing
lenghts occur in the description, I use Harris' convention of writing
[h(n):M] to denote the vector h of length M.

Assume I specify the desired magnitude response of a linear-phase FIR
filter at N/2+1 evenly spaced points [H(n):N/2+1], starting with 0,
ending at Fs/2, where Fs is the sampling rate. Then use the inverse DFT
to get the N+1 tap FIR impulse response coefficients [h(n):N+1]. I am
now required to reduce the order N of the FIR to M < N (both M and N
are even). I do this using some (as yet unspecified) window [w(n):M+1].
This results in a new FIR with coefficients [h'(n):M+1].

I would like to use fast convolution with an N point FFT to implement
the new FIR [h'(n):M+1]. This means I have to zero pad to length N and
transform the zero padded impulse response with the DFT, let's write
[H'(n):N/2+1] to denote the first half of the DFT of h'(n).

I have the following questions:

1. Due to real-time constraints, I would like to calculate
[H'(n):N/2+1] directly from [H(n):N/2+1], without having to switch to
time domain in between. If we denote by [W(n):N/2+1] the DFT of the
zero padded window [w(n):M+1], we should have H' = H * W (* is
convolution) from the fact that h' = h . w (where . is point wise
multiplication). However, the convolution of two N/2+1 length vectors
resuluts in a length N vector, whereas I want the length of H' to be
N/2+1.

2. How should the window [w(n):M+1] be chosen, such that H' can easily
be calculated from H? I am thinking of some interpolation scheme to
compute H' from H. The problem is that H' in general has no points in
common with H, so I need to apply some smoothing to H before
interpolation. How can this be interpreted in the time domain, ie. how
can I guarantee that the impulse response of the FIR h' is indeed
limited in the range -M/2 to M/2?

Best regards,
Andor

Andor,

Responses are interspersed below

Fred

"Andor" <an2or@mailcircuit.com> wrote in message 
news:1114597721.830492.52220@g14g2000cwa.googlegroups.com...
> Friends,
>
> I face the following problem. Because several vectors of differing
> lenghts occur in the description, I use Harris' convention of writing
> [h(n):M] to denote the vector h of length M.
>
> Assume I specify the desired magnitude response of a linear-phase FIR
> filter at N/2+1 evenly spaced points [H(n):N/2+1], starting with 0,
> ending at Fs/2, where Fs is the sampling rate.

Unless you are going to do something a bit more complex, the filter has to 
be defined from zero to fs-fs/(N+1) not Fs/2.  When there are N+1 samples 
starting at zero, there will be N intervals spanned to the last sample and 
one more interval to get to fs - so N+1 intervals (assumed) between zero and 
fs.  Each of the intervals spans fs/(N+1).  Thus the above.  (But you know 
this, eh?)

>Then use the inverse DFT
> to get the N+1 tap FIR impulse response coefficients [h(n):N+1]. I am
> now required to reduce the order N of the FIR to M < N (both M and N
> are even). I do this using some (as yet unspecified) window [w(n):M+1].
> This results in a new FIR with coefficients [h'(n):M+1].

(I do wish that folks wouldn't mix "order" and "length" in the same 
discussion and I wish folks would stick with N and M for lengths!  Not a 
criticism, just a wish.)  :~)

> I would like to use fast convolution with an N point FFT to implement
> the new FIR [h'(n):M+1]. This means I have to zero pad to length N and
> transform the zero padded impulse response with the DFT, let's write
> [H'(n):N/2+1] to denote the first half of the DFT of h'(n).
>
> I have the following questions:
>
> 1. Due to real-time constraints, I would like to calculate
> [H'(n):N/2+1] directly from [H(n):N/2+1], without having to switch to
> time domain in between. If we denote by [W(n):N/2+1] the DFT of the
> zero padded window [w(n):M+1], we should have H' = H * W (* is
> convolution) from the fact that h' = h . w (where . is point wise
> multiplication). However, the convolution of two N/2+1 length vectors
> resuluts in a length N vector, whereas I want the length of H' to be
> N/2+1.

Consider the time-domain above first:
You start with N+1 points in frequency and IFFT.
You multiply with a length M+1 window where M<N, thus introducing zeros.
You FFT the result including the zeros using length N+1.
In frequency, the dual convolution is as you've observed and there is 
overlap.
In order to avoid the overlap, you have to have arrays of length M+N+2-1.
By doing the time domain multiplication (windowing to a shorter length) you 
have introduced frequency aliasing - mostly in the form of smearing.

It's curious that you would use frequency multiplication (FFT/x/IFFT) to 
apply the desired filter but are trying to avoid doing the dual process for 
generating the filter (IFFT/x/FFT).

> 2. How should the window [w(n):M+1] be chosen, such that H' can easily
> be calculated from H? I am thinking of some interpolation scheme to
> compute H' from H. The problem is that H' in general has no points in
> common with H, so I need to apply some smoothing to H before
> interpolation. How can this be interpreted in the time domain, ie. how
> can I guarantee that the impulse response of the FIR h' is indeed
> limited in the range -M/2 to M/2?

Smoothing to H before interpolation?
Interpolation is had by a smoothing process, no?
I'd just look at some transform pairs of length M+1.

One way to do what you want:
Define a rectangular window of length M+1.  The transform will be a 
Dirichlet.
*All* windows will transform into a sum of a shifted set of weighted 
Dirichlets.
The Dirichlets frequency scale depends on the temporal length.
They have periodic zeros directly related to the temporal length.
So, this may give you some insight into how to make that sum of Dirichlets 
into a "nice" convolution kernel.
[I'm thinking that having lots of zeros in the kernel may be useful].
For example:
Let's assume that M+1 = (N+1)/2.  Now the Dirichlets will have zeros at 
every other frequency sample (except at the peak).  You can construct a sum 
of shifted Dirichlets such that there are lots of zero-valued samples by 
forcing many of the shifted members to have zero weights.  Similar to a 
half-band filter's weights which have almost N/2 zeros in their coefficient 
set.

This really anwers your question because there are few "windows" that won't 
have lots of frequency domain samples - if "window" has the normal meaning / 
intent.

Fred

Fred, thanks for pitching in!

Fred Marshall wrote:

> > Assume I specify the desired magnitude response of a linear-phase
FIR
> > filter at N/2+1 evenly spaced points [H(n):N/2+1], starting with 0,
> > ending at Fs/2, where Fs is the sampling rate.
>
> Unless you are going to do something a bit more complex, the filter
has to
> be defined from zero to fs-fs/(N+1) not Fs/2.

If you don't supply a value at Fs/2, the frequency response
specification is incomplete. Consider N = 4, N/2+1 = 3. Then the
frequency resolution is (0, 1/4 Fs, 1/2 Fs, 3/4 Fs). To get a real
vector after IDFT (the FIR coefficients),
the frequency response vector must be of the form
(H(0),H(1),H(2),H(1)*), where * denotes complex conjugate, and  H(0)
and H(2) have to be real (didn't you know that? ;-).

> > I would like to use fast convolution with an N point FFT to
implement
> > the new FIR [h'(n):M+1]. This means I have to zero pad to length N
and
> > transform the zero padded impulse response with the DFT, let's
write
> > [H'(n):N/2+1] to denote the first half of the DFT of h'(n).
> >
> > I have the following questions:
> >
> > 1. Due to real-time constraints, I would like to calculate
> > [H'(n):N/2+1] directly from [H(n):N/2+1], without having to switch
to
> > time domain in between. If we denote by [W(n):N/2+1] the DFT of the
> > zero padded window [w(n):M+1], we should have H' = H * W (* is
> > convolution) from the fact that h' = h . w (where . is point wise
> > multiplication). However, the convolution of two N/2+1 length
vectors
> > resuluts in a length N vector, whereas I want the length of H' to
be
> > N/2+1.
>
> Consider the time-domain above first:
> You start with N+1 points in frequency and IFFT.
> You multiply with a length M+1 window where M<N, thus introducing
zeros.
> You FFT the result including the zeros using length N+1.
> In frequency, the dual convolution is as you've observed and there is

> overlap.
> In order to avoid the overlap, you have to have arrays of length
M+N+2-1.
> By doing the time domain multiplication (windowing to a shorter
length) you
> have introduced frequency aliasing - mostly in the form of smearing.

Agreed - for discrete vectors, the point-wise multiplication in "time
domain" is equal to the circular convolution in "frequency domain".

> It's curious that you would use frequency multiplication (FFT/x/IFFT)
to
> apply the desired filter but are trying to avoid doing the dual
process for
> generating the filter (IFFT/x/FFT).

Haha, good point! I was motivated by the fact that the DFT of a Hann
window has only three non-zero coefficients, and therefore the circular
convolution can be calculated a lot quicker than (IFFT/x/FFT). Hence my
question about a "good" window.


> > 2. How should the window [w(n):M+1] be chosen, such that H' can
easily
> > be calculated from H? I am thinking of some interpolation scheme to
> > compute H' from H. The problem is that H' in general has no points
in
> > common with H, so I need to apply some smoothing to H before
> > interpolation. How can this be interpreted in the time domain, ie.
how
> > can I guarantee that the impulse response of the FIR h' is indeed
> > limited in the range -M/2 to M/2?
>
> Smoothing to H before interpolation?

If you have an ordered set of values, smoothing means that you make the
set less wiggly (usually with lowpass filtering of some kind). OTH,
interpolation assumes that the set of values are samples of some
continuous function, and you want to approximate that function in
between the available samples.

> Interpolation is had by a smoothing process, no?

In this case, the true underlying function (the zero phase frequency
response) can be retrieved by a weighted sum of shifted Dirichlet
kernels. Interpolation is had by evaluating this function at the new
sampling grid.


> I'd just look at some transform pairs of length M+1.
>
> One way to do what you want:
> Define a rectangular window of length M+1.  The transform will be a
> Dirichlet.
> *All* windows will transform into a sum of a shifted set of weighted
> Dirichlets.
> The Dirichlets frequency scale depends on the temporal length.
> They have periodic zeros directly related to the temporal length.
> So, this may give you some insight into how to make that sum of
Dirichlets
> into a "nice" convolution kernel.

That's what I arrived at yesterday as well. By the way, what is a good
numeric approximation to the Dirichlet kernel? I tried Mathematica's
minmax, but it wouldn't converge ...

> [I'm thinking that having lots of zeros in the kernel may be useful].
> For example:
> Let's assume that M+1 = (N+1)/2.

Unfortunately, I cannot choose M - it is a given design parameter. All
I can do is choose the window [w(n):M+1].


> This really anwers your question because there are few "windows" that
won't
> have lots of frequency domain samples - if "window" has the normal
meaning /
> intent.

As I said, the DFT of the Hann window has only three non-zero
coefficients. I guess I can compute the N point DFT of the zero-padded
M point Hann window.

Thanks for the discussion,
Regards,
Andor

"Andor" <an2or@mailcircuit.com> wrote in message 
news:1114676599.443239.231560@f14g2000cwb.googlegroups.com...
> Fred, thanks for pitching in!
>
> Fred Marshall wrote:
>
>> > Assume I specify the desired magnitude response of a linear-phase
> FIR
>> > filter at N/2+1 evenly spaced points [H(n):N/2+1], starting with 0,
>> > ending at Fs/2, where Fs is the sampling rate.
>>
>> Unless you are going to do something a bit more complex, the filter
> has to
>> be defined from zero to fs-fs/(N+1) not Fs/2.
>
> If you don't supply a value at Fs/2, the frequency response
> specification is incomplete. Consider N = 4, N/2+1 = 3. Then the
> frequency resolution is (0, 1/4 Fs, 1/2 Fs, 3/4 Fs). To get a real
> vector after IDFT (the FIR coefficients),
>

Oh, I didn't mean that you didn't have to define it at fs/2.  I meant you 
had to define it all the way up to fs-fs/(N+1) and to not stop at fs/2.

> the frequency response vector must be of the form
> (H(0),H(1),H(2),H(1)*), where * denotes complex conjugate, and  H(0)
> and H(2) have to be real (didn't you know that? ;-).

Yes but I'm lazy.  That's what I meant by "unless you do something a bit 
more complex".  Saying H(3)=H(1)* is more complex :-) and it is "defining". 
I would rather talk about an IFFT on N points for purposes of discussion to 
avoid misconceptions.  Then you can do what you want in implementation.

>> One way to do what you want:
>> Define a rectangular window of length M+1.  The transform will be a
>> Dirichlet.
>> *All* windows will transform into a sum of a shifted set of weighted
>> Dirichlets.
>> The Dirichlets frequency scale depends on the temporal length.
>> They have periodic zeros directly related to the temporal length.
>> So, this may give you some insight into how to make that sum of
> Dirichlets
>> into a "nice" convolution kernel.
>
> That's what I arrived at yesterday as well. By the way, what is a good
> numeric approximation to the Dirichlet kernel? I tried Mathematica's
> minmax, but it wouldn't converge ...

I'm being lazy again....  :-)
A good numeric approximation to the Dirichlet should be:

FFT(gate) e.g. FFT([1 1 1 1 0 0 ........  0 0 1 1 1])

Once you have the central length of 1's (in this case 7), then put as many 
zeros as you like in the middle. You may find that the total N=10x or 20x 
the number of 1's is a good choice for a lot of things.  So, if there are 7 
1's as above N=140 and there are N-7=133 zeros.

>
>> [I'm thinking that having lots of zeros in the kernel may be useful].
>> For example:
>> Let's assume that M+1 = (N+1)/2.
>
> Unfortunately, I cannot choose M - it is a given design parameter. All
> I can do is choose the window [w(n):M+1].
>
>
>> This really anwers your question because there are few "windows" that
> won't
>> have lots of frequency domain samples - if "window" has the normal
> meaning /
>> intent.
>
> As I said, the DFT of the Hann window has only three non-zero
> coefficients. I guess I can compute the N point DFT of the zero-padded
> M point Hann window.

Exactly, and then it won't have only three non-zero coefficients any more! 
As above...
It's pretty easy to try it out:
Try a "window" in frequency that's 0.5 1 0.5 so to keep things real it would 
be:
[1 .5 0 0 ... 0 0 .5] defined over M+1 samples.
IFFT this window with suitable scaling.
You will get a raised cosine with peak 1.0 that is centered at zero time.
Perhaps you will want to shift this so that it starts at zero time.
Append zeros at the end to reach N+1 samples.
FFT the result.
You will see lots of nonzero samples in frequency - even though you only 
started with 3 nonzero samples.  That's because you have effectively 
interpolated the original frequency response.
And, that should be the effect of windowing N+1 samples down to M+1 
samples - there will be frequency smearing - thus nonzero samples in the FFT 
of the window.

Fred 

Fred Marshall wrote:
...
> Oh, I didn't mean that you didn't have to define it at fs/2.  I meant
you
> had to define it all the way up to fs-fs/(N+1) and to not stop at
fs/2.

Sorry, my misunderstanding. It seems to be quite hard (as noted in
other threads) to agree on one meaning of a written phrase ...

> > As I said, the DFT of the Hann window has only three non-zero
> > coefficients. I guess I can compute the N point DFT of the
zero-padded
> > M point Hann window.
>
> Exactly, and then it won't have only three non-zero coefficients any
more!
> As above...
> It's pretty easy to try it out:
> Try a "window" in frequency that's 0.5 1 0.5 so to keep things real
it would
> be:
> [1 .5 0 0 ... 0 0 .5] defined over M+1 samples.
> IFFT this window with suitable scaling.
> You will get a raised cosine with peak 1.0 that is centered at zero
time.
> Perhaps you will want to shift this so that it starts at zero time.
> Append zeros at the end to reach N+1 samples.
> FFT the result.
> You will see lots of nonzero samples in frequency - even though you
only
> started with 3 nonzero samples.  That's because you have effectively
> interpolated the original frequency response.
> And, that should be the effect of windowing N+1 samples down to M+1
> samples - there will be frequency smearing - thus nonzero samples in
the FFT
> of the window.

I think I will do it the "right" way: FFT/x/IFFT. Everything else is
approximative, and probably takes longer to compute.

However, I'm still interested in windows with few DFT coefficients
(there is another application where I might be interested in windowing
_after_ the DFT). Is there any treatsie on windows with maximum
sidelobe roll-off given some fixed number of non-zero DFT coefficients?

Regards,
Andor

I asked:

"
Is there any treatsie on windows with maximum
sidelobe roll-off given some fixed number of non-zero DFT coefficients?

"

Well, I found just the article:

A H Nuttall:
"Some Windows with Very Good Sidelobe Behaviour"
IEEE Transactions on Acoustics, Speech, and Signal Processing,
Vol. ASSP-29, No. 1, Feb. 1981

"Andor" <an2or@mailcircuit.com> wrote in message 
news:1114779891.017133.244780@o13g2000cwo.googlegroups.com...
>I asked:
>
> "
> Is there any treatsie on windows with maximum
> sidelobe roll-off given some fixed number of non-zero DFT coefficients?
>
> "
>
> Well, I found just the article:
>
> A H Nuttall:
> "Some Windows with Very Good Sidelobe Behaviour"
> IEEE Transactions on Acoustics, Speech, and Signal Processing,
> Vol. ASSP-29, No. 1, Feb. 1981
> 

"Andor" <an2or@mailcircuit.com> wrote in message 
news:1114779891.017133.244780@o13g2000cwo.googlegroups.com...
>I asked:
>
> "
> Is there any treatsie on windows with maximum
> sidelobe roll-off given some fixed number of non-zero DFT coefficients?
>
> "
>
> Well, I found just the article:
>
> A H Nuttall:
> "Some Windows with Very Good Sidelobe Behaviour"
> IEEE Transactions on Acoustics, Speech, and Signal Processing,
> Vol. ASSP-29, No. 1, Feb. 1981

Yes, you found one of the classics!  Another is:
harris, fred, "On the Use of Windows for Harmonic Analysis with the
Discrete Fourier Transform", Proceedings of the IEEE, Vol. 66 No. 1, January
1978 pp 51-83

Something I did is in:
Temes, Barcilon, Marshall "The Optimization of Bandlimited Sytems" Proc IEEE 
Vol 61 No. 2 Feb 1973 pp 196-234.

What is particularly pertinent about the last paper is that we used sincs 
(or Dirichlets) as basis functions instead of sinusoids.  Depending on how 
you look at it you may think there is really no difference - because sincs 
and Dirichlets are just sums of sinusoids - so its a mapping from one basis 
function set to another.

However, if you approximate with sincs / Dirichlets then you can decide how 
many of them will be used while keeping the full set of sinusoids.  This 
means that you can design an optimum function with 3 sincs or 5 sincs only. 
This means that you will have only 3 or 5 or .... samples in frequency that 
are nonzero.  Just what you want it seems.

I wrote a Remez-based program that approximates with sincs.  So, it can be 
done.

Here's how it works:

Let's say you want to design a lowpass filter (which is a lot like the 
kernel that you want only you want it to be quite narrow, right?)
Let's say you will use a minimax approach like the Remez algorithm.
Let's say you don't mind if the stopband / sidelobes decay as you get some 
distance from the passband - as compared to being equiripple.
You realize that sincs decay as 1/f and Dirichlets' sidelobes go flat at 
fs/2.
You realize that sinc^2 decays as 1/f^3 and Dirichlets' sidelobes get pretty 
low at fs/2 also as long as there are enough sample points in frequency.

So, imagine a lowpass filter that's made up of a contiguous sequence of 
sincs.
And, imagine that there are nonzero sincs in the stopband that's adjacent to 
the passband.
And, imagine that there are all zero sincs in the remaining stopband.
So, there's a "cluster" of sincs on the stopband and its edges and there are 
no sincs elsewhere.
The peak of each sinc is the value of the frequency sample.

In your case you want the cluster of non-zero sincs to be small.  The peak 
values are the samples you will use for convolution.  Obviously your sincs 
will have periodicity related to M+1.  So, you may still have problems 
making this work over your N+1 points.....  The method works best if the new 
window is the same length as the original.  So, in retrospect, I don't think 
this solves your problem very well.  I'm still stuck on the notion that 
you're after non-obtainium.

Fred

On 29 Apr 2005 06:04:51 -0700, "Andor" <an2or@mailcircuit.com> wrote:

>I asked:
>
>"
>Is there any treatsie on windows with maximum
>sidelobe roll-off given some fixed number of non-zero DFT coefficients?
>
>"
>
>Well, I found just the article:
>
>A H Nuttall:
>"Some Windows with Very Good Sidelobe Behaviour"
>IEEE Transactions on Acoustics, Speech, and Signal Processing,
>Vol. ASSP-29, No. 1, Feb. 1981

Hi Andor,

 If you experiment with Nuttall's equations, 
pay attention to his Eq. (5).  It means that alternating 
a_k coefficients are negative.  For example, the coefficients 
for his Eq. (30), on pp. 88, are: a_0 = 10/32, a_1 = -15/32, 
a_2 = 6/32, & a_3 = -1/32.

See ya',
[-Rick-]