1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
|
//phasez Phase response of digital filter
//Calling Syntax
//[phi,w] = phasez(b,a,n)
//[phi,w] = phasez(sos,n)
///[phi,w] = phasez(b,a,n) returns
//the n-point unwrapped phase response vector, phi,
//in radians and the frequency vector, w, in radians/sample
//for the filter coefficients specified in b and a.
//The values of the frequency vector, w, range from
//0 to π. If n is omitted,
//the length of the phase response vector defaults to 512. For best
//results, set n to a value greater than the filter
//order.
//[phi,w] = phasez(sos,n) returns the unwrapped
//phase response for the second order sections matrix, sos. sos is
//a K-by-6 matrix, where the number of sections, K,
//must be greater than or equal to 2. If the number of sections is less
//than 2, phasez considers the input to be the
//numerator vector, b. Each row of sos corresponds
//to the coefficients of a second-order (biquad) filter. The ith
//row of the sos matrix corresponds to [bi(1)
//bi(2) bi(3) ai(1) ai(2) ai(3)].
////Author: Parthasarathi Panda
//parthasarathipanda314@gmail.com
function [phi, varargout]=phasez(varargin)
//cas variable is 2 if sos form is involved and 1 if direct rational form is given
//(sos,n) or (sos,w) or (sos,'whole')or (b,a) is the input
//cas variable is 2 if sos form is involved and 1 if direct rational form is given
//cas1 variable is 1 if f is to be given as output, 2 other wise
[nargout,nargin]=argn();
//do not forget to execute 'phaseInputParseAs_sos' and 'phaseInputParseAs_ab' before running
v=size(varargin(1));
if size(v)>2 then
error ('invalid input dimension');
end
[n,k]=size(varargin(1));
if type(varargin(1))~=1 then
error ('check the input type');
end
if (n==1 & k==6) then //not clear if sos or (a,b)
v=size(varargin(2));
if (nargin==1) //(sos) is the input
cas=2;
[sos,w,cas1,fs]=phaseInputParseAs_sos(varargin,nargin);
elseif (varargin(2)=='whole') //(sos,'whole')is the input
cas=2;
[sos,w,cas1,fs]=phaseInputParseAs_sos(varargin,nargin);
else //taking it as (a,b)
cas=1;
[a,b,w,cas1,fs]=phaseInputParseAs_ab(varargin,nargin);
end
elseif (n==1 | k==1) then
cas=1;
[a,b,w,cas1,fs]=phaseInputParseAs_ab(varargin,nargin);
elseif k==6 then //first variable is sos
cas=2;
[sos,w,cas1,fs]=phaseInputParseAs_sos(varargin,nargin);
end
//cas,cas1,fs,w,[(a,b),sos]
if cas==1 then
[m,n]=size(a);
N=[0:n-1];
M=N'*w;//computing matrix Mij=(i-1)*wj
ph_num=phasemag(a*exp(%i*M));//the operation computes phase of sum(ak*exp(i*w*k))
[m,n]=size(b);
N=[0:n-1];
M=N'*w;
ph_den=phasemag(b*exp(%i*M));//similar result for denominator
[m,n]=size(w);
phi=pmodulo(ph_num-ph_den,360);//takes the difference in phase modulo 360
else
N=[0,1,2];
M=N'*w;
ph_num=phasemag(sos(:,4:6)*exp(%i*M));
ph_den=phasemag(sos(:,1:3)*exp(%i*M));//the numerator phases for each second order componenet
phi_mat=ph_num-ph_den;
[m,n]=size(w);
phi=pmodulo(sum(phi_mat,1),360);//summing each of the componenet second order system phases
end
if cas1==1 then
varargout(2)=w*fs/(2*%pi);
if nargout>2 then
varargout(3)=struct('plot', 'both', 'fvflag', 0, 'yunits','degrees','xunit','Hz','fs',fs);
end
else
varargout(2)=w;
if nargout>2 then
varargout(3)=struct('plot', 'both', 'fvflag', 0, 'yunits','degrees','xunit','radian/sample','fs',[]);
end
end
endfunction
|
