s-plane to z-plane transformation

Hi,

I am working in a project where I need to filter a sound signal. The filter is called G-filter and it is specified in the ISO 7196 standard (Frequency-weighing characteristics for infrasound measurements).

The filter is specified in s-plane by it's poles and zeros. My signal is digital (time discrete) and thus I need to estimate a digital filter. I have tried yule-walker and bilinear transform but neither gives an acceptable estimate. Is there a better way?

clear all
close all

%%Define G-filter in s-plane according to ISO 7196:1995E
%poles
p=[ -0.707+j* 0.707
    -0.707-j* 0.707
   -19.27+j* 5.16
   -19.27-j* 5.16
   -14.11+j*14.11
   -14.11-j*14.11
    -5.16+j*19.27
    -5.16-j*19.27];
%zeros
z=[0 0 0 0]';
%scalar gain
k=10^(116/20);
%get analog filter coefficients
[bG,aG] = zp2tf(z,p,k);
%get filter response
[hG,fG] = freqs(bG,aG,[0 logspace(log10(0.2),log10(200),1000)]);

%%Try to get digital coefficients from yule-walker method
fs=400;
[b,a] = yulewalk(8,fG/max(fG),abs(hG));
%get digital filter response
[h,f] = freqz(b,a,100*fs,fs);

%%Try to get digital coefficients from bilinear transform
[bz,az]=bilinear(bG,aG,fs);
%get digital filter response
[Hz, fz]= freqz(bz,az,100*fs,fs);

%%plot
figure
semilogx(fG,20*log10(abs(hG)),'--'...
    ,f,20*log10(abs(h))...
    ,fz,20*log10(abs(Hz)))
xlim([0.2 fs])
xlabel('frequency (Hz)'), ylabel('response (dB)')
legend('G in s-plane','G in z-plane, yule-walker','G in z-plane, bilinear');
set(gca,'XTick',[0.2 0.5 1 2 5 10 20 50 100 200])
set(gca,'XTickLabel',{'0,2','0,5','1','2','5','10','20','50','100','200'})

The blue curve is the desired filter response. According to the standard, the response is most important in the 1-20Hz range where errors of 1dB are allowed. Below 1Hz and above 20Hz the error is allowed to be -infinity to +1dB. Clearly the estimated filters are not good enough. Any advice is appreciated! Thanks, Markus

Answer 1

[Teja Muppirala]

Have you tried c2d?

clear all
close all
%%Define G-filter in s-plane according to ISO 7196:1995E
%poles
p=[ -0.707+j* 0.707
    -0.707-j* 0.707
   -19.27+j* 5.16
   -19.27-j* 5.16
   -14.11+j*14.11
   -14.11-j*14.11
    -5.16+j*19.27
    -5.16-j*19.27];
%zeros
z=[0 0 0 0]';
%scalar gain
k=10^(116/20);

G_0 = zpk(z,p,k)

fs = 400;
G_1 = c2d(G_0,1/fs,'tustin');

bode(G_0);
hold all;
bode(G_1);
legend({'continuous' 'discrete filter'})

% Extract the filter coefficients
G_1_tf = tf(G_1);
B = G_1_tf.num{1}
A = G_1_tf.den{1}

% Compare the time series output using a random input signal
figure;
r = randn(1,10000);
t = 1/fs * (0:9999);
y1 = lsim(G_0,r,t);

y2 = filter(B,A,r);

plot(t,y1,t,y2,'r:');
legend({'continuous' 'discrete filter'})

Answer 2

[Markus]

Hi Teja, thanks for your answer.

Unfortunately I don't have the system toolboxes so I can't use c2d( ) and tf( ).

However, I found that if I divide the frequency by 2pi in the bilinear function I get a much better output :-)

Answer 3

[Markus]

The only difference is

[bz,az]=bilinear(bG,aG,fs/(2*pi));

instead of

[bz,az]=bilinear(bG,aG,fs);

The code:

clear all
close all

%%G-filter in s-plane according to ISO 7196:1995E
%poles
p=[ -0.707+j* 0.707
    -0.707-j* 0.707
   -19.27+j* 5.16
   -19.27-j* 5.16
   -14.11+j*14.11
   -14.11-j*14.11
    -5.16+j*19.27
    -5.16-j*19.27];
%zeros
z=[0 0 0 0]';
%scalar gain
k=10^(116/20);
%get analog filter coefficients
[bG,aG] = zp2tf(z,p,k);
%get filter response
[hG,fG] = freqs(bG,aG,[0 logspace(log10(0.2),log10(200),1000)]);

%%Try to get digital coefficients from yule-walker method
fs=200;
[b,a] = yulewalk(8,fG/max(fG),abs(hG));
%get digital filter response
[h,f] = freqz(b,a,100*fs,fs);

%%Try to get digital coefficients from bilinear transform
[bz,az]=bilinear(bG,aG,fs/(2*pi));
%get digital filter response
[Hz, fz]= freqz(bz,az,100*fs,fs);

%%plot
figure
semilogx(fG,20*log10(abs(hG)),'--'...
    ,f,20*log10(abs(h))...
    ,fz,20*log10(abs(Hz)))
xlim([0.2 fs])
xlabel('frequency (Hz)'), ylabel('response (dB)')
legend('G in s-plane','G in z-plane, yule-walker','G in z-plane, bilinear');
set(gca,'XTick',[0.2 0.5 1 2 5 10 20 50 100 200])
set(gca,'XTickLabel',{'0,2','0,5','1','2','5','10','20','50','100','200'})