s-plane to z-plane transformation

23 次查看（过去 30 天）
Markus
Markus 2011-4-28
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

请先登录，再回答此问题。

采纳的回答

Teja Muppirala
Teja Muppirala 2011-4-28
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'})

更多回答（2 个）

Markus
Markus 2011-4-28
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 :-)
Markus
Markus 2011-4-28
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'})

类别

在 Help Center 和 File Exchange 中查找有关 Analog Filters 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by