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 allclose 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
figuresemilogx(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
Best Answer