What is the physical interpretation of sos in Zp2sos function?

5 次查看（过去 30 天）
Andi
Andi 2022-11-17
回答： Paul 2022-11-23
Hi everyone,
I am trying to filtering my signal using the following script.
N=length(t);
dt=(t(end)-t(1))/(N-1); %sampling interval (days)
fs=1/dt; fNyq=fs/2; %sampling frequency, Nyquist frequency (cycles/day)
f1=1; f2=4; %bandpass corner frequencies (cycles/day)
[z,p,k]=butter(4,[f1,f2]/fNyq,'bandpass'); %4th order Butterworth bandpass filter coefficients
[sos,g]=zp2sos(z,p,k); %sos representation of the filter
xde=filtfilt(sos,g,y); %apply zero phase filter to ser1
But, i am still not sure about the physical interpretation of few parameteres.
For example:
[z,p,k] > are poles, zeroes and gain are converted to sos and g. here what exactly sos shows and why we need to do this.
In my case, i get this as sos matrix, but i am not clear how to interpret this
1 2 1 1 -1.9836 0.9838
1 2 1 1 -1.9918 0.9921
1 -2 1 1 -1.9922 0.9922
1 -2 1 1 -1.9979 0.9979

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

回答（2 个）

Amey Waghmare
Amey Waghmare 2022-11-23
Hi,
As per my understanding, you have used ‘zp2sos’ to obtain Second order section (sos) representation of the filter, and want to understand the interpretation of the resulting sos matrix.
Please refer to the attached image:
The sos representation is returned as a L-by-6 matrix, where L is the number of sections, which depends upon the number of poles and zeros of the filter. Each rows contains the numerator and denominator coefficients 'bik' and 'aik' of the second order sections of the filter. The structure of the filter is shown by H(z) in the above image.
For more information, please refer to the following documentation on zp2sos: https://in.mathworks.com/help/signal/ref/zp2sos.html
Hope this resolves the issue.
Paul
Paul 2022-11-23
Ran in:
Hi Andi,
Addresing the second part of your question: " why we need to do this"
filtfilt will accept either transfer function (b,a) or sos (sos,g) form of the the filter, but not zpk form. So we have to choose one or the other. sos form has better numerical properties for computation.
Example filter parameters
n = 6;
Wn = [2.5e6 29e6]/500e6;
ftype = 'bandpass';
% Zero-Pole-Gain design
[z,p,k] = butter(n,Wn,ftype);
Note that all of the poles are inside the unit circle, as expected.
sort(abs(p).')
ans = 1×12
0.8662 0.8662 0.8985 0.8985 0.9613 0.9613 0.9823 0.9823 0.9893 0.9893 0.9965 0.9965
Convert to transfer function form
[b,a] = zp2tf(z,p,k);
Because of floating point issues, some of the poles are now outside the unit circle, which would be an unstable filter.
sort(abs(roots(a)).')
ans = 1×12
0.8655 0.8655 0.8988 0.8988 0.9336 0.9336 0.9614 0.9614 1.0019 1.0019 1.0358 1.0358
Convert to sos form
[sos,g] = zp2sos(z,p,k);
The poles of the sections are inside the unit circle, and they should be very, very close matches to p.
C = cellfun(@(x) roots(x(4:6).'),mat2cell(sos,ones(size(sos,1),1)),'UniformOutput',false);
sort(abs(vertcat(C{:})).')
ans = 1×12
0.8662 0.8662 0.8985 0.8985 0.9613 0.9613 0.9823 0.9823 0.9893 0.9893 0.9965 0.9965
Do filtfilt with both forms, observe very different results
rng(101)
x = rand(1000,1);
y1 = filtfilt(b,a,x);
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.181425e-17.
y2 = filtfilt(sos,g,x);
plot(y1)
plot(y2)
Note that filter, which is not part of the signal processing toolbox, does not accept the sos form.

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by