bilinear transform

Hi again,

Thank's for the help and the good book tips.
Here is my new problem:
My goal is to create an inverse filter with which I'd like to filter my
digital seismograms. With the inverse filtering I like to remove the effects
of the seismometer's transfer function on the seismogram. The transfer
function of the seismometer is given as analog poles and zeros.
The transfer function has got a lower cutoff frequency of 1/30 Hertz and an
upper cutoff frequency of 50 Hz. My seismograms are sampled at 100 Hz.

I have managed to transform the analog bandpass filter to an according
digital filter.
It works quite well, but the filter slope is too steep at the upper cutoff
frequency.
I now that this is because the upper cutoff frequency is at 50Hz which is
the Nyquist frequency for my choosen sampling frequency of 100 Hz, but is
ther a possibility that I can avoid this steep slope?


Here is my code:
---------------------------------
z = [0; 0];
p = [-0.02365+0.02365i; -0.02365-0.02365i; -180; -160; -80];
k = 2304000;

% Conversion from Hertz to Rad/sec.
z = z*2*pi;
p = p*2*pi;
k = k*(2*pi)^(length(p)-length(z));

fs = 100;       % Time series sampling frequency


% Display of the analog filter
[b, a] = zp2tf(z, p, k);
[h, w] = freqs(b,a);
semilogx(w/(2*pi), 20*log10(abs(h)));
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
title('blue: analog filter; red: digital filter using bilinear transform');
grid on;
hold on;

% Analog to digital conversion using bilinear transform
[b_d, a_d] = bilinear(b, a, fs, 1/30);
[h, f] = freqz(b_d, a_d, 100000, fs);
semilogx(f, 20*log10(abs(h)), 'color', 'r');
grid on;
hold off;
---------------------------------------

Hi.

It depends. If you use a sampling frequency of 100Hz for a bandpass
filter with a upper cutoff frequency of 50Hz, this makes no sense.
Consider using a highpass filter that uses the lower cutoff frequency
of your bandpass and no upper cutoff.
Otherwise you might always consider using a higher sampling
frequency...

Hope this helps, Heiko

mankei wrote:
> Hi again,
>
> Thank's for the help and the good book tips.
> Here is my new problem:
> My goal is to create an inverse filter with which I'd like to filter my
> digital seismograms. With the inverse filtering I like to remove the effec
ts
> of the seismometer's transfer function on the seismogram. The transfer
> function of the seismometer is given as analog poles and zeros.
> The transfer function has got a lower cutoff frequency of 1/30 Hertz and a
n
> upper cutoff frequency of 50 Hz. My seismograms are sampled at 100 Hz.

You are in trouble! You mentin below that the 100 Hz sampling frequency

is "chosen" to be 100 Hz, and 100 Hz seem a bit low for seismological
work. If you can, resample your data to a higher sampling frequency.
As far as I know, seismometers these days use sampling rates on the
order of milliseconds. You need that oversampling to clean up your
data.

> I have managed to transform the analog bandpass filter to an according
> digital filter.
> It works quite well, but the filter slope is too steep at the upper cutoff
> frequency.

The inverse filter? How did you do that?

> I now that this is because the upper cutoff frequency is at 50Hz which is
> the Nyquist frequency for my choosen sampling frequency of 100 Hz, but is
> ther a possibility that I can avoid this steep slope?

What you need to do, is to transform the seismometer's transfer
function
into z domain before starting playing with the inverse filter. There
is a non-linear frequency mapping between s domain and z domain.
I suspect you will find that the stop-band slope of your z domain
seismometer transfer function is a lot steeper than you would expect
it to be.

But all that is irrelevant unless you get data at a much higher
sampling
frequency.

Rune

Hi Rune,

100 Hz sampling frequency is ok for my needs. The seismometer's transfer
function is roughly a bandpass filter from 0.0333 Hz to 50 Hz. So any
frequency content of higher frequencies than 50 Hz is not interesting for
me.

I have posted the code how I calculate the seismometer's transfer function
in the z-domain. I'm doing that with the bilinear transform. And yes, the
filter roll off at 50 Hz is very steep.

My question is, if there is a way to avoid this steep slope at the upper
cutoff frequency?

thanks,
Stefan.


"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:1126087010.761748.213040@g43g2000cwa.googlegroups.com...
>
> mankei wrote: 
effects 
an 
>
> You are in trouble! You mentin below that the 100 Hz sampling frequency
>
> is "chosen" to be 100 Hz, and 100 Hz seem a bit low for seismological
> work. If you can, resample your data to a higher sampling frequency.
> As far as I know, seismometers these days use sampling rates on the
> order of milliseconds. You need that oversampling to clean up your
> data.
> 
cutoff 
>
> The inverse filter? How did you do that?
> 
is 
is 
>
> What you need to do, is to transform the seismometer's transfer
> function
> into z domain before starting playing with the inverse filter. There
> is a non-linear frequency mapping between s domain and z domain.
> I suspect you will find that the stop-band slope of your z domain
> seismometer transfer function is a lot steeper than you would expect
> it to be.
>
> But all that is irrelevant unless you get data at a much higher
> sampling
> frequency.
>
> Rune
>

mankei wrote:
> Hi Rune,
>
> 100 Hz sampling frequency is ok for my needs.

No, it is not.

> The seismometer's transfer
> function is roughly a bandpass filter from 0.0333 Hz to 50 Hz. So any
> frequency content of higher frequencies than 50 Hz is not interesting for
> me.

You need that leeway to get anything useful from the data, if the
original data are sampled at a higher rate. If the ADC actually
samples at 100 Hz, the system is very badly designed.

> I have posted the code how I calculate the seismometer's transfer function
> in the z-domain. I'm doing that with the bilinear transform. And yes, the
> filter roll off at 50 Hz is very steep.
>
> My question is, if there is a way to avoid this steep slope at the upper
> cutoff frequency?

One trivial way would be to transform the 50 Hz limit to a lower
relative frequency. That would involve resampling the data to
a higher sampling frequency, assuming the data was sampled by
a properly designed system in the first place. Read up on "Frequency
Warping" in your text on the BLT.

Rune

Hi Rune,

The seismometer's ADC actually samples at 2000 Hz. The data is than filtered
and decimated by a dsp software to provide sample rates from 200 Hz
downwards. This reduction of the sample rate is needed to reduce the
filesize. The seismometer is a standalone instrument and stores the data in
a flash memory, therefore filesize is a big issue.

I think I will go for resampling to avoid the steep slope. Have you got any
experiences with Matlab's resample function? Are there any bugs?

Stefan.

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:1126100879.770334.232790@g44g2000cwa.googlegroups.com...
>
> mankei wrote: 
>
> No, it is not.
> 
for 
>
> You need that leeway to get anything useful from the data, if the
> original data are sampled at a higher rate. If the ADC actually
> samples at 100 Hz, the system is very badly designed.
> 
function 
the 
>
> One trivial way would be to transform the 50 Hz limit to a lower
> relative frequency. That would involve resampling the data to
> a higher sampling frequency, assuming the data was sampled by
> a properly designed system in the first place. Read up on "Frequency
> Warping" in your text on the BLT.
>
> Rune
>

mankei wrote:
> Hi Rune,
>
> The seismometer's ADC actually samples at 2000 Hz. The data is than filter
ed
> and decimated by a dsp software to provide sample rates from 200 Hz
> downwards. This reduction of the sample rate is needed to reduce the
> filesize. The seismometer is a standalone instrument and stores the data i
n
> a flash memory, therefore filesize is a big issue.
>
> I think I will go for resampling to avoid the steep slope. Have you got an
y
> experiences with Matlab's resample function? Are there any bugs?

There is little you can do with software, once the data are
decimated to 100 Hz sampling frequency. You need to access the
data files from the seismometer at an as high sampling frequency
as possible, preferably 200 Hz.

Depending on the purpose of your application, I would have tried to
get access to some of the raw data sampled at 2 kHz. If you can
make a good inverse filter, it would make sense to apply it as
one of the first steps in the processing chain, preferably before
decimation and storage.

Rune