Hi, As explained in Chapter 32 of dspguide, "However, it is very important to remember that the (Laplace) values in the s-plane along the y-axis are exactly equal to the Fourier transform."

If I want to filter a time series using a Laplace equation of the form wc/(s+wc), for example, how can I do it? I think I need to begin by creating a filter kernel for this filter in the time domain. Anyone know how to do that? 

My goal is to perform frequency convolution with this filter, but to avoid circular convolution, I want to zero-pad the filter kernel in the time domain appropriately (unless it can be done in the freq domain) first, then take an FFT of this zero padded filter kernel. Finally, I'd then multiply this FFT by an appropriately zero-padded signal FFT, and take an IFFT of the result.

Or, can I just compute the Laplace equation at each frequency bin using the equation wc/(s+wc), and multiply it by the signal's FFT (which I get by transforming it into the frequency domain without zero padding), then IFFT the result? (but wouldn't this have circular convolution?)

Would really (really) appreciate anyone's help on this. Thanks in advance.

• You can't filter sampled-time data with a Laplace-domain filter.  They're in two different domains.
• You can approximate sampling by that filter with sampled-time data, however.
• To avoid circular convolution you need to zero-pad your data in the time domain, then take the FFT, then do your frequency-domain magic, then do the IFFT.
• Yes, you can do a point-by-point multiplication by \(\frac{\omega_c}{j\ \omega + \omega_c}\) for \(\omega\) pertinent to each frequency bin.  Remember to get the negative frequencies, too.

Thank you Tim. 

Can you help me flush out the procedure as a list of steps? 

Suppose my time series is 8 points, all real (to keep it simple). 

1. Zero-pad the time series. But how much is needed? (Usually I can work out the required padding knowing the filter kernel length in the time domain, but I don't have that here). Should I use another 8 zero bins to produce a total of 16 time-series points?

2. Take FFT of padded (16-point) time series.

3. Isolate single-sideband FFT complex series. This will be points 1 through 9.

4. Compute the Laplace function (complex) values for each of these 9 frequency bins, and multiply it by the FFT coefficients for these 9 points.

5. Create a double-sideband FFT spectrum by copying the data from step 4 for the positive sideband (e.g. points 1 to 9), and compute points 10...16 as the complex conjugate of points 8...2, respectively (is this OK?).

6. Take IFFT of spectrum from step 5.

7. In the resulting time series, the filter ramps up in points 1 through 4, and ramps down in points 13 to 16. Discard these points, retaining the signal only in points 5 to 12. 

(1) For a filter with transfer function \( H \left ( s \right ) = \frac{\omega_c}{s + \omega_c} \), the impulse response is \(h \left ( t \right ) = \omega_c e^{-\omega_c t}\).  You should pad by three to five time constants (\(\tau = \frac{1}{\omega_c}\)).

(2) Yes.

(3-6) Yes.

(7) If everything is done right the resulting time series will not have any "pre-rise", because \( H \left ( s \right ) \) is causal.

Thanks Tim,

Regarding (1), the Laplace equation I wrote is a 1st order LPF, but in reality it could be any 1st, 2nd, or 3rd order, low-pass, high-pass, or band-pass filter. I'd like to define a default worst-case time constant based on the total signal acquisition time and the lowest-frequency of interest that I care about in my signal. For example, let's say the time difference between the first and last points in my signal (without 0-padding) is always 30 milliseconds (e.g. 33 Hz), but I only care about frequencies in the signal higher than 1 kHz (e.g. 1 ms). Can I pad an extra 5 ms for all cases and be safe? 

Regarding (7), so no pre-rise, but I still can't use the first (and last) X points, where X is half the number of points used for padding, right, because the filter isn't fully "on"?

Also, if I add a zero to the end of my time series for every point in my time series, so the last half of the padded time series are all zeros, going through this procedure results in the real-part of the filtered time series to not have a zero-padded region. There is non-zero signal in all points. Is this to be expected because the Laplace filter effects everywhere? 

There are several things that sound off, here -- starting with Laplace when you should be in the z-domain, and a transfer function of an infinite impulse response filter when you obviously want finite limits to the filtering in the time domain.  Your general method is always going to result in an infinite impulse response filter -- that's not from the fact that it's in the Laplace domain (you can find Laplace transforms of finite impulse responses, after all), but rather because of the polynomial form of the transfer function.

Rather than us (or at least me) finding ever more elaborate ways to disappoint you, why don't you tell us what you really want to do -- not "how do I do this specific little thing with Laplace", but instead "I want to filter this kind of signal to get this sort of result".

Perhaps then we can set your feet on a path of success.

Thanks Tim, I have a software application where the user enters parameters a0, a1, a2, b0, b1, b2, and b3 for the Laplace equation 

H= (b3*s^3 + b2*s^2 + b1*s + b0) / (s^3 + a2*s^2 + a1*s + a0)

in order to create a 1st, 2nd, or 3rd order filter (I don't have any flexibility on this). I don't need to use Laplace, it's just the way the data is input to me, so it seemed logical to try to compute things as described above. 

I need to use this information to filter a 30 ms time series of evenly-spaced real numbers. The signal consists of one or more spurious tones, each above 1 kHz, and no random noise. The end goal is to get a histogram of the time-domain filtered waveform. In that sense, I don't necessarily need the full 30 ms of filtered data (one third of it, say, 10 ms, would be OK, as long as it's contiguous; of course the longer the better, everything else being equal). Does that help?

That's an incredibly bad way to do it, but you're stuck with what you've got.  On the plus side, using your technique one could (I think) specify a filter that would be unstable if realized in physical components, but would just end up being noncausal in your formulation.  Go figure...

There will be issues with numerical stability.  What determines the behavior of a filter expressed that way are the roots of the numerator and denominator.  Polynomial roots get ever-more sensitive to small variations in the coefficients as the polynomial order goes up -- the rule of thumb is to break your transfer functions into a cascade of 1st- and 2nd-order transfer functions, and the only reason we use 2nd-order transfer functions is because you can't deal with resonant poles (or roots) otherwise.

If you can do anything in the UI, you might want to factor the denominator after it's entered and post some warnings if the pole locations are wonky.

Too late at night for me to be sensible, there may be more tomorrow.

>There will be issues with numerical stability.

Agreed. I've implemented checks to report potential issues to the user based on their input parameters, to ensure stability. The users work with Laplace equations often, thus the requirement. The software can use z-transforms or any other method, it's just that the user specifies the type of filtering they want with Laplace. 

So, given that you're locked into that specific way of representing things...

Yes, given your constraints, giving the filter's transients time to die down is probably a good idea.  How much time they need depends on how the filter is specified (basically, how close the filter poles approach the imaginary axis), so giving a fixed interval there depends on your users using the application wisely.

If there is a DC component to the signal (it sounds like there isn't) it would probably be wise to pad with that average, or to remove the DC component before zero-padding.  Windowing may or may not be a good idea -- it's often done with FFTs, but that's more for the purpose of analysis, not so much for filtering.  Someone else may weigh in on that point.

And yes, you could zero pad forever, and never have the response die off completely.  That's the "infinite" part of "infinite impulse response" filter.  You can only hope to pad things out enough so that the transients die off enough to not get in the way.

Thanks Tim, 

There's no DC component. 

I initially apply a window to the signal, because the signal's start and finish points don't connect, and I need to analyze the features of the FFT spectrum. When I get the filtered time trend, I divide it by the window (to undo it). That's another reason why I throw out the start and end of the filtered time series (since the window goes to zero at its edges). 

From your discussion, it sounds like it's going to be difficult to compute a fixed interval for padding for all possible filters. Perhaps a better approach is for me to fix the padding, say 5 ms on a 30 ms signal, and then specify the bandwidth (or other criteria) the software will perform well in. Working backwards, if there's 5 ms padding, we divide this by 5 time constants to get a time constant of 1 ms, which corresponds to bandwidth of 1 kHz. Can we say as long as the user equation has a bandwidth above 1 kHz, the software should work fine? Or any other meaningful way to quantify the limits of this algorithm?

I would use a flat-topped window rather than divide by the window after the fact.  Either a trapezoidal window or a raised cosine up to 1, then a raised cosine back to 0.  Then just chop off the part with the window.

Your notion of just using a constant padding and a bandwidth restriction should work.  The user needs to understand that the bandwidth restriction applies to notch and bandpass filters as well as to low-pass filters, though.

For a high-pass filter, would the low-pass cutoff frequency be the "bandwidth" of interest, that I'd want to be restrict to be above 1 kHz (using my example above)? 

For a notch (band pass?) filter, would it simply be the low-pass cutoff frequency again?

For all of them it's the absolute value of the real part of the pole.  For a high-pass filter that's roughly the corner frequency.  For a notch filter it's roughly the width of the notch.  For a bandpass filter it's roughly the width of the peak.

Thanks so much Tim, you've been a tremendous help (really).