Sampling a transfer function to insure causality

Suppose I have a continuous frequency transfer function H(w).  This
transfer function has an impulse response function, h(t), that is
continuous, real, and most importantly causal.

Now suppose I sample H(w) using N points to form a discrete transfer
function H(n).  I then invert H(n) to find the discrete impulse
response function h(n).

The problem is that h(n) is no longer causal.  There typically seems to
be some form of Gibbs like effect at the end of the h(n).  Sometimes
this effect contains a significant amount of energy, so I worry about
truncating it.  Increasing the number of sample points doesn't help.

Is there a way to pad or window H(n), or sample H(w) differently to
eliminate this problem?

Thanks in advance!

cdunn@metisdesign.com wrote:

   ...

> Now suppose I sample H(w) using N points to form a discrete transfer
> function H(n).  I then invert H(n) to find the discrete impulse
> response function h(n).

Does inverting H(n) really do that?

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

cdunn@metisdesign.com wrote:

> Suppose I have a continuous frequency transfer function H(w).  This
> transfer function has an impulse response function, h(t), that is
> continuous, real, and most importantly causal.
>
> Now suppose I sample H(w) using N points to form a discrete transfer
> function H(n).  I then invert H(n) to find the discrete impulse
> response function h(n).
>
> The problem is that h(n) is no longer causal.  There typically seems to
> be some form of Gibbs like effect at the end of the h(n).

I think the ringing you are seeing is the symmetric main impulse of the
filter. The reason that it is at the ends of the filter is because you
computed the inverse of zero-phase frequency response, instead of the
linear-phase. You can easily correct this by cyclically rotating the
vector h[n] by half its length (or by changing the sign of every second
sample of the frequency response before IDFT).

> Sometimes
> this effect contains a significant amount of energy, so I worry about
> truncating it.  Increasing the number of sample points doesn't help.

Yup. The main impulse.

>
> Is there a way to pad or window H(n), or sample H(w) differently to
> eliminate this problem?

Windowing comes afterwards, to reduce the real Gibbs effect (in the
frequency domain) that comes from truncating the Fourier series of the
frequency response (ie, the impulse response). You typically sample the
frequency response at a very large number, say N, and then, after IDFT,
you window the impulse response to the required filter length M. For
googling, this is called the windowing technique of FIR filter design.

Regards,
Andor

Sorry I simplified my description somewhat.  I actually sample H(w) so
that h(n) is real.  This is done in the following manner:

H(n) = H(w = w(n))                               n = 0 to N / 2 - 1
       = Conjugate[H(n = N - i)]               n = N / 2 + 1 to N - 1

Where w(i) is the i th radial frequency data point

I'm not quite sure what to with the data point at n = N / 2.  To
insure h(n) is real, H(n = N / 2) must be real.  I have tried the
following methods:

1.	Set H(n = N / 2) to zero
2.	Set H(n = N / 2) to the real part of H(w = w(N / 2)
3.	Set H(n = N / 2) to the absolute value of H(w = w(N / 2)
4.	Window H(n) such that  it smoothly goes to zero at n = N / 2 when
approached from either side

None of these methods have helped remove the non-causality in h(n).

Thanks for the question,
Chris

cdunn@metisdesign.com wrote:
> Sorry I simplified my description somewhat.  I actually sample H(w) so
> that h(n) is real.  This is done in the following manner:
>
> H(n) = H(w = w(n))                               n = 0 to N / 2 - 1
>        = Conjugate[H(n = N - i)]               n = N / 2 + 1 to N - 1
>
> Where w(i) is the i th radial frequency data point
>
> I'm not quite sure what to with the data point at n = N / 2.  To
> insure h(n) is real, H(n = N / 2) must be real.  I have tried the
> following methods:
>
> 1.	Set H(n = N / 2) to zero
> 2.	Set H(n = N / 2) to the real part of H(w = w(N / 2)

Try twice the real part, i.e., the sum of a conjugate pair of
coefficients is twice the real.

> 3.	Set H(n = N / 2) to the absolute value of H(w = w(N / 2)
> 4.	Window H(n) such that  it smoothly goes to zero at n = N / 2 when
> approached from either side
>
> None of these methods have helped remove the non-causality in h(n).

You say you have a frequency response that is causal.  There are a
couple of possibilities tat come to mind.

First, what is true for a continuous function is not necessarily true
for a sampled version of that function.  This is similar to the fact
that Fourier transform pairs are not discrete Fourier transform pairs,
though, if the sampling rate is high enough, they are close enough.

Second, examine the solution in the time domain.  If it extends off the
end of the array, then that could be the problem.  FFT is misnamed, in
my opinion, it should be called FFS, for fast Fourier series.
Functions are always periodic.  You may need a wider array, so that the
time response has enough time to fall to near zero.

Third, a causal response has real and imaginary parts that are related
by the Hilbert transform.  You can always write in the real part and
use a fast Hilbert transform to get the imaginary part, so it frequency
response would have to be causal.

Fourth, write the response in the time domain, making sure it has
decayed to near zero in the time period.

I hope this helps.  Have fun.

Hi,

I would like to add sixth: Please read and understand Andors answer.
That should clarify the main topics. A nice introduction can be found
in Smith's "The Scientist and Engineers Guide to DSP", page 299.

> Fourth, write the response in the time domain, making sure it has
> decayed to near zero in the time period.

To rtonyreeder:
I do not argue your answer, just the question: How do you decide that
the time response did decay close enough to zero? Is -60dB enough? -80?
-100dB? ...?

stereo

> First, what is true for a continuous function is not necessarily true
> for a sampled version of that function.  This is similar to the fact
> that Fourier transform pairs are not discrete Fourier transform pairs,
> though, if the sampling rate is high enough, they are close enough.

If the sampling frequency is high enough, it shouldn't matter wheather
it is a discrete or continuous transform pair. Since, if the Nyquist
criterion is fulfilled, the sampled description of the (bandlimited)
signal in time and frequency >>exactly<< represent the continuous
signal.
In my opinion. (?)

regards

stereo

stereo wrote:
> To rtonyreeder:
> I do not argue your answer, just the question: How do you decide that
> the time response did decay close enough to zero? Is -60dB enough? -80?
> -100dB? ...?

You have to decide what is close enough to zero, i.e., how much it is
affecting your answer.

I didn't go into how to make a periodic (cyclic) convolution give the
same answer as non-periodic (linear) convolution.  You don't have to be
concerned if the time domain kernel (impulse response) has decayed to
zero, if you are only interested in the answer over a certain temporal
range.  You have to pad with zeros properly.

stereo also wrote:
> > First, what is true for a continuous function is not necessarily true
> > for a sampled version of that function.  This is similar to the fact
> > that Fourier transform pairs are not discrete Fourier transform pairs,
> > though, if the sampling rate is high enough, they are close enough.
>
> If the sampling frequency is high enough, it shouldn't matter wheather
> it is a discrete or continuous transform pair. Since, if the Nyquist
> criterion is fulfilled, the sampled description of the (bandlimited)
> signal in time and frequency >>exactly<< represent the continuous
> signal.
> In my opinion. (?)

Real responses don't decay all that fast in the frequency domain.
Gaussians, Lorentzian's, filters, etc. not analytically bandlimited.
It just depends on how much of an error is caused by aliasing low
amplitude high frequencies.  If you sample at a high enough rate, you
get a representation that is close enough to the continuous function.

I said earlier.
> Fourth, write the response in the time domain, making sure it has
> decayed to near zero in the time period.

And.
> I didn't go into how to make a periodic (cyclic) convolution give the
> same answer as non-periodic (linear) convolution.  You don't have to be
> concerned if the time domain kernel (impulse response) has decayed to
> zero, if you are only interested in the answer over a certain temporal
> range.  You have to pad with zeros properly.

I should have indicated that making the impulse go to near zero is not
enough, or even necessary, depending on what you are looking for.  You
could still be getting cyclic convolutions rather than linear
convolutions if you don't zero pad properly.  Of course, the sampling
rate is still important.  The second quote is correct.

Generally, linear convolution of an input with an impulse response
requires an output array that is equal in length to the sum of the
lengths of the two arrays being convolved.  Zero padding allows all of
the arrays to be the same length.

I hope this isn't too confusing.

On 26 Jun 2006 11:51:37 -0700, cdunn@metisdesign.com wrote:

>Suppose I have a continuous frequency transfer function H(w).  This
>transfer function has an impulse response function, h(t), that is
>continuous, real, and most importantly causal.

Hi cdunn,

  ah, you're entering the "Twilight Zone" of 
considering the differences between continuous 
and discrete signals. 

>Now suppose I sample H(w) using N points to form a discrete transfer
>function H(n).  I then invert H(n) to find the discrete impulse
>response function h(n).
>
>The problem is that h(n) is no longer causal. 
 
Humm, ... I wonder exactly what you mean by 
"no longer causal".  Do you mean that the 
time-domain (integer) index value of zero is 
in the center of your impulse response.

>There typically seems to
>be some form of Gibbs like effect at the end of the h(n).  Sometimes
>this effect contains a significant amount of energy, so I worry about
>truncating it.  Increasing the number of sample points doesn't help.
>
>Is there a way to pad or window H(n), or sample H(w) differently to
>eliminate this problem?

When we perform a Fourier transform (or an inverse
Fourier transform) of a continuous function that has 
very abrupt changes (discontinuities), the transform 
results show ripples (this rippling effect is called the 
"Gibbs phenomenon").

Does your H(w) have very abrupt changes discontinuities)?

[And no, regardless of what you read in DSP books,
the "Gibbs phenomenon" is *NOT* named after 
Barry Gibb of the rock-n-roll group The GeeGees.]

Good Luck,
[-Rick-]