Skip to content

Commit 4497f28

Browse files
somefunAgbasomefunAgba
committed
first commit
0 parents  commit 4497f28

15 files changed

+786
-0
lines changed

‎.gitignore

+42
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
# Created by .ignore support plugin (hsz.mobi)
2+
### MATLAB template
3+
# Windows default autosave extension
4+
*.asv
5+
6+
# OSX / *nix default autosave extension
7+
*.m~
8+
9+
# Compiled MEX binaries (all platforms)
10+
*.mex*
11+
12+
# Packaged app and toolbox files
13+
*.mlappinstall
14+
*.mltbx
15+
16+
# Generated helpsearch folders
17+
helpsearch*/
18+
19+
# Simulink code generation folders
20+
slprj/
21+
sccprj/
22+
23+
# Matlab code generation folders
24+
codegen/
25+
26+
# Simulink autosave extension
27+
*.autosave
28+
29+
# Simulink cache files
30+
*.slxc
31+
32+
# Octave session info
33+
octave-workspace
34+
35+
### Example user template template
36+
### Example user template
37+
38+
# IntelliJ project files
39+
.idea
40+
*.iml
41+
out
42+
gen

‎funcs/bpc_skern.m

+92
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
function bpck = bpc_skern(n,k,varargin)
2+
% BPC_KERN
3+
% full interface: bpc_kern(n,type,rho)
4+
% Generating Function for Binomial Polynomial Coefficients
5+
% of any real non-zero order n
6+
7+
8+
% Generates Pascal-Triangle's Binomial Coefficients
9+
% for integer order n
10+
11+
% type = 0
12+
% Coefficients of the normalized Binomial Transform
13+
% - useful in denominator polynomials of bilinear transform from s2z
14+
% type = 1
15+
% Coefficients of the normalized Inverse Binomial Transform
16+
% - useful in numerator polynomials of bilinear transform from s2z
17+
% Damping factor inclusion for Somefun Filters
18+
% - useful for coefficients of somefun filters.
19+
20+
assert(nargin > 0 && nargin < 5, "min of 1, max of 4 arguments")
21+
if nargin == 2
22+
type = 0; % 0 - default bin, 1 - bin. inverse
23+
rho = 1; % default damping factor
24+
elseif nargin == 3
25+
type = varargin{1};
26+
rho = 1; % default damping factor
27+
elseif nargin > 3
28+
type = varargin{1};
29+
rho = varargin{2}; % default damping factor
30+
end
31+
% input test
32+
assert(n > 0, "should be a positive natural number ");
33+
% we set the limit of n as (2^16) which is a big positive number
34+
% for the sake of safety and repeatability
35+
% with respect to computer numerical constraints
36+
assert(n < 65536, "must not exceed the system's maximum integer");
37+
38+
assert (k <= n && k >= 0, "k must be within limits");
39+
40+
% initialize a vector of ones with size (n+1)
41+
bpck = 1;
42+
43+
flag = 0; % is it even or odd
44+
% compute symmetry stop point
45+
if mod(n,2) == 0
46+
% if even
47+
ns = n/2;
48+
flag = 0; % even
49+
elseif mod(n,2) ~= 0
50+
% if odd
51+
ns = (n-1)/2;
52+
flag = 1; % odd
53+
end
54+
55+
if k >= 0 && k <= ns
56+
57+
% for k== 0
58+
if (k == 0)
59+
bpck = 1;
60+
elseif (k == 1)
61+
% for k == 1
62+
bpck = n*rho;
63+
elseif (k > 1)
64+
% for k > 1
65+
num = n;
66+
for i = 1:(k-1)
67+
num = num * (n - i);
68+
end
69+
den = 1;
70+
for i = 2:k
71+
den = den * i;
72+
end
73+
val = num / den;
74+
bpck = val*rho;
75+
end
76+
77+
if type == 1
78+
bpck = (-1)^k * bpck;
79+
end
80+
81+
end
82+
83+
% symmetry for k > n-k
84+
if k >= (ns+1) && k <=n % k = ns+1: n
85+
if type == 0 || flag == 0
86+
bpck = bpc_skern(n,n-k,type,rho);
87+
elseif (type == 1) && (flag == 1)
88+
bpck = (-1) * bpc_skern(n,n-k,type,rho);
89+
end
90+
end
91+
92+
end

‎funcs/dbpoly_kern.m

+100
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
function dbpoly = dbpoly_kern(n,varargin)
2+
% BPC_KERN
3+
% Generating Function for Coefficients of a
4+
% Damped Binomial Polynomial in vector form
5+
% of any real non-zero order n
6+
7+
% Interface: dbpoly_kern(n,type,mode)
8+
9+
10+
% Generates Damped Pascal-Triangle's Binomial Coefficients
11+
% for integer order n
12+
13+
% type = 0
14+
% Coefficients of the normalized Binomial Transform
15+
% - useful in denominator polynomials of bilinear transform from s2z
16+
% type = 1
17+
% Coefficients of the normalized Inverse Binomial Transform
18+
% - useful in numerator polynomials of bilinear transform from s2z
19+
% Damping factor inclusion for Somefun Filters
20+
% - useful for coefficients of somefun filters.
21+
22+
% mode = 1, rho = damped.
23+
% mode = 0, rho = 1;
24+
25+
assert(nargin > 0 && nargin < 4, "min of 1, max of 3 arguments")
26+
if nargin == 1
27+
type = 0; % 0 - default bin, 1 - bin. inverse
28+
mode = 1; % default damping mode
29+
elseif nargin == 2
30+
type = varargin{1};
31+
mode = 1; % default damping mode
32+
else
33+
type = varargin{1};
34+
mode = varargin{2}; % default damping factor
35+
end
36+
% input test
37+
assert(n > 0, "should be a positive natural number ");
38+
% we set the limit of n as (2^16) which is a big positive number
39+
% for the sake of safety and repeatability
40+
% with respect to computer numerical constraints
41+
assert(n < 65536, "must not exceed the system's maximum integer");
42+
43+
rho = 1;
44+
if mode == 1
45+
% Damping value
46+
rho = rhoval(n);
47+
end
48+
% initialize a vector of ones with size (n+1)
49+
dbpoly = ones(1,n+1);
50+
51+
flag = 0; % is it even or odd
52+
% compute symmetry stop point
53+
if mod(n,2) == 0
54+
% if even
55+
ns = n/2;
56+
flag = 0; % even
57+
elseif mod(n,2) ~= 0
58+
% if odd
59+
ns = (n-1)/2;
60+
flag = 1; % odd
61+
end
62+
63+
for k = 0:ns
64+
65+
% for k== 0
66+
if (k == 0)
67+
dbpoly(k+1) = 1;
68+
elseif (k == 1)
69+
% for k == 1
70+
dbpoly(k+1) = n*rho;
71+
elseif (k > 1)
72+
% for k > 1
73+
num = n;
74+
for i = 1:(k-1)
75+
num = num * (n - i);
76+
end
77+
den = 1;
78+
for i = 2:k
79+
den = den * i;
80+
end
81+
val = num / den;
82+
dbpoly(k+1) = val*rho;
83+
end
84+
85+
if type == 1
86+
dbpoly(k+1) = (-1)^k * dbpoly(k+1);
87+
end
88+
89+
end
90+
91+
% symmetry for k > n-k
92+
for k = (ns+1):n % k = ns+1: n
93+
if type == 0 || flag == 0
94+
dbpoly(k+1) = dbpoly(n-k+1);
95+
elseif (type == 1) && (flag == 1)
96+
dbpoly(k+1) = (-1) * dbpoly(n-k+1);
97+
end
98+
end
99+
100+
end

‎funcs/grpdel.m

+28
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
function tau_g = grpdel(freqw,wn,n)
2+
magsq_udbf = magresp_udbf(freqw,wn,n);
3+
fp_term = magsq_udbf;
4+
5+
fl_term = n + (freqw./wn).^((2*n) - 2);
6+
7+
mid_term = 0;
8+
rho = rhoval(n);
9+
10+
for id=(n-2):-1:1
11+
t = n-id-1; lambda_t = 0;
12+
if ( (id >= (n/2) ) && isintegereven(n) ) ...
13+
|| ( (id >= ((n-1)/2) ) && ~isintegereven(n) )
14+
r_bar = n-id;
15+
else
16+
r_bar=id+1;
17+
end
18+
for r=1:1:r_bar
19+
j = t+1-r; k = t+r;
20+
lambda_t = lambda_t + ((-1)^(r-1))*(bpc_skern(n,j,0,rho))*(bpc_skern(n,k,0,rho));
21+
end
22+
mid_term = mid_term + lambda_t.*(freqw./wn).^(2*id);
23+
end
24+
25+
sump_term = fl_term + mid_term;
26+
tau_g = fp_term.*sump_term;
27+
28+
end

‎funcs/magresp_udbf.m

+60
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
function [magsq_udbf, mag_udbf, magsq_but, mag_but, magsq_bin, mag_bin] = magresp_udbf(freqw,wn,n)
2+
3+
fl_term = 1 + (freqw./wn).^(2*n);
4+
magsq_but = 1./fl_term;
5+
mag_but = 1./sqrt(fl_term);
6+
7+
% uniformly-damped BINOMIAL
8+
9+
mid_term = 0;
10+
o_term=0; % origin term;
11+
rho = rhoval(n);
12+
13+
for id=(n-1):-1:1
14+
t = n-id; alpha_t = 0;
15+
if ( (id >= (n/2) ) && isintegereven(n) ) ...
16+
|| ( (id >= ((n+1)/2) ) && ~isintegereven(n) )
17+
r_bar = n-id;
18+
else
19+
r_bar=id;
20+
end
21+
for r=1:1:r_bar
22+
j = t-r; k = t+r;
23+
alpha_t = alpha_t + ((-1)^(r))*(bpc_skern(n,j,0,rho))*(bpc_skern(n,k,0,rho));
24+
end
25+
alpha_t = (bpc_skern(n,t,0,rho))^2 + 2*alpha_t;
26+
mid_term = mid_term + alpha_t.*(freqw./wn).^(2*id);
27+
% o_term = o_term + alpha_t;
28+
29+
end
30+
% disp(o_term)
31+
den = fl_term + mid_term;
32+
magsq = 1./den;
33+
magsq_udbf = magsq;
34+
mag_udbf = 1./sqrt(den);
35+
36+
% STANDARD BINOMIAL
37+
mid_term = 0;
38+
rho = 1;
39+
for id=(n-1):-1:1
40+
t = n-id; alpha_t = 0;
41+
if ( (id >= (n/2) ) && isintegereven(n) ) ...
42+
|| ( (id >= ((n+1)/2) ) && ~isintegereven(n) )
43+
r_bar = n-id;
44+
else
45+
r_bar=id;
46+
end
47+
for r=1:1:r_bar
48+
j = t-r; k = t+r;
49+
alpha_t = alpha_t + ((-1)^(r))*(bpc_skern(n,j,0,rho))*(bpc_skern(n,k,0,rho));
50+
end
51+
alpha_t = (bpc_skern(n,t,0,rho))^2 + 2*alpha_t;
52+
mid_term = mid_term + alpha_t.*(freqw./wn).^(2*id);
53+
% o_term = o_term + alpha_t;
54+
end
55+
% disp(o_term);
56+
den = fl_term + mid_term;
57+
magsq_bin = 1./den;
58+
mag_bin = 1./sqrt(den);
59+
60+
end

‎funcs/phasedel.m

+37
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
function [phase, tau_p, inum, rden] = phasedel(freqw,wn,n)
2+
3+
tau_p = wn./freqw;
4+
rho = rhoval(n);
5+
inum = 0; rden = 0;
6+
7+
for id = 0:1:n
8+
9+
if (~isintegereven(id))
10+
x = (id-1)/2;
11+
w = ((-1).^(x)).*(freqw./wn).^(id);
12+
c = (bpc_skern(n,id,0,rho));
13+
inum = inum + (w.*c);
14+
end
15+
if (isintegereven(id))
16+
x = (id)/2;
17+
w = ((-1).^(x)).*(freqw./wn).^id;
18+
c = (bpc_skern(n,id,0,rho));
19+
rden = rden + (w.*c);
20+
end
21+
22+
end
23+
24+
% x = inum./rden;
25+
phase = -(atan2(inum,rden));
26+
27+
% correct the radian phase angles in a phase angle vector by adding
28+
% multiples of ±2pi when absolute jumps between consecutive elements of the vector
29+
% are greater than or equal to the default jump tolerance of pi radians.
30+
phase = unwrap(phase);
31+
32+
% convert to deg.
33+
% phase = phase*180/pi;
34+
35+
tau_p = -tau_p .* phase;
36+
37+
end

‎funcs/rhoval.m

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
function rho = rhoval(n)
2+
% RHOVAL
3+
% full interface: rhoval(n)
4+
5+
% Damping factor coefficient generation for Damped Binomial Filters
6+
% - useful for synthesis of damped binomial filters.
7+
8+
% input test
9+
assert(n > 0, "n should be a positive natural number ");
10+
% we set the limit of n as (2^16) which is a big positive number
11+
% for the sake of safety and repeatability
12+
% with respect to computer numerical constraints
13+
assert(n < 65536, "n must not exceed the system's maximum integer");
14+
15+
rho = sqrt( (n*(n-1))-(n-2) ) / n;
16+
% if n == 1
17+
% sfunc = 1;
18+
% elseif n == 2
19+
% sfunc = sqrt(2)/2;
20+
% elseif n == 3
21+
% sfunc = sqrt(5)/3;
22+
% elseif n > 3
23+
% sfunc = sqrt((n*(n-1))-(n-2))/n; % n^2 - n - 2
24+
% end
25+
26+
end

0 commit comments

Comments
 (0)