Frequency filtering from non stable IIR.

Hi pals, 

My post title is maybe very weird, but let me expose the situation:
I want to filter a real infinite signal x(k) by block convolution.
To do this I have the the frequency response of an IIR filter (called h).
The frequency response correspond to : H(k) = B(k)./A(k) (matlab notation
./ => pointwise division) where B(jk) is the 1024pts DFT of a sequence
b(k) which has 128pts, same for A(k).

When I plot the frequency response (DFT), plot(abs(H(k)) or freqz(b,a) i
have the desired curve...
My first idea to implement it, was to compute get the impulse response of
that IIR, then truncate it to approximate it to an FIR filter, and then
perform block convolution in freq domain with Overlap and Add method (With
the correct windows length of course ie (fft_size - fir_len + 1)).

BUT, I get stuck because the IIR filter (b,a) is not stable, 'a' has many
|roots| > 1 ... and so the impulse response is not converging to zero but
to infinite.

So my question, is first what I am missing with that method... and what is
the solution to implement such a frequency filtering ?

Thanks a lot for you help.

Nach

Nach skrev:
> Hi pals,
>
> My post title is maybe very weird, but let me expose the situation:
> I want to filter a real infinite signal x(k) by block convolution.

Do you? Are you sure you really want to stick with "infinite"?

> To do this I have the the frequency response of an IIR filter (called h).
> The frequency response correspond to : H(k) = B(k)./A(k) (matlab notation
> ./ => pointwise division) where B(jk) is the 1024pts DFT of a sequence
> b(k) which has 128pts, same for A(k).
>
> When I plot the frequency response (DFT), plot(abs(H(k)) or freqz(b,a) i
> have the desired curve...
> My first idea to implement it, was to compute get the impulse response of
> that IIR, then truncate it to approximate it to an FIR filter, and then
> perform block convolution in freq domain with Overlap and Add method (With
> the correct windows length of course ie (fft_size - fir_len + 1)).
>
> BUT, I get stuck because the IIR filter (b,a) is not stable, 'a' has many
> |roots| > 1 ... and so the impulse response is not converging to zero but
> to infinite.
>
> So my question, is first what I am missing with that method... and what is
> the solution to implement such a frequency filtering ?

First, you have not stated it explicitly, but I assume you want a
causal, stable IIR filter.

You miss the phase of the frequency response. In order to ensure
causality, you need a certain relationship betwen the imaginary and
real parts of the transfer function. It seems you don't have that, so
your IIR filter is under-specified.

That's the diagnosis, next the cure.

You have two choises: Either use the available spec as a giudeline
for the IIR filter and use one of the usual methods to design an IIR
filter that approximates the spec. Or you can use the spec as is
and design an FIR filter. 

Rune

Nach wrote:


[...]

> BUT, I get stuck because the IIR filter (b,a) is not stable, 'a' has many
> |roots| > 1 ... and so the impulse response is not converging to zero but
> to infinite.

Why would you need this? Classwork?

Here is the trick:

Despite the filter impulse response is diverging, you can design an 
artificial signal on which your filter will converge to zero. Then, you 
have to convert your data so it will be referenced to this artificial 
signal and process it with the finite impulse response of the filter.

> 
> So my question, is first what I am missing with that method... and what is
> the solution to implement such a frequency filtering ?

Hope that will impress the professor.

Vladimir Vassilevsky

DSP and Mixed Signal Design Consultant

http://www.abvolt.com

First of all, thank you for your fast reply Rune,
Let me explain a bit more clearly what I wanted to do,

>
>Nach skrev:
>> Hi pals,
>>
>> My post title is maybe very weird, but let me expose the situation:
>> I want to filter a real infinite signal x(k) by block convolution.
>
>Do you? Are you sure you really want to stick with "infinite"?
>

Hum, by real an infinite I was referring to a real time sequence, hence I
need to block-process it. We can forget infinite !

>> To do this I have the the frequency response of an IIR filter (called
h).
>> The frequency response correspond to : H(k) = B(k)./A(k) (matlab
notation
>> ./ => pointwise division) where B(jk) is the 1024pts DFT of a sequence
>> b(k) which has 128pts, same for A(k).
>>
>> When I plot the frequency response (DFT), plot(abs(H(k)) or freqz(b,a)
i
>> have the desired curve...
>> My first idea to implement it, was to compute get the impulse response
of
>> that IIR, then truncate it to approximate it to an FIR filter, and
then
>> perform block convolution in freq domain with Overlap and Add method
(With
>> the correct windows length of course ie (fft_size - fir_len + 1)).
>>
>> BUT, I get stuck because the IIR filter (b,a) is not stable, 'a' has
many
>> |roots| > 1 ... and so the impulse response is not converging to zero
but
>> to infinite.
>>
>> So my question, is first what I am missing with that method... and what
is
>> the solution to implement such a frequency filtering ?
>
>First, you have not stated it explicitly, but I assume you want a
>causal, stable IIR filter.
>

Not necessarily, as I understood block-convolution process only work with
FIR filters, hence I would better want an FIR from the frequency
specification.

>You miss the phase of the frequency response. In order to ensure
>causality, you need a certain relationship betwen the imaginary and
>real parts of the transfer function. It seems you don't have that, so
>your IIR filter is under-specified.
>

I have the relation between the imaginary and real parts of H(k) because
H(k) is computed like this:
H(k) = B(k)./A(k) 
where B(k) is the 1024pts DFT of a sequence b(k) of 128 pts.
and A(k) is the 1024pts DFT of a sequence a(k) of 128 pts.
Hence H(k) is perfectly known... but impz(b,a) lead to infinite impulse
response... so it is not stable.
The big issue is: How to filter a signal by block-convolution in frequency
domain with an aritrary frequency response (not so arbitrary because we
know it amplitude and phase/real and imaginary part)? 


>That's the diagnosis, next the cure.
>
>You have two choises: Either use the available spec as a giudeline
>for the IIR filter and use one of the usual methods to design an IIR
>filter that approximates the spec. Or you can use the spec as is
>and design an FIR filter. 
>

I would like to follow a guideline, can you tell me one ? and especially a
fast and simple method to deign FIR filter from frequency response.
(Can I simply do an ifft of my response and then truncate the result ?)

>Rune

Thanks Rune.
>
>

Nach skrev:
> First of all, thank you for your fast reply Rune,
> Let me explain a bit more clearly what I wanted to do,
>
> >
> >Nach skrev:
> >> Hi pals,
> >>
> >> My post title is maybe very weird, but let me expose the situation:
> >> I want to filter a real infinite signal x(k) by block convolution.
> >
> >Do you? Are you sure you really want to stick with "infinite"?
> >
>
> Hum, by real an infinite I was referring to a real time sequence, hence I
> need to block-process it. We can forget infinite !

Good. Everything becomes quite a bit easier that way...

> >First, you have not stated it explicitly, but I assume you want a
> >causal, stable IIR filter.
> >
>
> Not necessarily, as I understood block-convolution process only work with
> FIR filters, hence I would better want an FIR from the frequency
> specification.

OK, you relax the IIR part of the spec making everything even
easier to implement...

> >You miss the phase of the frequency response. In order to ensure
> >causality, you need a certain relationship betwen the imaginary and
> >real parts of the transfer function. It seems you don't have that, so
> >your IIR filter is under-specified.
> >
>
> I have the relation between the imaginary and real parts of H(k) because
> H(k) is computed like this:
> H(k) = B(k)./A(k)
> where B(k) is the 1024pts DFT of a sequence b(k) of 128 pts.
> and A(k) is the 1024pts DFT of a sequence a(k) of 128 pts.
> Hence H(k) is perfectly known... but impz(b,a) lead to infinite impulse
> response... so it is not stable.

Where do the sequences a(k) and b(k) come from? The spectrum
division method is generally a very bad way of designing filters.
Unless the poles of the transfer function lie strictly inside the
unit circle, the impulse response becomes unstable, as you have
learned.

> The big issue is: How to filter a signal by block-convolution in frequency
> domain with an aritrary frequency response (not so arbitrary because we
> know it amplitude and phase/real and imaginary part)?
>
>
> >That's the diagnosis, next the cure.
> >
> >You have two choises: Either use the available spec as a giudeline
> >for the IIR filter and use one of the usual methods to design an IIR
> >filter that approximates the spec. Or you can use the spec as is
> >and design an FIR filter.
> >
>
> I would like to follow a guideline, can you tell me one ?

In most DSP-related cases, one can use a design method based on
analog prototypes. The method is described in most DSP books,
look for the term "Butterworth filter".

> and especially a
> fast and simple method to deign FIR filter from frequency response.
> (Can I simply do an ifft of my response and then truncate the result ?)

Yes, you can. Make sure that no coefficients have very small
magnitudes, though. With this caveat, the backbone of the
procedure becomes:

- Choose a filter length N
- Sample the magnitude at N equidistant points in frequency domain
  such that Fs = (N+1)/N*df where df is the frequency separation
- IFFT the samples
- If necessary, shift the impulse response[*]

[*] if you use matlab, call FFTSHIFT.

And that's basically it. There may be some more details
to take care of, look for the Parks/McClellan algorithm
(also known as the Remez algorithm) in the literature.

If you have matlab and the signal processing toolbox availble,
you may want to use the Parks-McClellan function there.
It is called either FIRPM or PMFIR; I don't remember.

Rune

Rune Allnor wrote:

> Nach skrev:
...
> > I have the relation between the imaginary and real parts of H(k) because
> > H(k) is computed like this:
> > H(k) = B(k)./A(k)
> > where B(k) is the 1024pts DFT of a sequence b(k) of 128 pts.
> > and A(k) is the 1024pts DFT of a sequence a(k) of 128 pts.
> > Hence H(k) is perfectly known... but impz(b,a) lead to infinite impulse
> > response... so it is not stable.
>
> Where do the sequences a(k) and b(k) come from? The spectrum
> division method is generally a very bad way of designing filters.
> Unless the poles of the transfer function lie strictly inside the
> unit circle, the impulse response becomes unstable, as you have
> learned.

H[k] is a 1024-point vector, so the impulse response h[n] will also be
a 1024-point vector (via IDFT), and can never be unstable (as long as
A[k] =/= 0 for all k, else H[k] doesn't exist).

Where exactly is the problem?

Andor skrev:
> Rune Allnor wrote:
>
> > Nach skrev:
> ...
> > > I have the relation between the imaginary and real parts of H(k) because
> > > H(k) is computed like this:
> > > H(k) = B(k)./A(k)
> > > where B(k) is the 1024pts DFT of a sequence b(k) of 128 pts.
> > > and A(k) is the 1024pts DFT of a sequence a(k) of 128 pts.
> > > Hence H(k) is perfectly known... but impz(b,a) lead to infinite impulse
> > > response... so it is not stable.
> >
> > Where do the sequences a(k) and b(k) come from? The spectrum
> > division method is generally a very bad way of designing filters.
> > Unless the poles of the transfer function lie strictly inside the
> > unit circle, the impulse response becomes unstable, as you have
> > learned.
>
> H[k] is a 1024-point vector, so the impulse response h[n] will also be
> a 1024-point vector (via IDFT), and can never be unstable (as long as
> A[k] =/= 0 for all k, else H[k] doesn't exist).
>
> Where exactly is the problem?

The problem is a[k] (or maybe a[n] would be a better notation).
If A[k] = DFT{a[n]} then all the problems with stability and poles
and unit circles kick in.

Rune