diff options
34 files changed, 3813 insertions, 1526 deletions
diff --git a/macros/besselap.sci b/macros/besselap.sci index 2ed0047..1e7d952 100644 --- a/macros/besselap.sci +++ b/macros/besselap.sci @@ -1,17 +1,14 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Original contribution: FOSSEE, IIT Bombay // Original Source : https://octave.sourceforge.io/signal/ -// Modifieded by: Sonu Sharma, RGIT Mumbai +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - - function [zero, pole, gain]=besselap(n) //Bessel analog filter prototype. @@ -44,38 +41,44 @@ function [zero, pole, gain]=besselap(n) // zero = // // [] - + // Dependencies + // prepad funcprot(0); [nargout, nargin] = argn() ; - if (nargin>1 | nargin<1) error("besselap : wrong number of input argument") end - // interpret the input parameters if (~(length(n)==1 & n == round(n) & n > 0)) error ("besselap: filter order n must be a positive integer"); end - p0=1; p1=[1 1]; for nn=2:n px=(2*nn-1)*p1; py=[p0 0 0]; - px=prepad(px,max(length(px),length(py)),0); + px=prepad(px,max(length(px),length(py))); py=prepad(py,length(px)); p0=p1; p1=px+py; end // p1 now contains the reverse bessel polynomial for n - // scale it by replacing s->s/w0 so that the gain becomes 1 p1=p1.*p1(length(p1)).^((length(p1)-1:-1:0)/(length(p1)-1)); - zero=[]; pole=roots(p1); gain=1; - endfunction + +/* +Note : The function is tested with Octave's outputs as a reference. +# all passed +[zero, pole, gain] = besselap (1) +[zero, pole, gain] = besselap (2) +[zero, pole, gain] = besselap (7) +[zero, pole, gain] = besselap (13) +[zero, pole, gain] = besselap (43) + +*/
\ No newline at end of file diff --git a/macros/besself.sci b/macros/besself.sci index c52f3ac..a6eef80 100644 --- a/macros/besself.sci +++ b/macros/besself.sci @@ -1,15 +1,14 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt // Original Source : https://octave.sourceforge.io/signal/ -// Modifieded by:Sonu Sharma, RGIT Mumbai +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - function [a, b, c, d] = besself (n, w, varargin) //Bessel filter design. @@ -20,22 +19,17 @@ function [a, b, c, d] = besself (n, w, varargin) //[b, a] = besself (n, [Wl, Wh], "stop") //[z, p, g] = besself (…) //[…] = besself (…, "z") - - //Parameters //n: positive integer value (order of filter) //Wc: positive real value, // 1).Analog 3dB cutoff frequency/frequencies for analog filter, in the range [0, Inf] {rad/sec} // 2).Normalised digital 3dB cutoff frequency/frequencies for digital filter, in the range [0, 1] {dimensionless} - //Description //This function generates a Bessel filter. The default is a Laplace space (s) or analog filter. //If second argument is scalar the third parameter takes in high or low, the default value being low. The cutoff is Wc rad/sec. //If second argument is vector of length 2 ie [Wl Wh] then third parameter may be pass or stop default is pass for bandpass and band reject filter respectively //[z,p,g] = besself(...) returns filter as zero-pole-gain rather than coefficients of the numerator and denominator polynomials. //[...] = besself(...,’z’) returns a discrete space (Z) filter. Wc must be less than 1 {dimensionless}. - - //Examples //[b, a]=besself(2,.3,"high","z") //Output : @@ -46,14 +40,14 @@ function [a, b, c, d] = besself (n, w, varargin) // // 0.4668229 - 0.9336457 0.4668229 // - + // Dependencies + // besselap bilinear sftrans zp2tf funcprot(0); [nargout nargin] = argn(); if (nargin > 4 | nargin < 2 | nargout > 4 | nargout < 2) error("besself: invalid number of inputs") end - // interpret the input parameters if (~ (isscalar (n) & (n == fix (n)) & (n > 0))) error ("besself: filter order N must be a positive integer"); @@ -79,14 +73,13 @@ function [a, b, c, d] = besself (n, w, varargin) error ("besself: expected [high|stop] or [s|z]"); end end - // FIXME: Band-pass and stop-band currently non-functional, remove // this check once low-pass to band-pass transform is implemented. if (~ isscalar (w)) error ("besself: band-pass and stop-band filters not yet implemented"); end - [rows_w colums_w] = size(w); + [rows_w columns_w] = size(w); if (~ ((length (w) <= 2) & (rows_w == 1 | columns_w == 1))) error ("besself: frequency must be given as WC or [WL, WH]"); @@ -115,7 +108,6 @@ function [a, b, c, d] = besself (n, w, varargin) if (digital) [zero, pole, gain] = bilinear (zero, pole, gain, T); end - // convert to the correct output form if (nargout == 2) // a = real (gain * poly (zero)); @@ -127,7 +119,28 @@ function [a, b, c, d] = besself (n, w, varargin) c = gain; else // output ss results - // [a, b, c, d] = zp2ss (zero, pole, gain); + // FIXME : test zp2ss + //[a, b, c, d] = zp2ss (zero, pole, gain); error("besself: yet not implemented in state-space form OR invalid number of o/p arguments") end endfunction + +/* +Note : This function is tested with Octave's outputs as a reference. +# Test input validation +[a, b] = besself () // error passed +[a, b] = besself (1) // error passed +[a, b] = besself (1, 2, 3, 4, 5) // error passed +[a, b] = besself (.5, .2) // error passed +[a, b] = besself (3, .2, "invalid") // error passed + //[b, a] = besself(n, Wc) +[b a] = besself(2,1000) // passed +//[b, a] = besself (n, Wc, "high") +[b a] = besself(7,0.3,"high") // passed +//[b, a] = besself (n, [Wl, Wh]) +[b a] = besself(4,[1000 5000]) // error: besself: band-pass and stop-band filters not yet implemented +//[b, a] = besself (n, [Wl, Wh], "stop") +[b a] = besself(4,[1000 5000],"stop") // besself: band-pass and stop-band filters not yet implemented +//[z, p, g] = besself (…) +[z,p,g] = besself(9,0.8,"high") // passed +*/
\ No newline at end of file diff --git a/macros/bilinear.sci b/macros/bilinear.sci index d1ee3f7..1aa408a 100644 --- a/macros/bilinear.sci +++ b/macros/bilinear.sci @@ -1,15 +1,14 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt // Original Source : https://octave.sourceforge.io/signal/ -// Modifieded by:Sonu Sharma, RGIT Mumbai +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T) //Transforms a s-plane filter (Analog) into a z-plane filter (Digital) using Bilinear transformation @@ -48,13 +47,15 @@ function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T) // b = // // 0. - 0.1666667 - 0.3333333 2.5 - + // Dependencies + // tf2zp postpad zp2tf prepad funcprot(0); [nargout nargin] = argn(); ieee(2); if nargin==3 T = Sg; + // FIXME : tf2zp is not tested yet [Sz, Sp, Sg] = tf2zp(Sz, Sp); elseif nargin~=4 error("bilinear: invalid number of inputs") @@ -77,9 +78,6 @@ function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T) if Zg == 0 & nargout == 3 then error("bilinear: invalid value of gain due to zero(s) at infinity avoid z-p-g form and use tf form ") end - - - Zp = (2+Sp*T)./(2-Sp*T); SZp = size(Zp); if isempty(Sz) @@ -88,8 +86,6 @@ function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T) Zz = [(2+Sz*T)./(2-Sz*T)]; Zz = postpad(Zz, p, -1); end - - if nargout==2 // zero at infinity Zz1 = []; @@ -115,3 +111,17 @@ function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T) end ieee(0); endfunction +/* +// FIXME- not working with three argument + +Note : This function is tested with Octave's outputs as a reference. + +[Zb,Za] = bilinear([1 0],[1 1],1,0.5) // passed +[Zb, Za] = bilinear([], [], 1, 1) // error PASSED +[Zb, Za] = bilinear([0], [], 1, 0.5) // error PASSED +[Zb, Za] = bilinear([], [0], 1, 0.5) // PASSED + +[Zb, Za] = bilinear([2], [1], 1, 0.5) // PASSED +[Zb, Za] = bilinear([1; -1], [0.5; -0.5], 2, 1) //PASSED +[Zb, Za] = bilinear([0], [1], 1, 1) //PASSED +*/
\ No newline at end of file diff --git a/macros/bitrevorder.sci b/macros/bitrevorder.sci index cf27cd8..7ffe2ec 100644 --- a/macros/bitrevorder.sci +++ b/macros/bitrevorder.sci @@ -1,54 +1,48 @@ -// Returns input data in bit-reversed order - -// Calling Sequence -//[y,i] = bitrevorder(x) -//y = bitrevorder(x) - -// Parameters -//x: Vector of real or complex values -//y: input vector in bit reverse order -//i: indices - -// Description -//This function returns the input data after reversing the bits of the indices and reordering the elements of the input array. - -// Examples -//x = [%i,1,3,6*%i] ; -//[y i]=bitrevorder(x) -//Output : -// i = -// -// 1. 3. 2. 4. -// y = -// -// i 3. 1. 6.i - - -//************************************************************************************* -//-------------------version1 (using callOctave / errored)----------------------------- -//************************************************************************************* - -//function [y,i]=bitrevorder(x) -//funcprot(0); -//[lhs,rhs]=argn(0); -//if (rhs<1) then -// error ("Wrong number of input arguments.") -//end -//[y,i]=callOctave("bitrevorder",x) -// -//endfunction - -//************************************************************************************* -//-----------------------------version2 (pure scilab code)----------------------------- -//************************************************************************************* +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in function [y, i] = bitrevorder (x) + // Returns input data in bit-reversed order + // + // Calling Sequence + //[y,i] = bitrevorder(x) + //y = bitrevorder(x) + + // Parameters + //x: Vector of real or complex values + //y: input vector in bit reverse order + //i: indices + + // Description + //This function returns the input data after reversing the bits of the indices and reordering the elements of the input array. + + // Examples + //x = [%i,1,3,6*%i] ; + //[y i]=bitrevorder(x) + //Output : + // i = + // + // 1. 3. 2. 4. + // y = + // + // i 3. 1. 6.i + // Dependencies + // digitrevorder funcprot(0); [nargout, nargin] = argn() ; if (nargin ~= 1) - print_usage (); - elseif (~ isvector (x)) + error("bitrevorder: Usage : [ y , i ] = bitrevorder(x) "); + elseif ( ~isvector(x) & ~isscalar(x) ) error ("bitrevorder: X must be a vector"); elseif (fix (log2 (length (x))) ~= log2 (length (x))) error ("bitrevorder: X must have length equal to an integer power of 2"); @@ -57,3 +51,22 @@ function [y, i] = bitrevorder (x) [y, i] = digitrevorder (x, 2); endfunction +/* +tests +assert_checkequal (bitrevorder (0), 0); //passed +assert_checkequal (bitrevorder (0:1), 0:1); //passed +assert_checkequal (bitrevorder ([0:1]'), [0:1]'); //passed +assert_checkequal (bitrevorder (0:7), [0 4 2 6 1 5 3 7]); //passed +assert_checkequal (bitrevorder ([0:7]'), [0 4 2 6 1 5 3 7]'); // passed +assert_checkequal (bitrevorder ([0:7]*%i), [0 4 2 6 1 5 3 7]*%i); // passed +assert_checkequal (bitrevorder ([0:7]'*%i), [0 4 2 6 1 5 3 7]'*%i); // passed +assert_checkequal (bitrevorder (0:15), [0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15]); //passed + + Test input validation + // error testing +assert_checkerror (bitrevorder ()); +assert_checkerror bitrevorder (1, 2); +assert_checkerror bitrevorder ([]); +assert_checkerror bitrevorder (0:2); + +*/
\ No newline at end of file diff --git a/macros/cell2sos.sci b/macros/cell2sos.sci index 39402b3..7cd657d 100644 --- a/macros/cell2sos.sci +++ b/macros/cell2sos.sci @@ -1,95 +1,110 @@ -function [s,g] = cell2sos(c) -//Converts a cell array to a second order section matrix -//Calling Sequences -//s=cell2sos(c) -//[s,g]=cell2sos(c) -//Parameters -//c -//A cell array -//g -//The scalar gain -//Description -//s=cell2sos(c) converts a a cell array c = { {B1},{A1}, {B2},{A2}, ... {BL},{AL}} -//to an L-by-6 second-order-section matrix s given by: -// s = [B1 A1 -// B2 A2 -// ... -// BL AL] -//numerator vector Bi and denominator vector Ai contains the coefficients of a -//linear or quadratic polynomial. If the polynomial is linear, the coefficients -//zero-padded on the right. -//[s,g]=cell2sos(c) estimates the gain from the leading term of the cell array -//c={ {[g1,g2]},{B1},{A1}, {B2},{A2}, ... {BL},{AL}} to give g=g1/g2 as the gain -//Example -//c=cell(1,5); -// -//c(1,1).entries=[2, 1]; -// -//c(1,2).entries=rand(1,3); -// -//c(1,3).entries=rand(1,3); -// -//c(1,4).entries=rand(1,3); -// -//c(1,5).entries=rand(1,3); -// -// c = -// column 1 to 3 -// -//![2,1] [0.2113249,0.7560439,0.0002211] [0.3303271,0.6653811,0.6283918] ! -// -// column 4 to 5 -// -//![0.8497452,0.6857310,0.8782165] [0.0683740,0.5608486,0.6623569] ! -//[s,g]=cell2sos(c); -//s = -// -// column 1 to 5 -// -// 0.2113249 0.7560439 0.0002211 0.3303271 0.6653811 -// 0.8497452 0.6857310 0.8782165 0.0683740 0.5608486 -// -// column 6 -// -// 0.6283918 -// 0.6623569 -// -//g = -// -// 2. -//Author -//Ankur Mallick - if(argn(2)~=1) then - error("Wrong number of input arguments"); - end +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Author: Abinash Singh Under FOSSEE Internship +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +/* +Calling Sequence : + sos = cell2sos(cll) + [sos,g] = cell2sos(cll) +Description + sos = cell2sos(cll) generates a matrix sos containing the coefficients of the filter system described by the second-order section cell array cll. + [sos,g] = cell2sos(cll) also returns the scale gain g. + + Second-order section cell-array representation, specified as a cell array. + +Input Argument: + For a filter system with L sections, specify 'cll' using this structure: + + * Cell array with L elements — For unity-gain filter systems. Each element of the cell + array corresponds to a second-order section. The kth cell array element of 'cll' + + cll{k} = {[b_0k b_1k b_2k] [1 a_1k a_2k]} + + contains the coefficients from the kth second-order-section of the filter system H(z): + + H(z) = product(k=1 to L) H_k(z) + = product(k=1 to L) (b_0k + b_1k*z^(-1) + b_2k*z^(-2))/(1 + a_1k*z^(-1) + a_2k*z^(-2)) + + * Cell array with L+1 elements — If the gain of the filter system is different from 1. + The first element of 'cll' contains the system gains at the numerator (g_n) and at + the denominator (g_d). Then, the function appends each element of the cell array for + the corresponding second-order section. + + The first and the k+1th cell array element of 'cll' - L=prod(size(c)); - for i=1:L - if(type(c(i))~=17) - error('Cell contents must themselves be cell objects'); + cll{1} = {g_n g_d} + cll{k+1} = {[b_0k b_1k b_2k] [1 a_1k a_2k]} + + contain the system gain and the coefficients from the kth second-order section of + the filter system H(z), respectively, such that: + + H(z) = (g_n/g_d) * product(k=1 to L) H_k(z) + = (g_n/g_d) * product(k=1 to L) (b_0k + b_1k*z^(-1) + b_2k*z^(-2))/(1 + a_1k*z^(-1) + a_2k*z^(-2)) + +Output Argument: + Second-order section representation, returned as an L-by-6 matrix, where L is the + number of second-order sections. The matrix + + sos = [b_01 b_11 b_21 1 a_11 a_21] + [b_02 b_12 b_22 1 a_12 a_22] + [ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ] + [b_0L b_1L b_2L 1 a_1L a_2L] + + represents the second-order sections of H(z): + + H(z) = g * product(k=1 to L) H_k(z) + = g * product(k=1 to L) (b_0k + b_1k*z^(-1) + b_2k*z^(-2))/(1 + a_1k*z^(-1) + a_2k*z^(-2)) +*/ + +function [s,g] = cell2sos(c) + if(argn(2)~=1) then + error("cell2sos: Wrong number of input arguments"); end - end - if (argn(1)==2) - d=c(1).entries; - if(length(d)==2) - g1=d(1); - g2=d(2); - g=g1/g2; - c=c(2:L); - else - g=1; + L=prod(size(c)); + k=1 + for i=1:L + if(type(c(i))~=17) + error('cell2sos: Cell contents must themselves be cell objects'); + end + end + gain_p = 0 + if length(cell2mat(c{1}))== 2 then + gain_vec=cell2mat(c{1}) + k=gain_vec(1)/gain_vec(2) + L=L-1 + gain_p=1 + end + s = zeros(L,6); + for i=1:L + if gain_p + s(i,:)=cell2mat(cll{i+1}) + else + s(i,:)=cell2mat(cll{i}) + end end - end - L=prod(size(c)); - s=zeros(L/2,6); - for i=1:2:L-1 - j=ceil(i/2) - b=c(i).entries; - a=c(i+1).entries; - b=b(:).'; - a=a(:).'; - b=[b,zeros(1,3-length(b))]; - a=[a,zeros(1,3-length(b))]; - s(j,:)=[b,a]; - end + if nargout < 2 then + s(1,1:3)= k * s(1,1:3) + else + g=k; + end endfunction + +/* +cll = {{[3 6 7] [1 1 2]} + {[1 4 5] [1 9 3]} + {[2 7 1] [1 7 8]}}; +sos = cell2sos(cll) // passed + +cll = {{1 2} {[3 6 7] [1 1 2]} + {[1 4 5] [1 9 3]} + {[2 7 1] [1 7 8]}}; +sos = cell2sos(cll) // passed + +*/
\ No newline at end of file diff --git a/macros/cummax.sci b/macros/cummax.sci index e9a68e5..465b515 100644 --- a/macros/cummax.sci +++ b/macros/cummax.sci @@ -1,4 +1,15 @@ -function M = cummax(varargin) +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 13 March 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +function [M , iM] = cummax(varargin) // Cumulative maximum // // Calling Sequence @@ -11,7 +22,8 @@ function M = cummax(varargin) // The operation is performed along the dimension specified by dim // M = cummax(_,direction) // direction specifies as the direction of operation - // + // [M , iM] = cummax(..) + // If called with two output arguments the index of the maximum value is also returned. // Parameters // A - real|complex numbers - vector|matrix // Input Array @@ -32,15 +44,6 @@ function M = cummax(varargin) // M = cummax(v) // // Expected output: [8 9 9 10 10 10 10 10 10 10] - // - // Authors - // Ayush Baid - // - // See Also - // cummax | cumprod | cumsum | max | max - - - [numOutArgs,numInArgs] = argn(0); // ** Checking number of arguments @@ -50,8 +53,8 @@ function M = cummax(varargin) error(77,msg); end - if numOutArgs~=1 then - msg = "cummax: Wrong number of output argument; 1 expected"; + if numOutArgs > 2 then + msg = "cummax: Wrong number of output argument; 1 or 2 expected"; error(78,msg); end @@ -103,11 +106,11 @@ function M = cummax(varargin) end // extracting direction - if strcmpi(directionArg,"reverse")==0 then + if strcmp(directionArg,"reverse")==0 then isForward = %f; - elseif strcmpi(directionArg,"forward")==0 then + elseif strcmp(directionArg,"forward")==0 then isForward = %t; - elseif strcmpi(directionArg,"")~=0 then + elseif strcmp(directionArg,"")~=0 then msg = "cummax: Wrong value for argument #3 (direction)"; error(53,msg); end @@ -132,7 +135,18 @@ function M = cummax(varargin) end M = matrix(M_,sizeA); - + + if numOutArgs == 2 then + // calculating the index + // for vectors + iM = zeros(sizeA(1),sizeA(2)); + for i=1:sizeA(1) + for j=1:sizeA(2) + index = find (M(i,j) == A(i,:) ) + iM(i,j) = index(1) + end + end + end endfunction @@ -191,6 +205,16 @@ function out = cummaxVec(inp,isForward) end end end - - endfunction +/* +# tests +[w, iw] = cummax ([1 3 2 6 4 5]); // passed +x = [1 2 3; 4 1 2; 3 5 1]; +w = cummax(x); //passsed + +x = [1 2 3; 4 1 2; 3 5 1]; +w = cummax(x, 2); // passed + +x = [1 2 3; 4 1 2; 3 5 1]; +[w,iw] = cummax(x, 2); // passed +*/
\ No newline at end of file diff --git a/macros/cummin.sci b/macros/cummin.sci index fccbd52..7a4db33 100644 --- a/macros/cummin.sci +++ b/macros/cummin.sci @@ -1,4 +1,15 @@ -function M = cummin(varargin) +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 13 March 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +function [M ,iM ]= cummin(varargin) // Cumulative minimum // // Calling Sequence @@ -12,6 +23,8 @@ function M = cummin(varargin) // M = cummin(_,direction) // direction specifies as the direction of operation // + // [M , iM] = cummin(..) + // If called with two output arguments the index of the minimum value is also returned. // Parameters // A - real|complex numbers - vector|matrix // Input Array @@ -33,25 +46,17 @@ function M = cummin(varargin) // // Expected output: [8 8 1 1 1 1 1 1 1 1] // - // Authors - // Ayush Baid - // - // See Also - // cummax | cumprod | cumsum | max | min - - - [numOutArgs,numInArgs] = argn(0); // ** Checking number of arguments if numInArgs<1 | numInArgs>3 then - msg = "cummin: Wrong number of input argument; 1-6 expected"; + msg = "cummin: Wrong number of input argument; 1-3 expected"; error(77,msg); end - if numOutArgs~=1 then - msg = "cummin: Wrong number of output argument; 1 expected"; + if numOutArgs > 2 then + msg = "cummin: Wrong number of output argument; 1 or 2 expected"; error(78,msg); end @@ -103,11 +108,11 @@ function M = cummin(varargin) end // extracting direction - if strcmpi(directionArg,"reverse")==0 then + if strcmp(directionArg,"reverse")==0 then isForward = %f; - elseif strcmpi(directionArg,"forward")==0 then + elseif strcmp(directionArg,"forward")==0 then isForward = %t; - elseif strcmpi(directionArg,"")~=0 then + elseif strcmp(directionArg,"")~=0 then msg = "cummin: Wrong value for argument #3 (direction)"; error(53,msg); end @@ -132,7 +137,17 @@ function M = cummin(varargin) end M = matrix(M_,sizeA); - + if numOutArgs == 2 then + // calculating the index + // for vectors + iM = zeros(sizeA(1),sizeA(2)); + for i=1:sizeA(1) + for j=1:sizeA(2) + index = find (M(i,j) == A(i,:) ) + iM(i,j) = index(1) + end + end + end endfunction @@ -202,3 +217,19 @@ function out = cumminVec(inp,isForward) endfunction + +/* +# tests +w = cummin ([5 4 6 2 3 1]) // passed + +[w,iw] = cummin ([5 4 6 2 3 1]) // passed + +x = [1 2 3; 4 1 2; 3 5 1]; +result = cummin(x) // passed + +x = [1 2 3; 4 1 2; 3 5 1]; +result = cummin(x, 2) //passsed + +x = [1 2 3; 4 1 2; 3 5 1]; +[w,iw] = cummin(x, 2) //passsed +*/
\ No newline at end of file diff --git a/macros/deconv.sci b/macros/deconv.sci index 1c4d8fa..08ea92b 100644 --- a/macros/deconv.sci +++ b/macros/deconv.sci @@ -1,104 +1,106 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in function [b, r] = deconv (y, a) -// calling sequence: -// [b,r]= deconv (y, a) -// Deconvolve two vectors. -// -// [b, r] = deconv (y, a) solves for b and r such that -// y = conv (a, b) + r. -// -// If y and a are polynomial coefficient vectors, b will -// contain the coefficients of the polynomial quotient and r will be -// a remainder polynomial of lowest order. -//Test cases: -//1. -//[b, r] = deconv ([3, 6, 9, 9], [1, 2, 3]) -//Output: -//b=[3, 0] -//r=[0, 0, 0, 9] + if (nargin ~= 2) + error("deconv : Two arguments are required "); + end -//2. -//[b, r] = deconv ([3, 6], [1; 2; 3]) -//Output: -//b = 0. -//r= [- 2. ; 8] + if (~isvector(y) && ~isscalar(y)) + error ("deconv: Y must be vector"); + end +if ( ~isvector(a) && ~isscalar(a)) + error ("deconv: A must be vector"); +end + // Ensure A is oriented as Y. + if ((size(y,1)==1 && size(a,2)==1) || (size(y,2)==1 && size(a,1)==1)) + a = a.'; + end + + la = length (a); + ly = length (y); + + lb = ly - la + 1; + + if (ly > la) + x = zeros (size (y,1) - size (a,1) + 1,size(y,2)-size(a,2)+1); + x(1) = 1; + [b, r] = filter (y, a, x); + r = r * a(1); + elseif (ly == la) + [b, r] = filter (y, a, 1); + r = r * a(1); + else + b = 0; + r = y; + end + + if (nargout() > 1) + if (ly >= la) + r = [zeros(ly - la + 1, 1); r(1:la - 1)]; + // Respect the orientation of Y + r = matrix (r, size (y,1),size(y,2)); + end + end +endfunction +/* + //test + [b, r] = deconv ([3, 6, 9, 9], [1, 2, 3]); + assert_checkequal (b, [3, 0]); + assert_checkequal (r, [0, 0, 0, 9]); + //test + [b, r] = deconv ([3, 6], [1, 2, 3]); + assert_checkequal (b, 0); + assert_checkequal (r, [3, 6]); + //test + [b, r] = deconv ([3, 6], [1; 2; 3]); + assert_checkequal (b, 0); + assert_checkequal (r, [3, 6]); - [nargout,nargin]=argn(); +//test + [b,r] = deconv ([3; 6], [1; 2; 3]); + assert_checkequal (b, 0); + assert_checkequal (r, [3; 6]); - if (nargin ~= 2) - error ("wrong number of input arguments"); - end +//test + [b, r] = deconv ([3; 6], [1, 2, 3]); + assert_checkequal (b, 0); + assert_checkequal (r, [3; 6]); - if (~ (isvector (y) & isvector (a))) - error ("deconv: both arguments must be vectors"); - end + assert_checkequal (deconv ((1:3)',[1, 1]), [1; 1]) - la = length (a); - ly = length (y); +// Test input validation +// error deconv (1) +// error deconv (1,2,3) +// error <Y .* must be vector> deconv ([3, 6], [1, 2; 3, 4]) +// error <A must be vector> deconv ([3, 6], [1, 2; 3, 4]) - lb = ly - la + 1; +//test + y = (10:-1:1); + a = (4:-1:1); + [b, r] = deconv (y, a); + assert_checkequal (conv (a, b) + r, y); - // Ensure A is oriented as Y. - if (diff (size (y)) * diff (size (a)) < 0) - a = permute (a, [2, 1]); - end - if (ly > la) - o=size (y) - size (a) + 1; - x = zeros (o(1),o(2)); - x(1) = 1; - b = filter (real(y), real(a), x); - elseif (ly == la) - b = filter (real(y), real(a), 1); - else - b = 0; - end +//test + [b, r] = deconv ([1, 1], 1); + assert_checkequal (r, [0, 0]); - lc = la + length (b) - 1; - if (ly == lc) - if (length(a)==length(b) | length(a)>length(b)) - if isrow(a) - q=conv(a,b); - - u=[]; - for i=1:length(y) - u=[u;q]; - end - - w=[]; - if (isrow(y)) - for i=1:length(q) - w=[w;y]; - end - else - for i=1:length(q) - w=[w,y]; - end - end - - r = w-u; - - r=diag(r); - end - end - r=y-conv(a,b); - - -elseif(la~=lc) - // Respect the orientation of Y" - if (size (y,"r") <= size (y,"c")) - r = [(zeros (1, lc - ly)), y] - conv (a, b); - else - r = [(zeros (lc - ly, 1)); y] - conv (a, b); - end - if (ly < la) - // Trim the remainder is equal to the length of Y. - r = r($-(length(y)-1):$); - end -end +//test + [b, r] = deconv ([1; 1], 1); + assert_checkequal (r, [0; 0]); -endfunction + */ diff --git a/macros/digitrevorder.sci b/macros/digitrevorder.sci index 886c43f..0137542 100644 --- a/macros/digitrevorder.sci +++ b/macros/digitrevorder.sci @@ -1,36 +1,43 @@ -// Returns input data in digit-reversed order - -// Calling Sequence -//[y,i] = digitrevorder(x,r) -//y = digitrevorder(x,r) - -// Parameters -//x: Vector of real or complex values -//r: radix / base -//y: input vector in digit reverse order -//i: indices - -// Description -//This function returns the input data after reversing the digits of the indices and reordering the elements of the input array. - -// Examples -//x = [%i,1,3,6*%i] ; -//r = 2 ; -//[y i]=digitrevorder(x, r) -//Output : -// i = -// -// 1. 3. 2. 4. -// y = -// -// i 3. 1. 6.i +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in function [y, i] = digitrevorder (x, r) - funcprot(0); - [nargout, nargin] = argn() ; + // Returns input data in digit-reversed order + // Calling Sequence + //[y,i] = digitrevorder(x,r) + //y = digitrevorder(x,r) + // Parameters + //x: Vector of real or complex values + //r: radix / base + //y: input vector in digit reverse order + //i: indices + // Description + //This function returns the input data after reversing the digits of the indices and reordering the elements of the input array. + // Examples + //x = [%i,1,3,6*%i] ; + //r = 2 ; + //[y i]=digitrevorder(x, r) + //Output : + // i = + // + // 1. 3. 2. 4. + // y = + // + // i 3. 1. 6.i + funcprot(0); + [nargout, nargin] = argn() ; if (nargin > 2 | nargin <= 1) error("digitrevorder : invalid number of inputs") - elseif (~ isvector (x)) + elseif ( ~isvector(x) & ~isscalar(x)) error ("digitrevorder : X must be a vector"); elseif (~ (isscalar (r) & r == fix (r) & r >= 2 & r <= 36)) error ("digitrevorder : R must be an integer between 2 and 36"); @@ -41,27 +48,37 @@ function [y, i] = digitrevorder (x, r) end end - old_ind = 0:length(x) - 1; - - //new_ind = base2dec(mtlb_fliplr(dec2base(old_ind, r)), r); //it works only on octave - old_ind_base = dec2base(old_ind, r) ; - new_ind_base = [] ; - b = [] ; - for i=1:length(x) - new_ind_base = [new_ind_base mtlb_fliplr(old_ind_base(i))]; - end - new_ind = base2dec(new_ind_base,r) ; - //end of index conversion - + old_ind = 0:max(size(x)) - 1; + new_ind = base2dec (strrev(dec2base (old_ind, r)), r); i = new_ind + 1; y(old_ind + 1) = x(i); - - [rows_x columns_x] = size(x) ; - - if (columns_x == 1) - y = y'; + if (size(x,2)== 1) + y = y(:); else - i = i; + i = i.'; end - endfunction + +/* +// tests base 0 2 OK +assert_checkequal (digitrevorder (0, 2), 0); // passed +assert_checkequal (digitrevorder (0, 36), 0); //passed +assert_checkequal (digitrevorder (0:3, 4), 0:3); //passed +assert_checkequal (digitrevorder ([0:3]', 4), [0:3]');//passed +assert_checkequal (digitrevorder (0:7, 2), [0 4 2 6 1 5 3 7]);//passed +assert_checkequal (digitrevorder ([0:7]', 2), [0 4 2 6 1 5 3 7]');//passed +assert_checkequal (digitrevorder ([0:7]*%i, 2), [0 4 2 6 1 5 3 7]*%i); //passed +assert_checkequal (digitrevorder ([0:7]'*%i, 2), [0 4 2 6 1 5 3 7]'*%i);//passed +assert_checkequal (digitrevorder (0:15, 2), [0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15]); // passed + +//FIXME : debug the failed test +assert_checkequal (digitrevorder (0:15, 4), [0 4 8 12 1 5 9 13 2 6 10 14 3 7 11 15]); // failed + +//Test input validation - passed +error digitrevorder (); +error digitrevorder (1); +error digitrevorder (1, 2, 3); +error digitrevorder ([], 1); +error digitrevorder ([], 37); +error digitrevorder (0:3, 8); +*/
\ No newline at end of file diff --git a/macros/ellip.sci b/macros/ellip.sci index ab34458..1dc8435 100644 --- a/macros/ellip.sci +++ b/macros/ellip.sci @@ -1,15 +1,14 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt // Original Source : https://octave.sourceforge.io/signal/ -// Modifieded by:Sonu Sharma, RGIT Mumbai +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - function [a, b, c, d] = ellip (n, rp, rs, w, varargin) //Elliptic or Cauer filter design with rp dB of passband ripple and rs dB of stopband attenuation. @@ -59,7 +58,8 @@ function [a, b, c, d] = ellip (n, rp, rs, w, varargin) // column 5 // // 0.0202774 - + // Dependencies + // ellipap sftrans bilinear zp2tf funcprot(0); [nargout nargin] = argn(); @@ -85,7 +85,7 @@ function [a, b, c, d] = ellip (n, rp, rs, w, varargin) case "stop" stop = %T; case "low" - stop = %T; + stop = %F; case "pass" stop = %F; else @@ -93,7 +93,7 @@ function [a, b, c, d] = ellip (n, rp, rs, w, varargin) end end - [rows_w colums_w] = size(w); + [rows_w columns_w] = size(w); if (~ ((length (w) <= 2) & (rows_w == 1 | columns_w == 1))) error ("ellip: frequency must be given as WP or [WL, WH]"); @@ -119,7 +119,7 @@ function [a, b, c, d] = ellip (n, rp, rs, w, varargin) // Prewarp the digital frequencies if (digital) T = 2; // sampling frequency of 2 Hz - w = 2 / T * tan (%pi * w / T); + w = (2 / T )* tan (%pi * w / T); end // Generate s-plane poles, zeros and gain @@ -139,8 +139,8 @@ function [a, b, c, d] = ellip (n, rp, rs, w, varargin) if (nargout == 2) [a b] = zp2tf(zero, pole, gain); elseif (nargout == 3) - a = zero; - b = pole; + a = conj(zero'); + b = conj(pole'); c = gain; else // output ss results @@ -148,3 +148,49 @@ function [a, b, c, d] = ellip (n, rp, rs, w, varargin) error("ellip: yet not implemented in state-space form OR invalid number of o/p arguments") end endfunction + +/* +%% Test input validation // all passed +error [a, b] = ellip () +error [a, b] = ellip (1) +error [a, b] = ellip (1, 2) +error [a, b] = ellip (1, 2, 3) +error [a, b] = ellip (1, 2, 3, 4, 5, 6, 7) +error [a, b] = ellip (.5, 2, 40, .2) +error [a, b] = ellip (3, 2, 40, .2, "invalid") + +test + +//[b, a] = ellip (n, rp, rs, [wl, wh]) +[b, a] = ellip (5, 1, 90, [.1 .2]); // passed + +//[b, a] = ellip (n, rp, rs, wp) +[A, B] = ellip (6, 3, 50, .6); //passed + +//[b, a] = ellip (n, rp, rs, [wl, wh], "stop") +[b, a] = ellip (5, 1, 90, [.1 .2],"stop"); // passed + +//[b, a] = ellip (n, rp, rs, wp, "high") +[A, B] = ellip (6, 3, 50, .6,"high"); // passed + +//[b, a] = ellip (n, rp, rs, [wl, wh]) +[b, a] = ellip (5, 1, 90, [.1 .2],"s"); // passed + +//[b, a] = ellip (n, rp, rs, wp) +[A, B] = ellip (6, 3, 50, .6,"z"); // passed + +//[b, a] = ellip (n, rp, rs, [wl, wh], "stop") +[b, a] = ellip (5, 1, 90, [.1 .2],"pass"); // passed + +//[b, a] = ellip (n, rp, rs, wp, "high") +[A, B] = ellip (6, 3, 50, .6,"low"); // passed + + +//[z, p, g] = ellip (…) +[z, p, g] = ellip (6, 3, 50, .6); // passed + +[z,p,g] = ellip( 10,15,4,.9,"s") // passed + +test + [a, b, c, d] = ellip (6, 3, 50, .6); // not implemented yet + */
\ No newline at end of file diff --git a/macros/fftfilt.sci b/macros/fftfilt.sci index ab171ce..1f696ab 100644 --- a/macros/fftfilt.sci +++ b/macros/fftfilt.sci @@ -1,3 +1,14 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in function y = fftfilt(b, x, varargin) // Performs FFT-based FIR filtering using overlap-add method // @@ -31,13 +42,7 @@ function y = fftfilt(b, x, varargin) // x = sin(1:2000); // b = [1 1/4;1/3 1/5]; // y = fftfilt(b,x); - // - // Authors - // Ayush Baid - - - - + [numOutArgs,numInArgs] = argn(0); // ** Checking number of arguments @@ -201,4 +206,4 @@ function y = fftfilt(b, x, varargin) end -endfunction +endfunction
\ No newline at end of file diff --git a/macros/findpeaks.sci b/macros/findpeaks.sci index 3e7ce47..a96dff7 100644 --- a/macros/findpeaks.sci +++ b/macros/findpeaks.sci @@ -1,58 +1,330 @@ -function [PKS, LOC, EXTRA] = findpeaks(DATA, varargin) -//This function find peaks on DATA. -//Calling Sequence -//[PKS, LOC, EXTRA] = findpeaks(DATA) -//[PKS, LOC, EXTRA] = findpeaks(..., PROPERTY, VALUE) -//[PKS, LOC, EXTRA] = findpeaks(..., "DoubleSided") -//Description -//Peaks of a positive array of data are defined as local maxima. For double-sided data, they are maxima of the positive part and minima of the negative part. DATA is expected to be a single column vector. -// -//The function returns the value of DATA at the peaks in PKS. The index indicating their position is returned in LOC. -// -//The third output argument is a structure with additional information: -// -//"parabol": -// A structure containing the parabola fitted to each returned peak. The structure has two fields, "x" and "pp". The field "pp" contains the coefficients of the 2nd degree polynomial and "x" the extrema of the intercal here it was fitted. -// -//"height": -// The estimated height of the returned peaks (in units of DATA). -// -//"baseline": -// The height at which the roots of the returned peaks were calculated (in units of DATA). -// -//"roots": -// The abscissa values (in index units) at which the parabola fitted to each of the returned peaks crosses the "baseline" value. The width of the peak is calculated by 'diff(roots)'. -// -//This function accepts property-value pair given in the list below: -// -//"MinPeakHeight": -// Minimum peak height (positive scalar). Only peaks that exceed this value will be returned. For data taking positive and negative values use the option "DoubleSided". Default value '2*std (abs (detrend (data,0)))'. -// -//"MinPeakDistance": -// Minimum separation between (positive integer). Peaks separated by less than this distance are considered a single peak. This distance is also used to fit a second order polynomial to the peaks to estimate their width, therefore it acts as a smoothing parameter. Default value 4. -// -//"MinPeakWidth": -// Minimum width of peaks (positive integer). The width of the peaks is estimated using a parabola fitted to the neighborhood of each peak. The neighborhood size is equal to the value of "MinPeakDistance". The width is evaluated at the half height of the peak with baseline at "MinPeakHeight". Default value 2. -// -//"DoubleSided": -// Tells the function that data takes positive and negative values. The base-line for the peaks is taken as the mean value of the function. This is equivalent as passing the absolute value of the data after removing the mean. -funcprot(0); -rhs=argn(2); -lhs=argn(1) -if(rhs<1 | rhs>2) then - error("Wrong number of input arguments."); -end -if(lhs<3 | lhs>3) then - error("Wrong number of output arguments."); -end +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +function [pks ,idx, varargout] = findpeaks (data, varargin) +/* + + [pks, loc, extra] = findpeaks (data) + … = findpeaks (…, property, value) + … = findpeaks (…, "DoubleSided") + Finds peaks on data. + + Peaks of a positive array of data are defined as local maxima. For double-sided data, they are maxima of the positive part and minima of the negative part. data is expected to be a single column vector. + + The function returns the value of data at the peaks in pks. The index indicating their position is returned in loc. + + The third output argument is a structure with additional information: + + "parabol" + A structure containing the parabola fitted to each returned peak. The structure has two fields, "x" and "pp". The field "pp" contains the coefficients of the 2nd degree polynomial and "x" the extrema of the interval where it was fitted. + + "height" + The estimated height of the returned peaks (in units of data). + + "baseline" + The height at which the roots of the returned peaks were calculated (in units of data). + + "roots" + The abscissa values (in index units) at which the parabola fitted to each of the returned peaks realizes its width as defined below. + + This function accepts property-value pair given in the list below: + + "MinPeakHeight" + Minimum peak height (non-negative scalar). Only peaks that exceed this value will be returned. For data taking positive and negative values use the option "DoubleSided". Default value eps. + + "MinPeakDistance" + Minimum separation between (positive integer). Peaks separated by less than this distance are considered a single peak. This distance is also used to fit a second order polynomial to the peaks to estimate their width, therefore it acts as a smoothing parameter. The neighborhood size is equal to the value of "MinPeakDistance". Default value 1. + + "MinPeakWidth" + Minimum width of peaks (positive integer). The width of the peaks is estimated using a parabola fitted to the neighborhood of each peak. The width is calculated with the formula + + a * (width - x0)^2 = 1 + where a is the concavity of the parabola and x0 its vertex. Default value 1. + + "MaxPeakWidth" + Maximum width of peaks (positive integer). Default value Inf. + + "DoubleSided" + Tells the function that data takes positive and negative values. The base-line for the peaks is taken as the mean value of the function. This is equivalent as passing the absolute value of the data after removing the mean. +*/ + + if (nargin < 1) + error("findpeaks:InsufficientInputArguments\nfindpeaks: DATA must be given"); + end + + if (~(isvector (data) && length (data) >= 3)) + error ("findpeaks:InvalidArgument\nfindpeaks: DATA must be a vector of at least 3 elements"); + end + transpose = (size(data,1) == 1); + if (transpose) + data = data.'; + end + __data__ = abs (detrend (data, 'c')); + // --- Parse arguments --- // +[ dSided, minH, minD, minW, maxW ] = parser ( varargin(:) ); + + if (dSided) + temp = __data__ + __data__ = data + data = temp + temp = [] // free temp + elseif (min (data) < 0) + error ("findpeaks:InvalidArgument\nData contains negative values. You may want to DoubleSided option"); + end + // Rough estimates of first and second derivative + df1 = diff (data, 1) + df1=df1([1; [1:length(df1)].']); + df2 = diff (data, 2) + df2=df2([1; 1; [1:length(df2)].']); + + // check for changes of sign of 1st derivative and negativity of 2nd + // derivative. + // <= in 1st derivative includes the case of oversampled signals. + idx = find (df1.*[df1(2:$); 0] <= 0 & [df2(2:$); 0] < 0); + + // Get peaks that are beyond given height + tf = data(idx) > minH; + idx = idx(tf); + + // sort according to magnitude + [_, tmp] = gsort (data(idx)); + idx_s = idx(tmp); + + // Treat peaks separated less than minD as one + D = abs (bsminuseq (idx_s)); + D = D + diag(ones(1,size(D,1))*%nan); // eliminate diagonal cpmparison + + if (or(D(:) < minD)) // FIXME : this branch is not tested + i = 1; + peak = cell (); + node2visit = 1:size(D,1); + visited = []; + idx_pruned = idx_s; + + // debug +// h = plot(1:length(data),data,"-",idx_s,data(idx_s),'.r',idx_s,data(idx_s),'.g'); +// set(h(3),"visible","off"); + + while (~isempty (node2visit)) + + d = D(node2visit(1),:); + + visited = [visited node2visit(1)]; + node2visit(1) = []; + + neighs = setdiff (find (d < minD), visited); + if ( ~isempty (neighs)) + // debug +// set(h(3),"xdata",idx_s(neighs),"ydata",data(idx_s(neighs)),"visible","on") +// pause(0.2) +// set(h(3),"visible","off"); + + idx_pruned = setdiff (idx_pruned, idx_s(neighs)); + + visited = [visited neighs]; + node2visit = setdiff (node2visit, visited); + + // debug +// set(h(2),"xdata",idx_pruned,"ydata",data(idx_pruned)) +// pause + end + + end + idx = idx_pruned; + end + + extra = struct ("parabol", [], "height", [], "baseline", [], "roots", []); + + // Estimate widths of peaks and filter for: + // width smaller than given. + // wrong concavity. + // not high enough + // data at peak is lower than parabola by 1% + // position of extrema minus center is bigger equal than minD/2 + // debug +// h = plot(1:length(data),data,"-",idx,data(idx),'.r',... +// idx,data(idx),'og',idx,data(idx),'-m'); +// set(h(4),"linewidth",2) +// set(h(3:4),"visible","off"); + + idx_pruned = idx; + n = length (idx); + np = length (data); + struct_count = 0; + + for i=1:n + ind = (floor (max(idx(i)-minD/2,1)) : ... + ceil (min(idx(i)+minD/2,np))).'; + pp = zeros (1,3); + // If current peak is not local maxima, then fit parabola to neighbor + if or(data(idx(i)-1) == data(idx(i))) + // sample on left same as peak + xm = 0; + pp = ones (1,3); + elseif or(data(ind) > data(idx(i))) + pp = polyfit (ind, data(ind), 2); + xm = -pp(2)^2 / (2*pp(1)); // position of extrema + H = polyval (pp, xm); // value at extrema + else // use it as vertex of parabola + H = data(idx(i)); + xm = idx(i); + pp = zeros (1,3); + pp(1) = pinv((ind-xm).^2 ) * (data(ind)-H); + pp(2) = - 2 * pp(1) * xm; + pp(3) = H + pp(1) * xm^2; + end + + // debug +// x = linspace(ind(1)-1,ind(end)+1,10); +// set(h(4),"xdata",x,"ydata",polyval(pp,x),"visible","on") +// set(h(3),"xdata",ind,"ydata",data(ind),"visible","on") +// pause(0.2) +// set(h(3:4),"visible","off"); + +// thrsh = min (data(ind([1 end]))); +// rz = roots ([pp(1:2) pp(3)-thrsh]); +// width = abs (diff (rz)); + width = sqrt (abs(1 / pp(1))); + + if ( (width > maxW || width < minW) || pp(1) > 0 || H < minH || data(idx(i)) < 0.99*H ||abs (idx(i) - xm) > minD/2) then + idx_pruned = setdiff (idx_pruned, idx(i)); + elseif (nargout > 2) + struct_count=struct_count+1; + extra.parabol(struct_count).x = ind([1 $]); + extra.parabol(struct_count).pp = pp; + + extra.roots(struct_count,1:2)= xm + [-width width]/2; + extra.height(struct_count) = H; + extra.baseline(struct_count) = mean ([H minH]); + end + + // debug +// set(h(2),"xdata",idx_pruned,"ydata",data(idx_pruned)) +// pause(0.2) -select(rhs) -case 1 then - [PKS, LOC, EXTRA] = callOctave("findpeaks", DATA); -case 2 then - [PKS, LOC, EXTRA] = callOctave("findpeaks", DATA, varargin(1)); -case 3 then - [PKS, LOC, EXTRA] = callOctave("findpeaks", DATA, varargin(1), varargin(2)); end + + + idx = idx_pruned; + if (dSided) + pks = __data__(idx); + else + pks = data(idx); + end + + if (transpose) + pks = pks.'; + idx = idx.'; + end + idx = idx.'; + if (nargout() > 2) + varargout(1) = extra; + end + endfunction +/* +demo +pi = %pi ; + t = 2*pi*linspace(0,1,1024)'; + y = sin(3.14*t) + 0.5*cos(6.09*t) + 0.1*sin(10.11*t+1/6) + 0.1*sin(15.3*t+1/3); + data1 = abs(y); + [pks idx] = findpeaks(data1); + data2 = y; + [pks2 idx2] = findpeaks(data2,"DoubleSided"); + [pks3 idx3] = findpeaks(data2,"DoubleSided","MinPeakHeight",0.5); + + subplot(1,2,1) + plot(t,data1,t(idx),data1(idx),'xm') + subplot(1,2,2) + plot(t,data2,t(idx2),data2(idx2),"xm",t(idx3),data2(idx3),"or") + legend("Location","NorthOutside","Orientation","horizontal"); +/////////////////////////////// octave version ///////////////////////////// + t = 2*pi*linspace(0,1,1024)'; + y = sin(3.14*t) + 0.5*cos(6.09*t) + 0.1*sin(10.11*t+1/6) + 0.1*sin(15.3*t+1/3); + data1 = abs(y); + [pks idx] = findpeaks(data1); + data2 = y; + [pks2 idx2] = findpeaks(data2,"DoubleSided"); + [pks3 idx3] = findpeaks(data2,"DoubleSided","MinPeakHeight",0.5); + subplot(1,2,1) + plot(t,data1,t(idx),data1(idx),'xm') + axis tight + subplot(1,2,2) + plot(t,data2,t(idx2),data2(idx2),"xm",t(idx3),data2(idx3),"or") + axis tight + legend("Location","NorthOutside","Orientation","horizontal"); +////////////////////octvae version end //////////////////////////////////////////// +assert_checkequal(size(pks),[11 1]); +assert_checkequal(size(idx),[11 1]); +assert_checkequal(size(pks2),[11 1]); +assert_checkequal(size(idx2),[11 1]); +assert_checkequal(size(pks3),[8 1]); +assert_checkequal(size(idx3),[8 1]); + //---------------------------------------------------------------------------- + + Finding the peaks of smooth data is not a big deal! +// Not as accurate as Octave +demo + t = 2*pi*linspace(0,1,1024)'; + y = sin(3.14*t) + 0.5*cos(6.09*t) + 0.1*sin(10.11*t+1/6) + 0.1*sin(15.3*t+1/3); + data = abs(y + 0.1*rand(length(y),1,'normal')); + [pks idx] = findpeaks(data,"MinPeakHeight",1); + dt = t(2)-t(1); + [pks2 idx2] = findpeaks(data,"MinPeakHeight",1,"MinPeakDistance",round(0.5/dt)); + subplot(1,2,1) + plot(t,data,t(idx),data(idx),'or') + subplot(1,2,2) + plot(t,data,t(idx2),data(idx2),'or') + + ///octave version/////////////////////////////////////////////////////// + t = 2*pi*linspace(0,1,1024)'; + y = sin(3.14*t) + 0.5*cos(6.09*t) + 0.1*sin(10.11*t+1/6) + 0.1*sin(15.3*t+1/3); + data = abs(y + 0.1*randn(length(y),1)); + [pks idx] = findpeaks(data,"MinPeakHeight",1); + dt = t(2)-t(1); + [pks2 idx2] = findpeaks(data,"MinPeakHeight",1,"MinPeakDistance",round(0.5/dt)); + subplot(1,2,1) + plot(t,data,t(idx),data(idx),'or') + subplot(1,2,2) + plot(t,data,t(idx2),data(idx2),'or') + ///////////////////////////////////DO not run enclosed on scilab////////////////////////////////////////////////////// + assert_checkequal(size(pks),[86 1]); + assert_checkequal(size(pks2),[7 1]); + assert_checkequal(size(idx),[86 1]); + assert_checkequal(size(idx2),[7 1]); + + //---------------------------------------------------------------------------- + // Noisy data may need tuning of the parameters. In the 2nd example, + // MinPeakDistance is used as a smoother of the peaks. + +assert_checkequal (findpeaks ([1, 1, 1]),[]) +assert_checkequal (findpeaks ([1; 1; 1]),[]) + + +test + // Test input vector is an oversampled sinusoid with clipped peaks + x = min (3, cos (2*pi*[0:8000] ./ 600) + 2.01); + assert_checkequal (~isempty (findpeaks (x)),%t) + + +test + x = [1 10 2 2 1 9 1]; + [pks, loc] = findpeaks(x); + assert_checkequal (loc, [2 6]) + assert_checkequal (pks, [10 9]) + +// Test input validation +error findpeaks () +error findpeaks (1) +error findpeaks ([1, 2]) + +*/ diff --git a/macros/freqz.sci b/macros/freqz.sci index dc71b75..4579da3 100644 --- a/macros/freqz.sci +++ b/macros/freqz.sci @@ -1,71 +1,247 @@ -function [H, W] = freqz(B, varargin) -//This function returns the complex frequency response H of the rational IIR filter whose numerator and denominator coefficients are B and A, respectively. -//Calling Sequence -//[H, W] = freqz(B, A, N, "whole") -//[H, W] = freqz(B) -//[H, W] = freqz(B, A) -//[H, W] = freqz(B, A, N) -//H = freqz(B, A, W) -//[H, W] = freqz(..., FS) -//freqz(...) -//Parameters -//B, A, N: Integer or Vector -//Description -// Return the complex frequency response H of the rational IIR filter whose numerator and denominator coefficients are B and A, respectively. -// -//The response is evaluated at N angular frequencies between 0 and 2*pi. -// -//The output value W is a vector of the frequencies. -// -//If A is omitted, the denominator is assumed to be 1 (this corresponds to a simple FIR filter). -// -//If N is omitted, a value of 512 is assumed. For fastest computation, N should factor into a small number of small primes. -// -//If the fourth argument, "whole", is omitted the response is evaluated at frequencies between 0 and pi. -// -// 'freqz (B, A, W)' -// -//Evaluate the response at the specific frequencies in the vector W. The values for W are measured in radians. -// -// '[...] = freqz (..., FS)' -// -//Return frequencies in Hz instead of radians assuming a sampling rate FS. If you are evaluating the response at specific frequencies W, those frequencies should be requested in Hz rather than radians. -// -// 'freqz (...)' -// -//Plot the magnitude and phase response of H rather than returning them. -//Examples -//H = freqz([1,2,3], [4,3], [1,2,5]) -//ans = -// 0.4164716 - 0.5976772i - 0.4107690 - 0.2430335i 0.1761948 + 0.6273032i -funcprot(0); -rhs=argn(2); -lhs=argn(1); -if(rhs<2 | rhs>4) then - error("Wrong number of input arguments."); -end -if (lhs<1 | lhs>2) - error("Wrong number of output arguments."); -end -if (lhs==1) then -select(rhs) -case 1 then - H = callOctave("freqz",B); -case 2 then - H = callOctave("freqz",B, varargin(1)); -case 3 then - H = callOctave("freqz",B, varargin(1), varargin(2)); -end -elseif (lhs==2) then - select(rhs) -case 1 then - [H, W] = callOctave("freqz",B); -case 2 then - [H, W] = callOctave("freqz",B, varargin(1)); -case 3 then - [H, W] = callOctave("freqz",B, varargin(1), varargin(2)); -case 4 then - [H, W] = callOctave("freqz", B, varargin(1), varargin(2), varargin(3)); -end -end +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh , Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +// Example usage of freqz in different formats: +// [h, w] = freqz (b, a, n, "whole") +// [h, w] = freqz (b) +// [h, w] = freqz (b, a) +// [h, w] = freqz (b, a, n) +// h = freqz (b, a, w) +// [h, w] = freqz (…, Fs) +// freqz (...) + +/* +Return the complex frequency response `h` of the rational IIR filter whose +numerator and denominator coefficients are `b` and `a`, respectively. + +The response is evaluated at `n` angular frequencies between 0 and 2*pi. + +The output value `w` is a vector of the frequencies. + +If `a` is omitted, the denominator is assumed to be 1 (this corresponds to a simple FIR filter). + +If `n` is omitted, a value of 512 is assumed. For fastest computation, `n` should factor into a small number of small primes. + +If the fourth argument, "whole", is omitted, the response is evaluated at frequencies between 0 and pi. + +freqz (b, a, w) + +Evaluate the response at the specific frequencies in the vector `w`. The values for `w` are measured in radians. + +[...] = freqz (…, Fs) + +Return frequencies in Hz instead of radians assuming a sampling rate `Fs`. If you are evaluating the response at specific frequencies `w`, those frequencies should be requested in Hz rather than radians. + +freqz (...) + +Plot the magnitude and phase response of `h` rather than returning them. +*/ +// Dependencies +// fft1 unwrap2 postpad +function [h_r, f_r] = freqz (b, a, n, region, Fs) + + if (nargin < 1) + error("Invalid numbers of inputs"); + elseif (nargin == 1) + // Response of an FIR filter. + a = []; + n = []; + region =[]; + Fs = []; + elseif (nargin == 2) + // Response of an IIR filter + n = []; + region = []; + Fs = []; + elseif (nargin == 3) + region =[]; + Fs = []; + elseif (nargin == 4) + Fs = []; + if (~ (type(region)==10) && ~ isempty(region)) + Fs = region; + region = []; + end + end + + if (isempty (b)) + b = 1; + elseif (~ isvector (b)) + error ("freqz: B must be a vector"); + end + if (isempty (a)) + a = 1; + elseif (~ isvector (a)) + error ("freqz: A must be a vector"); + end + if (isempty (n)) + n = 512; + elseif (isscalar (n) && n < 1) + error ("freqz: N must be a positive integer"); + end + if (isempty (region)) + if (isreal (b) && isreal (a)) + region = "half"; + else + region = "whole"; + end + end + if (isempty (Fs)) + freq_norm = %t; + if (nargout == 1) + Fs = 2; + else + Fs = 2*%pi; + end + else + freq_norm = %f; + end + // FIXME : nargout != 0 even if no output parameter is used + // FIXME : problem in argn() or nargin function + //plot_output = (nargout == 0); + plot_output = (nargout == 1) // temp fix + whole_region = ~strcmp (region, "whole"); + + a = a(:); + b = b(:); + + if (~ isscalar (n)) + // Explicit frequency vector given + w = n; + f = n; + if (nargin == 4) + // Sampling rate Fs was specified + w = 2*%pi*f/Fs; + end + k = max (length (b), length (a)); + hb = polyval (postpad (b, k), exp (%i*w)); + ha = polyval (postpad (a, k), exp (%i*w)); + else + // polyval(fliplr(P),exp(jw)) is O(p n) and fft(x) is O(n log(n)), + // where p is the order of the polynomial P. For small p it + // would be faster to use polyval but in practice the overhead for + // polyval is much higher and the little bit of time saved isn't + // worth the extra code. + k = max (length (b), length (a)); + if (k > n/2 && nargout == 0) + // Ensure a causal phase response. + n = n * 2 .^ ceil (log2 (2*k/n)); + end + + if (whole_region) + N = n; + if (plot_output) + f = Fs * (0:n).' / N; // do 1 more for the plot + else + f = Fs * (0:n-1).' / N; + end + else + N = 2*n; + if (plot_output) + n = n + 1; + end + f = Fs * (0:n-1).' / N; + end + + pad_sz = N*ceil (k/N); + b = postpad (b, pad_sz); + a = postpad (a, pad_sz); + + hb = zeros (n, 1); + ha = zeros (n, 1); + + for i = 1:N:pad_sz + fftresult = fft1 (postpad (b(i:i+N-1), N))(1:n); + if size(fftresult,1) == 1 then fftresult = fftresult';end + hb = hb + fftresult ; + tempfftresult=fft1 (postpad (a(i:i+N-1), N))(1:n); + if size(tempfftresult,1) == 1 then tempfftresult = tempfftresult';end + ha = ha + tempfftresult; + end + + end + + h = hb ./ ha; + + if (plot_output) + // Plot and don't return values. + if (whole_region && isscalar (n)) + h($+1) = h(1); // Solution is periodic. Copy first value to end. + end + freqz_plot (f, h, freq_norm); + end + // Return values and don't plot. + h_r = h; + f_r = f; + + endfunction +function freqz_plot (w, h, freq_norm) + if (nargin < 2) + error("Invalid numbers of inputs"); + end + + if nargin < 3 then + freq_norm = %f + end + n = size(max(w)); + mag = 20 * log10 (abs (h)); + phase = unwrap2 (angle (h)); + + if (freq_norm) + x_label = 'Normalized Frequency (\times\pi rad/sample)'; + else + x_label = "Frequency (Hz)"; + end + subplot (2, 1, 1); + plot (w, mag); + xgrid; + xlabel (x_label); + ylabel ("Magnitude (dB)"); + + subplot (2, 1, 2); + plot (w, phase*360/(2*%pi)); + xgrid; + xlabel (x_label); + ylabel ("Phase (degrees)"); + +endfunction +/* +// passed +testif HAVE_FFTW # correct values and fft-polyval consistency + // butterworth filter, order 2, cutoff pi/2 radians + b = [0.292893218813452 0.585786437626905 0.292893218813452]; + a = [1 0 0.171572875253810]; + [h,w] = freqz (b,a,32); + +//passed +testif HAVE_FFTW # whole-half consistency + b = [1 1 1]/3; + [h,w] = freqz (b,1,32,"whole"); + + [h2,w2] = freqz (b,1,16,"half"); + + +//passed +testif HAVE_FFTW # Sampling frequency properly interpreted + b = [1 1 1]/3; a = [1 0.2]; + [h,f] = freqz (b,a,16,320); + + [h2,f2] = freqz (b,a,[0:15]*10,320); + + [h3,f3] = freqz (b,a,32,"whole",320); + + +// Test input validation +// FIXME: Need to put tests here and simplify input validation in the main code. + +*/ diff --git a/macros/h1_z_deriv.sci b/macros/h1_z_deriv.sci new file mode 100644 index 0000000..65e87e8 --- /dev/null +++ b/macros/h1_z_deriv.sci @@ -0,0 +1,47 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +function b = h1_z_deriv(n, p, ts) + + // Build the vector d that holds coefficients for all the derivatives of H1(z) + // The results reads d(n)*z^(1)*(d/dz)^(1)*H1(z) + d(n-1)*z^(2)*(d/dz)^(2)*H1(z) +...+ d(1)*z^(n)*(d/dz)^(n)*H1(z) + d = (-1)^n; // Vector with the derivatives of H1(z) + for i= (1:n-1) + d = [d 0]; // Shift result right (multiply by -z) + d = d + prepad(polyder(d), i+1, 0, 2); // Add the derivative + end + + // Build output vector + b = zeros (1, n + 1); + for i = (1:n) + b = b + d(i) * prepad(h1_deriv(n-i+1), n+1, 0, 2); + end + b = b * ts^(n+1)/factorial(n); + // Multiply coefficients with p^i, where i is the index of the coeff. + b = b.*p.^(n+1:-1:1); +endfunction + +// Find (z^n)*(d/dz)^n*H1(z), where H1(z)=ts*z/(z-p), ts=sampling period, +// p=exp(sm*ts) and sm is the s-domain pole with multiplicity n+1. +// The result is (ts^(n+1))*(b(1)*p/(z-p)^1 + b(2)*p^2/(z-p)^2 + b(n+1)*p^(n+1)/(z-p)^(n+1)), +// where b(i) is the binomial coefficient bincoeff(n,i) times n!. Works for n>0. + +function b = h1_deriv(n) + b = factorial(n)*nchoosek(n,0:n); // Binomial coefficients: [1], [1 1], [1 2 1], [1 3 3 1], etc. + b = b*(-1)^n; +endfunction + +function y = polyder(a) + y = poly(flipdim(a,2),'a','coeff') + y = derivat(y) + y = coeff(y) +endfunction + diff --git a/macros/iirlp2mb.sci b/macros/iirlp2mb.sci index 54f1a6b..09ec743 100644 --- a/macros/iirlp2mb.sci +++ b/macros/iirlp2mb.sci @@ -1,41 +1,345 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +// IIR Low Pass Filter to Multiband Filter Transformation +// +// [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt) +// [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt,Pass) +// +// Num,Den: numerator,denominator of the transformed filter +// AllpassNum,AllpassDen: numerator,denominator of allpass transform, +// B,A: numerator,denominator of prototype low pass filter +// Wo: normalized_angular_frequency/pi to be transformed +// Wt: [phi=normalized_angular_frequencies]/pi target vector +// Pass: This parameter may have values 'pass' or 'stop'. If +// not given, it defaults to the value of 'pass'. +// +// With normalized ang. freq. targets 0 < phi(1) < ... < phi(n) < pi radians +// +// for Pass == 'pass', the target multiband magnitude will be: +// -------- ---------- -----------... +// / \ / \ / . +// 0 phi(1) phi(2) phi(3) phi(4) phi(5) (phi(6)) pi +// +// for Pass == 'stop', the target multiband magnitude will be: +// ------- --------- ----------... +// \ / \ / . +// 0 phi(1) phi(2) phi(3) phi(4) (phi(5)) pi +// function [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(varargin) -//This function does IIR Low Pass Filter to Multiband Filter Transformation. -//Calling Sequence -//[Num, Den, AllpassNum, AllpassDen] = iirlp2mb(B, A, Wo, Wt) -//[Num, Den, AllpassNum, AllpassDen] = iirlp2mb(B, A, Wo, Wt, Pass) -//Parameters -//B: real or complex value -//A: real or complex value -//Wo: scalar or vector -//Wt: scalar or vector, elements must be monotonically increasing and >= 0 and <= 1. -//Description -//This is an Octave function. -//This function does IIR Low Pass Filter to Multiband Filter Transformation. -//The first two parameters give the numerator and denominator of the prototype low pass filter. -//The third parameter is the normalized angular frequency/pi to be transformed. -//The fourth parameter is the normalized angular frequency/pi target vector. -//The first two output variables are the numerator and denominator of the transformed filter. -//The third and fourth output variables are the numerator and denominator of the allpass transform. -//The fifth parameter can have values pass or stop, default value is pass. -//Examples -//iirlp2mb(5,9,0.3,0.4) -//ans = 0.55556 - -B = varargin(1) -A = varargin(2) -Wo = varargin(3) -Wt = varargin(4) -rhs = argn(2) -lhs = argn(1) -if(rhs < 4 | rhs > 5) -error("Wrong number of input arguments.") -end -select(rhs) -case 4 then -[Num,Den,AllpassNum,AllpassDen] = callOctave("iirlp2mb",varargin(1),varargin(2),varargin(3),varargin(4)) -case 5 then -[Num,Den,AllpassNum,AllpassDen] = callOctave("iirlp2mb",varargin(1),varargin(2),varargin(3),varargin(4),varargin(5)) -end + + usage = sprintf("iirlp2mb Usage: [Num,Den,AllpassNum,AllpassDen]=iirlp2mb(B,A,Wo,Wt[,Pass])\n"); + B = varargin(1); // numerator polynomial of prototype low pass filter + A = varargin(2); // denominator polynomial of prototype low pass filter + Wo = varargin(3); // (normalized angular frequency)/pi to be transformed + Wt = varargin(4); // vector of (norm. angular frequency)/pi transform targets + // [phi(1) phi(2) ... ]/pi + if(nargin < 4 || nargin > 5) + error(usage) + end + + // Checking proper input arguments + if ~isvector(B) then + if ~isscalar(B) then // if numerator is having only one coeffcient then it can be passed directly + error(sprintf("iirlp2mb : Argument B must be a vector containing numerator polynomial of the prototype low pass filter\n %s",usage)); + end + end + if ~isvector(A) then + if ~isscalar(A)then // if Denominator is having only one coeffcient then it can be passed directly + error(sprintf("iirlp2mb : A must be a vector containing denominator polynomial of the prototype low pass filter\n %s",usage)); + end + end + if(nargin == 5) + Pass = varargin(5); + switch(Pass) + case 'pass' + pass_stop = -1; + case 'stop' + pass_stop = 1; + otherwise + error(sprintf("Pass must be pass or stop\n%s",usage)) + end + else + pass_stop = -1; // Pass == 'pass' is the default + end + if(abs(Wo) <= 0) + error(sprintf("Wo is %f <= 0\n%s",abs(Wo),usage)); + end + if(abs(Wo) >= 1) + error(sprintf("Wo is %f >= 1\n%s",abs(Wo),usage)); + end + oWt = 0; + for i = 1 : length(Wt) + if(abs(Wt(i)) <= 0) + error(sprintf("Wt(%d) is %f <= 0\n%s",i,abs(Wt(i)),usage)); + end + if(abs(Wt(i)) >= 1) + error(sprintf("Wt(%d) is %f >= 1\n%s",i,abs(Wt(i)),usage)); + end + if(abs(Wt(i)) <= oWt) + error(sprintf("Wt(%d) = %f, not monotonically increasing\n%s",i,abs(Wt(i)),usage)); + else + oWt = Wt(i); + end + end + + // B(z) + // Inputs B,A specify the low pass IIR prototype filter G(z) = ---- . + // A(z) + // This module transforms G(z) into a multiband filter using the iterative + // algorithm from: + // [FFM] G. Feyh, J. Franchitti, and C. Mullis, "All-Pass Filter + // Interpolation and Frequency Transformation Problem", Proceedings 20th + // Asilomar Conference on Signals, Systems and Computers, Nov. 1986, pp. + // 164-168, IEEE. + // [FFM] moves the prototype filter position at normalized angular frequency + // .5*pi to the places specified in the Wt vector times pi. In this module, + // a generalization allows the position to be moved on the prototype filter + // to be specified as Wo*pi instead of being fixed at .5*pi. This is + // implemented using two successive allpass transformations. + // KK(z) + // In the first stage, find allpass J(z) = ---- such that + // K(z) + // jWo*pi -j.5*pi + // J(e ) = e (low pass to low pass transformation) + // + // PP(z) + // In the second stage, find allpass H(z) = ---- such that + // P(z) + // jWt(k)*pi -j(2k - 1)*.5*pi + // H(e ) = e (low pass to multiband transformation) + // + // ^ + // The variable PP used here corresponds to P in [FFM]. + // len = length(P(z)) == length(PP(z)), the number of polynomial coefficients + // + // len 1-i len 1-i + // P(z) = SUM P(i)z ; PP(z) = SUM PP(i)z ; PP(i) == P(len + 1 - i) + // i=1 i=1 (allpass condition) + // Note: (len - 1) == n in [FFM] eq. 3 + // + // The first stage computes the denominator of an allpass for translating + // from a prototype with position .5 to one with a position of Wo. It has the + // form: + // -1 + // K(2) - z + // ----------- + // -1 + // 1 - K(2)z + // + // From the low pass to low pass transformation in Table 7.1 p. 529 of A. + // Oppenheim and R. Schafer, Discrete-Time Signal Processing 3rd edition, + // Prentice Hall 2010, one can see that the denominator of an allpass for + // going in the opposite direction can be obtained by a sign reversal of the + // second coefficient, K(2), of the vector K (the index 2 not to be confused + // with a value of z, which is implicit). + + // The first stage allpass denominator computation + K = apd([%pi * Wo]); + + // The second stage allpass computation + phi = %pi * Wt; // vector of normalized angular frequencies between 0 and pi + P = apd(phi); // calculate denominator of allpass for this target vector + PP = flipdim(P,2); // numerator of allpass has reversed coefficients of P + + // The total allpass filter from the two consecutive stages can be written as + // PP + // K(2) - --- + // P P + // ----------- * --- + // PP P + // 1 - K(2)--- + // P + AllpassDen = P - (K(2) * PP); + AllpassDen = AllpassDen / AllpassDen(1); // normalize + AllpassNum = pass_stop * flipdim(AllpassDen,2); + [Num,Den] = transform(B,A,AllpassNum,AllpassDen,pass_stop); + endfunction +function [Num,Den] = transform(B,A,PP,P,pass_stop) + + // Given G(Z) = B(Z)/A(Z) and allpass H(z) = PP(z)/P(z), compute G(H(z)) + // For Pass = 'pass', transformed filter is: + // 2 nb-1 + // B1 + B2(PP/P) + B3(PP/P)^ + ... + Bnb(PP/P)^ + // ------------------------------------------------- + // 2 na-1 + // A1 + A2(PP/P) + A3(PP/P)^ + ... + Ana(PP/P)^ + // For Pass = 'stop', use powers of (-PP/P) + // + na = length(A); // the number of coefficients in A + nb = length(B); // the number of coefficients in B + // common low pass iir filters have na == nb but in general might not + n = max(na,nb); // the greater of the number of coefficients + // n-1 + // Multiply top and bottom by P^ yields: + // + // n-1 n-2 2 n-3 nb-1 n-nb + // B1(P^ ) + B2(PP)(P^ ) + B3(PP^ )(P^ ) + ... + Bnb(PP^ )(P^ ) + // --------------------------------------------------------------------- + // n-1 n-2 2 n-3 na-1 n-na + // A1(P^ ) + A2(PP)(P^ ) + A3(PP^ )(P^ ) + ... + Ana(PP^ )(P^ ) + + // Compute and store powers of P as a matrix of coefficients because we will + // need to use them in descending power order + global Ppower; // to hold coefficients of powers of P, access inside ppower() + np = length(P); + powcols = np + (np-1)*(n-2); // number of coefficients in P^(n-1) + // initialize to "Not Available" with n-1 rows for powers 1 to (n-1) and + // the number of columns needed to hold coefficients for P^(n-1) + Ppower = %nan * ones(n-1,powcols); + Ptemp = P; // start with P to the 1st power + for i = 1 : n-1 // i is the power + for j = 1 : length(Ptemp) // j is the coefficient index for this power + Ppower(i,j) = Ptemp(j); + end + Ptemp = conv(Ptemp,P); // increase power of P by one + end + + // Compute numerator and denominator of transformed filter + Num = []; + Den = []; + for i = 1 : n + // n-i + // Regenerate P^ (p_pownmi) + if((n-i) == 0) + p_pownmi = [1]; + else + p_pownmi = ppower(n-i,powcols); + end + // i-1 + // Regenerate PP^ (pp_powim1) + if(i == 1) + pp_powim1 = [1]; + else + pp_powim1 = flipdim(ppower(i-1,powcols),2); + end + if(i <= nb) + Bterm = (pass_stop^(i-1))*B(i)*conv(pp_powim1,p_pownmi); + Num = polysum(Num,Bterm); + end + if(i <= na) + Aterm = (pass_stop^(i-1))*A(i)*conv(pp_powim1,p_pownmi); + Den = polysum(Den,Aterm); + end + end + // Scale both numerator and denominator to have Den(1) = 1 + temp = Den(1); + for i = 1 : length(Den) + Den(i) = Den(i) / temp; + end + for i = 1 : length(Num) + Num(i) = Num(i) / temp; + end + +endfunction + +function P = apd(phi) // all pass denominator + // Given phi, a vector of normalized angular frequency transformation targets, + // return P, the denominator of an allpass H(z) + lenphi = length(phi); + Pkm1 = 1; // P0 initial condition from [FFM] eq. 22 + for k = 1:lenphi + P = pk(Pkm1, k, phi(k)); // iterate + Pkm1 = P; + end + +endfunction + +function Pk = pk(Pkm1, k, phik) // kth iteration of P(z) + + // Given Pkminus1, k, and phi(k) in radians , return Pk + // + // From [FFM] eq. 19 : k + // Pk = (z+1 )sin(phi(k)/2)Pkm1 - (-1) (z-1 )cos(phi(k)/2)PPkm1 + // Factoring out z + // -1 k -1 + // = z((1+z )sin(phi(k)/2)Pkm1 - (-1) (1-z )cos(phi(k)/2)PPkm1) + // PPk can also have z factored out. In H=PP/P, z in PPk will cancel z in Pk, + // so just leave out. Use + // -1 k -1 + // PK = (1+z )sin(phi(k)/2)Pkm1 - (-1) (1-z )cos(phi(k)/2)PPkm1 + // (expand) k + // = sin(phi(k)/2)Pkm1 - (-1) cos(phi(k)/2)PPkm1 + // + // -1 k -1 + // + z sin(phi(k)/2)Pkm1 + (-1) z cos(phi(k)/2)PPkm1 + Pk = zeros(1,k+1); // there are k+1 coefficients in Pk + sin_k = sin(phik/2); + cos_k = cos(phik/2); + for i = 1 : k + Pk(i) = Pk(i)+ sin_k * Pkm1(i) - ((-1)^k * cos_k * Pkm1(k+1-i)); + // + // -1 + // Multiplication by z just shifts by one coefficient + Pk(i+1) = Pk(i+1)+ sin_k * Pkm1(i) + ((-1)^k * cos_k * Pkm1(k+1-i)); + end + // now normalize to Pk(1) = 1 (again will cancel with same factor in PPk) + Pk1 = Pk(1); + for i = 1 : k+1 + Pk(i) = Pk(i) / Pk1; + end + +endfunction +function p = ppower(i,powcols) // Regenerate ith power of P from stored PPower + + global Ppower + if(i == 0) + p = 1; + else + p = []; + for j = 1 : powcols + if(isnan(Ppower(i,j))) + break; + end + p = [p, Ppower(i,j)]; + end + end + +endfunction + +function poly = polysum(p1,p2) // add polynomials of possibly different length + + n1 = length(p1); + n2 = length(p2); + if(n1 > n2) + // pad p2 + p2 = [p2, zeros(1,n1-n2)]; + elseif(n2 > n1) + // pad p1 + p1 = [p1, zeros(1,n2-n1)]; + end + poly = p1 + p2; + +endfunction + +/* +test passed +// Butterworth filter of order 3 with cutoff frequency 0.5 +B = [0.1667 0.5000 0.5000 0.1667]; +A = [1.0000e+00 -3.0531e-16 3.3333e-01 -1.8504e-17 ] ; +[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B, A, 0.5, [.2 .4 .6 .8]) // 5th argument default pass +[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B, A,[.2 .4 .6 .8],0.2) // 5th argument default pass +[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B, A, 0.5, [.2 .4 .6 .8],"stop") +[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B, A,[.2 .4 .6 .8],0.5,"stop") + +test passed +// scalar input for B and A +[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(0, 0, 0.5, [.2 .4 .6 .8]) +[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(1, 2, [.2 .4 .6 .8],0.2,"stop") + + +test passed +// complex B and A +[Num,Den,AllpassNum,AllpassDen] = iirlp2mb([%i,2+%i],[1 2*%i 3], [.2 .4 .6 .8],0.2,"stop") +*/
\ No newline at end of file diff --git a/macros/impinvar.sci b/macros/impinvar.sci index cc9aa6f..57fdb88 100644 --- a/macros/impinvar.sci +++ b/macros/impinvar.sci @@ -1,43 +1,146 @@ -function [b_out, a_out] = impinvar (b, a, fs, tol) -//This function converts analog filter with coefficients b and a to digital, conserving impulse response. -//Calling Sequence -//[b, a] = impinvar (b, a) -//[b, a] = impinvar (b, a, fs) -//[b, a] = impinvar (b, a, fs, tol) -//Parameters -//b: real or complex valued scalar or vector -//a: real or complex valued scalar or vector, order should be greater than b -//fs: real or complex value, default value 1Hz -//tol: real or complex value, default value 0.0001 -//Description -//This is an Octave function. -//This function converts analog filter with coefficients b and a to digital, conserving impulse response. -//This function does the inverse of impinvar. -//Examples -//b = 0.0081000 -//a = [2.0000000, 0.56435378, 0.4572792, 0.00705544, 0.091000] -//[ay,by] = impinvar(b,a,10) -//ay = -// 0.0000e+00 7.5293e-08 2.9902e-07 7.4238e-08 -//by = -// 1.00000 -3.96992 5.91203 -3.91428 0.97218 - -funcprot(0); -rhs = argn(2) -if(rhs<2) -error("Wrong number of input arguments.") -end - - - select(rhs) - case 2 then -// [b, a] = callOctave("impinvar",b,a) - [b_out, a_out] = callOctave("impinvar",b,a) - case 3 then -// [b, a] = callOctave("impinvar",b,a,fs) - [b_out, a_out] = callOctave("impinvar",b,a,fs) - case 4 then -// [b, a] = callOctave("impinvar",b,a,fs,tol) - [b_out, a_out] = callOctave("impinvar",b,a,fs,tol) - end +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh , SOE CUSAT +// Last Modified on : Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +/* +Calling Sequence +[b_out, a_out] = impinvar (b, a, fs, tol) +[b_out, a_out] = impinvar (b, a, fs) +[b_out, a_out] = impinvar (b, a) +Converts analog filter with coefficients b and a to digital, conserving impulse response. + +If fs is not specified, or is an empty vector, it defaults to 1Hz. + +If tol is not specified, it defaults to 0.0001 (0.1%) This function does the inverse of impinvar so that the following example should restore the original values of a and b. + +invimpinvar implements the reverse of this function. + +[b, a] = impinvar (b, a); +[b, a] = invimpinvar (b, a); +Reference: Thomas J. Cavicchi (1996) “Impulse invariance and multiple-order poles”. IEEE transactions on signal processing, Vol 44 (9): 2344–2347 +Dependencies + residue + + +*/ +function [b_out, a_out] = impinvar (b_in, a_in, fs , tol) + error("impinvar: Missing functionality .This function is not implemented yet. Will available in next release"); +endfunction +/* + if (nargin <2) + error ("impinvar: Insufficient input arguments"); + end + if nargin < 3 + fs = 1; + end + if nargin < 4 + tol = 0.0001; + end + // to be compatible with the matlab implementation where an empty vector can + // be used to get the default + if (isempty(fs)) + ts = 1; + else + ts = 1/fs; // we should be using sampling frequencies to be compatible with Matlab + end + + [r_in, p_in, k_in] = residue(b_in, a_in); // partial fraction expansion + + n = length(r_in); // Number of poles/residues + + if (length(k_in)>0) // Greater than zero means we cannot do impulse invariance + error("Order numerator >= order denominator"); + end + + r_out = zeros(1,n); // Residues of H(z) + p_out = zeros(1,n); // Poles of H(z) + k_out = 0; // Constant term of H(z) + + i=1; + while (i<=n) + m = 1; + first_pole = p_in(i); // Pole in the s-domain + while (i<n && abs(first_pole-p_in(i+1))<tol) // Multiple poles at p(i) + i=i+1; // Next residue + m=m+1; // Next multiplicity + end + [r, p, k] = z_res(r_in(i-m+1:i), first_pole, ts); // Find z-domain residues + k_out = k_out + k; // Add direct term to output + p_out(i-m+1:i) = p; // Copy z-domain pole(s) to output + r_out(i-m+1:i) = r; // Copy z-domain residue(s) to output + + i=i+1; // Next s-domain residue/pole + end + + [b_out, a_out] = inv_residue(r_out, p_out, k_out, tol); + a_out = to_real(a_out); // Get rid of spurious imaginary part + b_out = to_real(b_out); + + // Shift results right to account for calculating in z instead of z^-1 + b_out($)=[]; endfunction +*/ +// Convert residue vector for single and multiple poles in s-domain (located at sm) to +// residue vector in z-domain. The variable k is the direct term of the result. + +function [r_out, p_out, k_out] = z_res (r_in, sm, ts) + + p_out = exp(ts * sm); // z-domain pole + n = length(r_in); // Multiplicity of the pole + r_out = zeros(1,n); // Residue vector + + // First pole (no multiplicity) + k_out = r_in(1) * ts; // PFE of z/(z-p) = p/(z-p)+1; direct part + r_out(1) = r_in(1) * ts * p_out; // pole part of PFE + + for i=(2:n) // Go through s-domain residues for multiple pole + r_out(1:i) = r_out(1:i) + r_in(i) * polyrev(h1_z_deriv(i-1, p_out, ts)); // Add z-domain residues + end + +endfunction +function p_out = polyrev (p_in) + + p_out = p_in($:-1:1); + + +endfunction +function p_out = to_real(p_in) + + p_out = abs(p_in) .* sign(real(p_in)); + + endfunction +/* +//test passed +[b_out,a_out]=impinvar([1],[1 1],100); +assert_checkalmostequal(b_out,0.01,%eps,1e-4) +assert_checkalmostequal(a_out,[1 -0.99],%eps,1e-4) + +//test passed +[b_out,a_out]=impinvar([1],[1 2 1],100) +assert_checkalmostequal(b_out,[0 9.9005e-5],%eps,1e-4) +assert_checkalmostequal(a_out,[1 -1.9801 0.9802],%eps,1e-4) + +[b_out, a_out] = impinvar([1], [1 3 3 1], 100) // test passed + +[b_out, a_out] = impinvar([1 1], [1 3 3 1], 100) // test passed + +[b_out, a_out] = impinvar([1 1 1], [1 3 3 1], 100) // test passed + +// FIXME : builtin filter doesn't accepts complex parameters +// These test cases will through errors +// [b_out, a_out] = impinvar([1], [1 0 1], 100) +// [b_out, a_out] = impinvar([1 1], [1 0 1], 100) +// [b_out, a_out] = impinvar([1], [1 0 2 0 1], 100) +// [b_out, a_out] = impinvar([1 1], [1 0 2 0 1], 100) +// [b_out, a_out] = impinvar([1 1 1], [1 0 2 0 1], 100) +// [b_out, a_out] = impinvar([1 1 1 1], [1 0 2 0 1], 100) + +*/ diff --git a/macros/impz.sci b/macros/impz.sci index 1172b81..e4e1ba4 100644 --- a/macros/impz.sci +++ b/macros/impz.sci @@ -1,69 +1,125 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +/* +Calling Sequence + [x, t] = impz (b) ¶ + [x, t] = impz (b, a) ¶ + [x, t] = impz (b, a, n) ¶ + [x, t] = impz (b, a, n, fs) ¶ + impz (…) ¶ +Generate impulse-response characteristics of the filter. +The filter coefficients correspond to the the z-plane rational function with numerator b and denominator a. If a is not specified, it defaults to 1. +If n is not specified, or specified as [], it will be chosen such that the signal has a chance to die down to -120dB, or to not explode beyond 120dB, or to show five periods if there is no significant damping. +If no return arguments are requested, plot the results. +Dependencies + fftfilt +*/ function [x_r, t_r] = impz(b, a, n, fs) -// It gives Impulse response of digital filter -//Calling Sequence -//x_r = impz(b) -//x_r = impz(b, a) -//x_r = impz(b, a, n) -//x_r = impz(b, a, n, fs) -//[x_r, t_r] = impz(b, a, n, fs) + if nargin < 1 || nargin > 4 then + error(" impz : Incorrect number of input arguments ") + end + if nargin < 2 then a = [1]; end + if nargin < 3 then n = [] ; end + if nargin < 4 then fs = 1 ; end + + if (~isvector(b) && ~isscalar(b)) || (~isvector(a) && ~isscalar(a) ) then + error("impz: B and A must be vectors"); + end + if isempty(n) && length(a) > 1 then + precision = 1e-6; + r = roots(a); + maxpole = max(abs(r)); + if (maxpole > 1+precision) // unstable -- cutoff at 120 dB + n = floor(6/log10(maxpole)); + elseif (maxpole < 1-precision) // stable -- cutoff at -120 dB + n = floor(-6/log10(maxpole)); + else // periodic -- cutoff after 5 cycles + n = 30; -//Parameters -//x_r: impz chooses the number of samples and returns the response in the column vector, x_r. -//t_r : impz returns the sample times in the column vector, t_r -// b : numerator coefficients of the filter -// a : denominator coefficients of the filter -// n : samples of the impulse response t(by default ,n = length(t) and is computed automatically. -// fs : sampling frequency + // find longest period less than infinity + // cutoff after 5 cycles (w=10*pi) + rperiodic = r(find(abs(r)>=1-precision & abs(arg(r))>0)); + if ~isempty(rperiodic) + n_periodic = ceil(10*pi./min(abs(arg(rperiodic)))); + if (n_periodic > n) + n = n_periodic; + end + end -//Description -//[x_r,t_r] = impz(b,a) returns the impulse response of the filter with numerator coefficients, b, and denominator coefficients, a. impz chooses the number of samples and returns the response in the column vector, x_r, and the sample times in the column vector, t_r. t_r = [0:n-1]' and n = length(t) is computed automatically. + // find most damped pole + // cutoff at -60 dB + rdamped = r(find(abs(r)<1-precision)); + if ~isempty(rdamped) + n_damped = floor(-3/log10(max(abs(rdamped)))); + if (n_damped > n) + n = n_damped; + end + end + end + n = n + length(b); + elseif isempty(n) + n = length(b); + elseif ( ~isscalar(n)) + // TODO: missing option of having N as a vector of values to + // compute the impulse response. + error ("impz: N must be empty or a scalar"); + end -//Examples -//[x_r,t_r]=impz([0 1 1],[1 -3 3 -1],10) -//OUTPUT : -// t_r = 0. 1. 2. 3. 4. 5. 6. 7. 8. 9 -// x_r= 0. 1. 4. 9. 16. 25. 36. 49.....64......81 + if length(a) == 1 then + x = fftfilt(b/a, [1, zeros(1,n-1)]'); + else + x = filter(b, a, [1, zeros(1,n-1)]'); + end + t = [0:n-1]/fs; -//[x_r,t_r]=impz(1,[1 1],5) -//OUTPUT -// t_r = 0. 1. 2. 3. 4 -//x_r = 1. - 1. 1. - 1. 1. + if nargout >= 1 x_r = x; end + if nargout >= 2 t_r = t; end + //if nargout ~= 2 then //FIXME: fix nargout to detect 0 output arguments . till then plot it always + title("Impulse Response"); + if (fs > 1000) + t = t * 1000; + xlabel("Time (msec)"); + else + xlabel("Time (sec)"); + end + plot(t, x,); + //end -//This function is being called from Octave +endfunction +/* +test case 1 +assert_checkequal(size(impz (1, [1 -1 0.9], 100)), [100 1]) // passed -funcprot(0); -rhs = argn(2) -lhs = argn(1) -if(rhs<1 | rhs>4) -error("Wrong number of input arguments.") -end +test case 2 +// 7th order butterworth filter with fc = 0.5 //passed +B=[0.016565 0.115957 0.347871 0.579785 0.579785 0.347871 0.115957 0.016565]; +A=[1.0000e+00 -5.6205e-16 9.1997e-01 -3.6350e-16 1.9270e-01 -4.3812e-17 7.6835e-03 -4.2652e-19]; +impz(B, A) - select(rhs) - case 1 then - if(lhs==1) - [x_r] = callOctave("impz",b) - elseif(lhs==2) - [x_r,t_r] = callOctave("impz",b) - end - case 2 then - if(lhs==1) - [x_r] = callOctave("impz",b,a) - elseif(lhs==2) - [x_r,t_r] = callOctave("impz",b,a) - end - case 3 then - if(lhs==1) - [x_r] = callOctave("impz",b,a,n) - elseif(lhs==2) - [x_r,t_r] = callOctave("impz",b,a,n) - end - case 4 then - if(lhs==1) - [x_r] = callOctave("impz",b,a,n,fs) - elseif(lhs==2) - [x_r,t_r] = callOctave("impz",b,a,n,fs) - end - end -endfunction +test case 3 +// +[x_r,tr]=impz([0.4 0.2 8 2] ,1,10) +assert_checkalmostequal(x_r,[ 0.4000000 0.2000000 8. 2.0000000 1.943D-16 1.388D-17 0. 0 0. 0.]',%eps,1e-4); +assert_checkalmostequal(tr,0:9,%eps,1e-4); + +//test case 4 +[xr,tr]=impz([0.4 0.2],[4 5 6],3,10) + + +// test case 5 +B = [0.021895 0.109474 0.218948 0.218948 0.109474 0.021895]; +A = [1.0000 -1.2210 1.7567 -1.3348 0.7556 -0.2560]; +impz(B, A) + + */
\ No newline at end of file diff --git a/macros/inv_residue.sci b/macros/inv_residue.sci new file mode 100644 index 0000000..e7354a6 --- /dev/null +++ b/macros/inv_residue.sci @@ -0,0 +1,54 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +// Inverse of residue function +function [b_out, a_out] = inv_residue(r_in, p_in, k_in, tol) + + n = length(r_in); // Number of poles/residues + + k = 0; // Capture constant term + if (length(k_in)==1) // A single direct term (order N = order D) + k = k_in(1); // Capture constant term + elseif (length(k_in)>1) // Greater than one means non-physical system + error("Order numerator > order denominator"); + end + + a_out = octave_poly(p_in); + + b_out = zeros(1,n+1); + b_out = b_out + k*a_out; // Constant term: add k times denominator to numerator + i=1; + while (i<=n) + term = [1 -p_in(i)]; // Term to be factored out + p = r_in(i)*deconv(a_out,term); // Residue times resulting polynomial + p = prepad(p, n+1, 0, 2); // Pad for proper length + b_out = b_out + p; + + m = 1; + mterm = term; + first_pole = p_in(i); + while (i<n && abs(first_pole-p_in(i+1))<tol) // Multiple poles at p(i) + i=i+1; // Next residue + m=m+1; + mterm = conv(mterm, term); // Next multiplicity to be factored out + p = r_in(i) * deconv(a_out, mterm); // Resulting polynomial + p = prepad(p, n+1, 0, 2); // Pad for proper length + b_out = b_out + p; + end + i=i+1; + end + +endfunction +function ocpoly = octave_poly(A) + p = poly(A, 'x'); + c = coeff(p); + ocpoly = c($:-1:1); +endfunction
\ No newline at end of file diff --git a/macros/invfreq.sci b/macros/invfreq.sci index a720115..0fe4c81 100644 --- a/macros/invfreq.sci +++ b/macros/invfreq.sci @@ -1,43 +1,304 @@ -function [B,A] = invfreq(H,F,nB,nA,W,iter,tol, plane) -// Calculates inverse frequency vectors -// -// Calling Sequence -//[B,A] = invfreq(H,F,nB,nA) -//[B,A] = invfreq(H,F,nB,nA,W) -//[B,A] = invfreq(H,F,nB,nA,W,[],[],plane) -//[B,A] = invfreq(H,F,nB,nA,W,iter,tol,plane) -// -// Parameters -// H: desired complex frequency response,It is assumed that A and B are real polynomials, hence H is one-sided. -// F: vector of frequency samples in radians -// nA: order of denominator polynomial A -// nB: order of numerator polynomial B -// -// Description -//Fit filter B(z)/A(z) or B(s)/A(s) to complex frequency response at frequency points F. A and B are real polynomial coefficients of order nA and nB respectively. Optionally, the fit-errors can be weighted vs frequency according to the weights W. Also, the transform plane can be specified as either 's' for continuous time or 'z' for discrete time. 'z' is chosen by default. Eventually, Steiglitz-McBride iterations will be specified by iter and tol. -// -// Examples -// [B,A] = butter(12,1/4); -// [H,w] = freqz(B,A,128); -// [Bh,Ah] = invfreq(H,F,4,4); -// Hh = freqz(Bh,Ah); -// disp(sprintf('||frequency response error|| = %f',norm(H-Hh))); -// -funcprot(0); -lhs= argn(1); -rhs= argn(2); -if(rhs < 4 | rhs > 8 | rhs == 6 | rhs == 7 ) -error("Wrong number of input arguments"); -end +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +/* + : [B,A] = invfreq(H,F,nB,nA,W) + : [B,A] = invfreq(H,F,nB,nA,W,[],[],plane) + : [B,A] = invfreq(H,F,nB,nA,W,iter,tol,plane) + + Fit filter B(z)/A(z) or B(s)/A(s) to complex frequency response at frequency points F. + + A and B are real polynomial coefficients of order nA and nB respectively. Optionally, the fit-errors can be weighted vs frequency according to the weights W. Also, the transform plane can be specified as either ’s’ for continuous time or ’z’ for discrete time. ’z’ is chosen by default. Eventually, Steiglitz-McBride iterations will be specified by iter and tol. + + H: desired complex frequency response It is assumed that A and B are real polynomials, hence H is one-sided. + + F: vector of frequency samples in radians + + nA: order of denominator polynomial A + + nB: order of numerator polynomial B + + plane=’z’: F on unit circle (discrete-time spectra, z-plane design) + + plane=’s’: F on jw axis (continuous-time spectra, s-plane design) + + H(k) = spectral samples of filter frequency response at points zk, where zk=exp(sqrt(-1)*F(k)) when plane=’z’ (F(k) in [0,.5]) and zk=(sqrt(-1)*F(k)) when plane=’s’ (F(k) nonnegative) +*/ +// FIXME: implement Steiglitz-McBride iterations +// FIXME: improve numerical stability for high order filters (matlab is a bit better) +// FIXME: modify to accept more argument configurations +function [B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, plane,varargin) + + if nargin < 4 then + error("invfreq : Incorrect number of input arguments ") + end -select(rhs) - case 4 then - [B,A]= callOctave("invfreq", H,F,nB,nA); - case 5 then - [B,A]= callOctave("invfreq", H,F,nB,nA, W); - case 8 then - [B,A]= callOctave("invfreq", H,F,nB, nA,iter, tol, plane); + if ~isvector(H) && ~isscalar(H) then + error("invfreq : H is the desired frequency response , a vector expected") + end + + if ~isvector(F) && ~isscalar(F) then + error("invfreq : F is a vector of frequency samples in radians") + end + + if max(size(nB)) > 1 then zB = nB(2); nB = nB(1); else zB = 0; end + n = max(nA, nB); + m = n+1; mA = nA+1; mB = nB+1; + nF = max(size(F)); + if nargin < 5 || isempty(W) then W = ones(1, nF); end + if nargin < 6 then iter = []; end + if nargin < 7 then tol = []; end + if nargin < 8 || isempty(tr) then tr = ''; end + if nargin < 9 then plane = 'z'; end + if nargin < 10 then varargin = {}; end + if ( strcmp (plane, "s") && strcmp (plane, "z")) + error (sprintf("invfreq: invealid PLANE argument %s, expected s or z ", plane)) + end + + fname = ["invfreq", plane]; + + if (nF ~= max(size(H))) then + error ("%s: Length of H and F must be the same\n", fname) + end + + if (~ isempty (iter) || ~ isempty (tol)) then + warning (sprintf("%s: iterative algorithm not yet implemented, ", ... + "ITER and TOL arguments are ignored\n", fname)); + end + +////////////////////////////////////////////////////////////// +norm = %f ; // should we normalize freqs to avoid matrices with rank deficiency ? +method = 'LS'; // by default, use Ordinary Least Square to solve normal equations +prop = varargin; +if length(prop) > 0 then + indi = 1; + while indi < length(prop) + switch prop(indi) + case 'norm' + if indi < length(prop) && ~(type(prop(indi+1)) == 10) + norm = prop(indi+1); + indi = indi + 2; // Skip the processed element + else + norm = %t; // Default true + indi = indi + 1; + end + + case 'method' + if indi < length(prop) && type(prop(indi+1)) == 10 && strcmp(prop(indi+1), "norm") + method = prop(indi+1); + indi = indi + 2; // Skip the processed element + else + error("invfreq : incorrect/missing method argument"); + indi = indi + 1; + end + + otherwise + disp("Ignoring unknown or incomplete argument"); + indi = indi + 1; + end + end end + + +//////////////////////////////////////////////////////////////// + + + Ruu = zeros(mB, mB); Ryy = zeros(nA, nA); Ryu = zeros(nA, mB); + Pu = zeros(mB, 1); Py = zeros(nA,1); + if ~strcmp(tr,'trace') + disp(' ') + disp('Computing nonuniformly sampled, equation-error, rational filter.'); + disp(['plane = ',plane]); + disp(' ') + end + + s = sqrt(-1)*F; + switch plane + case 'z' + if max(F) > %pi || min(F) < 0 then + disp('hey, you frequency is outside the range 0 to %pi, making my own') + F = linspace(0, %pi, max(size(H))); + s = sqrt(-1)*F; + end + s = exp(-s); + case 's' + if max(F) > 1e6 && n > 5 then + if ~norm then + disp('Be careful, there are risks of generating singular matrices'); + disp('Call invfreqs as (..., norm, 1) to avoid it'); + else + Fmax = max(F); s = sqrt(-1)*F/Fmax; + end + end + end + + ////////////////////////////// + ///////////////////////////// + for k=1:nF, + Zk = (s(k).^[0:n]).'; + Hk = H(k); + aHks = Hk*conj(Hk); + Rk = (W(k)*Zk)*Zk'; + rRk = real(Rk); + Ruu = clean(Ruu + rRk(1:mB, 1:mB)); + Ryy = Ryy + aHks*rRk(2:mA, 2:mA); + Ryu = Ryu + real(Hk*Rk(2:mA, 1:mB)); + Pu = Pu + W(k)*real(conj(Hk)*Zk(1:mB)); + Py = Py + (W(k)*aHks)*real(Zk(2:mA)); + end + Rr = ones(max(size(s)), mB+nA); Zk = s; + for k = 1:min(nA, nB), + Rr(:, 1+k) = Zk; + Rr(:, mB+k) = -Zk.*H; + Zk = Zk.*s; + end + for k = 1+min(nA, nB):max(nA, nB)-1, + if k <= nB, Rr(:, 1+k) = Zk; end + if k <= nA, Rr(:, mB+k) = -Zk.*H; end + Zk = Zk.*s; + end + k = k+1; + if k <= nB then Rr(:, 1+k) = Zk; end + if k <= nA then Rr(:, mB+k) = -Zk.*H; end + + // complex to real equation system -- this ensures real solution + Rr = Rr(:, 1+zB:$); + Rr = [real(Rr); imag(Rr)]; Pr = [real(H(:)); imag(H(:))]; + // normal equations -- keep for ref + // Rn= [Ruu(1+zB:mB, 1+zB:mB), -Ryu(:, 1+zB:mB)'; -Ryu(:, 1+zB:mB), Ryy]; + // Pn= [Pu(1+zB:mB); -Py]; + //////////////////////////////////////////////// + switch method + case {'ls' 'LS'} + // avoid scaling errors with Theta = R\P; + // [Q, R] = qr([Rn Pn]); Theta = R(1:$, 1:$-1)\R(1:$, $); + [Q, R] = qr([Rr Pr]); Theta = pinv(R(1:$-1, 1:$-1)) * R(1:$-1, $); + ////////////////////////////////////////////////// + // SigN = R($, $-1); + SigN = R($, $); + case {'tls' 'TLS'} + // [U, S, V] = svd([Rn Pn]); + // SigN = S($, $-1); + // Theta = -V(1:$-1, $)/V($, $); + [U, S, V] = svd([Rr Pr], 0); + SigN = S($, $); + Theta = -V(1:$-1, $)/V($, $); + case {'mls' 'MLS' 'qr' 'QR'} + // [Q, R] = qr([Rn Pn], 0); + // solve the noised part -- DO NOT USE ECONOMY SIZE ~ + // [U, S, V] = svd(R(nA+1:$, nA+1:$)); + // SigN = S($, $-1); + // Theta = -V(1:$-1, $)/V($, $); + // unnoised part -- remove B contribution and back-substitute + // Theta = [R(1:nA, 1:nA)\(R(1:nA, $) - R(1:nA, nA+1:$-1)*Theta) + // Theta]; + // solve the noised part -- economy size OK as #rows > #columns + [Q, R] = qr([Rr Pr], 0); + eB = mB-zB; sA = eB+1; + [U, S, V] = svd(R(sA:$, sA:$)); + // noised (A) coefficients + Theta = -V(1:$-1, $)/V($, $); + // unnoised (B) part -- remove A contribution and back-substitute + Theta = [R(1:eB, 1:eB)\(R(1:eB, $) - R(1:eB, sA:$-1)*Theta) + Theta]; + SigN = S($, $); + otherwise + error(sprintf("invfreq : unknown method %s", method)); + end + + B = [zeros(zB, 1); Theta(1:mB-zB)].'; + A = [1; Theta(mB-zB+(1:nA))].'; + if ~strcmp(plane,'s') + B = B(mB:-1:1); + A = A(mA:-1:1); + if norm, // Frequencies were normalized -- unscale coefficients + Zk = Fmax.^[n:-1:0].'; + for k = nB:-1:1+zB, B(k) = B(k)/Zk(k); end + for k = nA:-1:1, A(k) = A(k)/Zk(k); end + end + end endfunction - +/* +// method - LS +test case 1 // passed + +[B,A,Sign] = invfreq(1,1,1,1,1,[],[],'','z','norm',1,'method','LS') +assert_checkequal(B,[0.6314 0.3411]) +assert_checkequal(A,[1 -0.3411]) +assert_checkequal(Sign,0) + +[B,A,Sign] = invfreq(1,1,1,1,1,[],[],'','s') +assert_checkequal(B,[0 1]) +assert_checkequal(A,[0 1]) +assert_checkequal(Sign,0) + + +test case 2 // passed +order = 6 +fc = 1/2 +n = 128 +B = [ 0.029588 0.177529 0.443823 0.591764 0.443823 0.177529 0.029588] ; +A = [ 1.0000e+00 -6.6613e-16 7.7770e-01 -2.8192e-16 1.1420e-01 -1.4472e-17 1.7509e-03]; +[H,w] = freqz(B,A,n) ; +[Bh , Ah] = invfreq(H,w,order,order); +[Hh,wh] = freqz(Bh,Ah,n); +plot(w,[abs(H), abs(Hh)]) +xlabel("Frequency (rad/sample)"); +ylabel("Magnitude"); +legend('Original','Measured'); +err = norm(H-Hh); +disp(sprintf('L2 norm of frequency response error = %f',err)); + +test case 3 // passed +// buttter worth filter of order 12 and fc=1/4 +B = [ 1.1318e-06 1.3582e-05 7.4702e-05 2.4901e-04 5.6026e-04 8.9642e-04 1.0458e-03 8.9642e-04 5.6026e-04 2.4901e-04 7.4702e-05 1.3582e-05 1.1318e-06]; +A = [ 1.0000e+00 -5.9891e+00 1.7337e+01 -3.1687e+01 4.0439e+01 -3.7776e+01 2.6390e+01 -1.3851e+01 5.4089e+00 -1.5296e+00 2.9688e-01 -3.5459e-02 1.9688e-03]; +[H,w] = freqz(B,A,128); +[Bh,Ah] = invfreq(H,w,4,4); +[Hh,wh] = freqz(Bh,Ah,128); +disp(sprintf('||frequency response error||= %f',norm(H-Hh))); + + +method TLS + +test case 1 // passed + +B = [ 1.1318e-06 1.3582e-05 7.4702e-05 2.4901e-04 5.6026e-04 8.9642e-04 1.0458e-03 8.9642e-04 5.6026e-04 2.4901e-04 7.4702e-05 1.3582e-05 1.1318e-06]; +A = [ 1.0000e+00 -5.9891e+00 1.7337e+01 -3.1687e+01 4.0439e+01 -3.7776e+01 2.6390e+01 -1.3851e+01 5.4089e+00 -1.5296e+00 2.9688e-01 -3.5459e-02 1.9688e-03]; +[H,w] = freqz(B,A,128); +[Bh,Ah] = invfreq(H,w,4,4,[],[],[],'','z','norm',1,'method','TLS'); +[Hh,wh] = freqz(Bh,Ah,128); +disp(sprintf('||frequency response error||= %f',norm(H-Hh))); + + + + +method MLS - passed + + +// elliptic filter with ellip (5, 1, 90, [.1 .2]) +n = 128 +B = [ 1.3214e-04 -6.6404e-04 1.4928e-03 -1.9628e-03 1.4428e-03 0 -1.4428e-03 1.9628e-03 -1.4928e-03 6.6404e-04 -1.3214e-04] ; +A = [ 1.0000 -8.6483 34.6032 -84.2155 137.9276 -158.7598 130.0425 -74.8636 29.0044 -6.8359 0.7456]; +[H,w] = freqz(B,A,n) ; + +[Bh,Ah] = invfreq(H,w,4,4,[],[],[],'','z','norm',1,'method','MLS'); +[Hh,wh] = freqz(Bh,Ah,n); +plot(w,[abs(H), abs(Hh)]) +xlabel("Frequency (rad/sample)"); +ylabel("Magnitude"); +legend('Original','Measured'); +err = norm(H-Hh); +disp(sprintf('L2 norm of frequency response error = %f',err)); + + +*/ diff --git a/macros/invfreqs.sci b/macros/invfreqs.sci index 32e732b..99c721b 100644 --- a/macros/invfreqs.sci +++ b/macros/invfreqs.sci @@ -1,96 +1,110 @@ -function [B,A,C] = invfreqs(H,F,nB,nA,W,iter,tol,trace) -//Fit filter B(s)/A(s)to the complex frequency response H at frequency points F. A and B are real polynomial coefficients of order nA and nB. -//Calling Sequence -//[B,A,C] = invfreqs(H,F,nB,nA,W,iter,tol,trace) -//[B,A,C] = invfreqs(H,F,nB,nA,W) -//[B,A,C] = invfreqs(H,F,nB,nA) -//Parameters -//H: desired complex frequency response. -//F: frequency (must be same length as H). -//nB: order of the numerator polynomial B. -//nA: order of the denominator polynomial A. -//W: vector of weights (must be same length as F). -//Description -//This is an Octave function. -//Fit filter B(s)/A(s)to the complex frequency response H at frequency points F. A and B are real polynomial coefficients of order nA and nB. -//Optionally, the fit-errors can be weighted vs frequency according to the weights W. -//Note: all the guts are in invfreq.m -//Examples -//B = [1/2 1]; -//A = [1 1]; -//w = linspace(0,4,128); -//H = freqs(B,A,w); -//[Bh,Ah, C] = invfreqs(H,w,1,1); -//Bh = -// -// 0.50000 1.00000 -// -//Ah = -// -// 1.0000 1.0000 -// -//C = -3.0964e-15 - -funcprot(0); -lhs = argn(1) -rhs = argn(2) -if (rhs < 4 | rhs > 8) -error("Wrong number of input arguments.") -end - -select(rhs) - - case 4 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA) - else - error("Wrong number of output argments.") - end - - case 5 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA,W) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA,W) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA,W) - else - error("Wrong number of output argments.") - end - case 6 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA,W,iter) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter) - else - error("Wrong number of output argments.") - end - case 7 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA,W,iter,tol) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol) - else - error("Wrong number of output argments.") - end - case 8 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace) - else - error("Wrong number of output argments.") - end - end +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +// FIXME: check invfreq.sci for todo's + +/* + :[B,A] = invfreqs(H,F,nB,nA) + :[B,A] = invfreqs(H,F,nB,nA,W) + :[B,A] = invfreqs(H,F,nB,nA,W,iter,tol,'trace') + + Fit filter B(s)/A(s)to the complex frequency response H at frequency points F. + + A and B are real polynomial coefficients of order nA and nB. + + Optionally, the fit-errors can be weighted vs frequency according to the weights W. + + Note: all the guts are in invfreq.m + + H: desired complex frequency response + + F: frequency (must be same length as H) + + nA: order of the denominator polynomial A + + nB: order of the numerator polynomial B + + W: vector of weights (must be same length as F) + +*/ +// Dependencies +// invfreq +function [B, A, SigN] = invfreqs(H,F,nB,nA,W,iter,tol,tr, varargin) + + if nargin < 9 + varargin = {}; + if nargin < 8 + tr = ''; + if nargin < 7 + tol = []; + if nargin < 6 + iter = []; + if nargin < 5 + W = ones(1,length(F)); + end + end + end + end + end + + // now for the real work + [B, A, SigN] = invfreq(H, F,nB, nA, W, iter, tol, tr, 's', varargin); + endfunction +/* +demo + B = [1 0 0]; + A = [1 6 15 15]/15; + w = linspace(0, 8, 128); + [H0 ,_ ] = freqz(B, A, w); + Nn = (rand(size(w,1),size(w,2),'normal')+%i*rand(size(w,1),size(w,2),'normal'))/sqrt(2); + order = length(A) - 1; + [Bh, Ah, Sig0] = invfreqs(H0, w, [length(B)-1 2], length(A)-1); + [Hh ,_ ] = freqz(Bh,Ah,w); + [BLS, ALS, SigLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "LS"); + [HLS,_] = freqz(BLS, ALS, w); + [BTLS, ATLS, SigTLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "TLS"); + [HTLS,_] = freqz(BTLS, ATLS, w); + [BMLS, AMLS, SigMLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "QR"); + [HMLS,_] = freqz(BMLS, AMLS, w); + plot(w,[abs(H0); abs(Hh)]) + xlabel("Frequency (rad/sec)"); + ylabel("Magnitude"); + legend('Original','Measured'); + err = norm(H0-Hh); + disp(sprintf('L2 norm of frequency response error = %f',err)); +*/ +/* Octave version + B = [1 0 0]; + A = [1 6 15 15]/15; + w = linspace(0, 8, 128); + [H0 ,_ ] = freqz(B, A, w); + Nn = (randn(size(w,1),size(w,2))+i*randn(size(w,1),size(w,2)))/sqrt(2); + order = length(A) - 1; + [Bh, Ah, Sig0] = invfreqs(H0, w, [length(B)-1 2], length(A)-1); + [Hh ,_ ] = freqz(Bh,Ah,w); + [BLS, ALS, SigLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "LS"); + [HLS,_] = freqz(BLS, ALS, w); + [BTLS, ATLS, SigTLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "TLS"); + [HTLS,_] = freqz(BTLS, ATLS, w); + [BMLS, AMLS, SigMLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "QR"); + [HMLS,_] = freqz(BMLS, AMLS, w); + plot(w,[abs(H0); abs(Hh)]) + xlabel("Frequency (rad/sec)"); + ylabel("Magnitude"); + legend('Original','Measured'); + err = norm(H0-Hh); + disp(sprintf('L2 norm of frequency response error = %f',err)); +*/ + + + diff --git a/macros/invfreqz.sci b/macros/invfreqz.sci index d051499..433100b 100644 --- a/macros/invfreqz.sci +++ b/macros/invfreqz.sci @@ -1,94 +1,107 @@ -function [B,A,C] = invfreqz(H,F,nB,nA,W,iter,tol,trace) -//Fit filter B(z)/A(z)to the complex frequency response H at frequency points F. A and B are real polynomial coefficients of order nA and nB. -//Calling Sequence -//[B,A,C] = invfreqz(H,F,nB,nA,W,iter,tol,trace) -//[B,A,C] = invfreqz(H,F,nB,nA,W) -//[B,A,C] = invfreqz(H,F,nB,nA) -//Parameters -//H: desired complex frequency response. -//F: frequency (must be same length as H). -//nB: order of the numerator polynomial B. -//nA: order of the denominator polynomial A. -//W: vector of weights (must be same length as F). -//Description -//This is an Octave function. -//Fit filter B(z)/A(z)to the complex frequency response H at frequency points F. A and B are real polynomial coefficients of order nA and nB. -//Optionally, the fit-errors can be weighted vs frequency according to the weights W. -//Note: all the guts are in invfreq.m -//Examples -//[B,A] = butter(4,1/4); -//[H,F] = freqz(B,A); -//[Bh,Ah,C] = invfreq(H,F,4,4) -//Bh = -// -// 0.010209 0.040838 0.061257 0.040838 0.010209 -// -//Ah = -// -// 1.00000 -1.96843 1.73586 -0.72447 0.12039 -// -//C = -7.7065e-15 +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +// FIXME: check invfreq.sci for todo's +/* + : [B,A] = invfreqz(H,F,nB,nA) ¶ + : [B,A] = invfreqz(H,F,nB,nA,W) ¶ + : [B,A] = invfreqz(H,F,nB,nA,W,iter,tol,'trace') ¶ -funcprot(0); -lhs = argn(1) -rhs = argn(2) -if (rhs < 4 | rhs > 8) -error("Wrong number of input arguments.") -end + Fit filter B(z)/A(z)to the complex frequency response H at frequency points F. -select(rhs) - - case 4 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA) - else - error("Wrong number of output argments.") - end + A and B are real polynomial coefficients of order nA and nB. Optionally, the fit-errors can be weighted vs frequency according to the weights W. - case 5 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA,W) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA,W) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA,W) - else - error("Wrong number of output argments.") - end - case 6 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA,W,iter) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter) - else - error("Wrong number of output argments.") - end - case 7 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA,W,iter,tol) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol) - else - error("Wrong number of output argments.") - end - case 8 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace) - else - error("Wrong number of output argments.") - end - end -endfunction + Note: all the guts are in invfreq.m + + H: desired complex frequency response + + F: normalized frequency (0 to pi) (must be same length as H) + + nA: order of the denominator polynomial A + + nB: order of the numerator polynomial B + W: vector of weights (must be same length as F) + +*/ +// Dependencies +// invfreq +function [B, A, SigN] = invfreqz(H, F, nB, nA, W, iter, tol, tr, varargin) + + if nargin < 9 + varargin = {}; + if nargin < 8 + tr = ''; + if nargin < 7 + tol = []; + if nargin < 6 + iter = []; + if nargin < 5 + W = ones(1,length(F)); + end + end + end + end + end + // now for the real work + [B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, 'z', varargin); + +endfunction +/* +demo +order = 9; //order of test filter +//going to 10 or above leads to numerical instabilities and large errors +fc = 1/2; // sampling rate / 4 +n = 128; // frequency grid size +// butterworth filter of order 9 and fc=0.5 +B0 = [5.1819e-03 4.6637e-02 1.8655e-01 4.3528e-01 6.5292e-01 6.5292e-01 4.3528e-01 1.8655e-01 4.6637e-02 5.1819e-03]; +A0 = [ 1.0000e+00 -8.6736e-16 1.2010e+00 -7.7041e-16 4.0850e-01 -1.7013e-16 4.2661e-02 -9.0155e-18 9.6666e-04 -5.3661e-20]; +[H0, w] = freqz(B0, A0, n); +Nn = (rand(size(w,1),size(w,2),'normal')+%i*rand(size(w,1),size(w,2),'normal'))/sqrt(2); +[Bh, Ah, Sig0] = invfreqz(H0, w, order, order); +[Hh, wh] = freqz(Bh, Ah, n); +[BLS, ALS, SigLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "LS"); +[HLS _ ] = freqz(BLS, ALS, n); +[BTLS, ATLS, SigTLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "TLS"); +[HTLS _ ]= freqz(BTLS, ATLS, n); +[BMLS, AMLS, SigMLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "QR"); +[HMLS _ ] = freqz(BMLS, AMLS, n); +plot(w,[abs(H0) abs(Hh)]) +xlabel("Frequency (rad/sample)"); +ylabel("Magnitude"); +legend('Original','Measured'); +err = norm(H0-Hh); +disp(sprintf('L2 norm of frequency response error = %f',err)); + +*/ +/* +order = 9; +fc = 1/2; +n = 128; +B0 = [5.1819e-03 4.6637e-02 1.8655e-01 4.3528e-01 6.5292e-01 6.5292e-01 4.3528e-01 1.8655e-01 4.6637e-02 5.1819e-03]; +A0 = [ 1.0000e+00 -8.6736e-16 1.2010e+00 -7.7041e-16 4.0850e-01 -1.7013e-16 4.2661e-02 -9.0155e-18 9.6666e-04 -5.3661e-20]; +[H0, w] = freqz(B0, A0, n); +Nn = (randn(size(w,1),size(w,2))+i*randn(size(w,1),size(w,2)))/sqrt(2); +[Bh, Ah, Sig0] = invfreqz(H0, w, order, order); +[Hh, wh] = freqz(Bh, Ah, n); +[BLS, ALS, SigLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "LS"); +[HLS _ ] = freqz(BLS, ALS, n); +[BTLS, ATLS, SigTLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "TLS"); +[HTLS _ ]= freqz(BTLS, ATLS, n); +[BMLS, AMLS, SigMLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "QR"); +[HMLS _ ] = freqz(BMLS, AMLS, n); +plot(w,[abs(H0) abs(Hh)]) +xlabel("Frequency (rad/sample)"); +ylabel("Magnitude"); +legend('Original','Measured'); +err = norm(H0-Hh); +disp(sprintf('L2 norm of frequency response error = %f',err)); +*/ diff --git a/macros/invimpinvar.sci b/macros/invimpinvar.sci index 4e8bd7b..c854324 100644 --- a/macros/invimpinvar.sci +++ b/macros/invimpinvar.sci @@ -1,42 +1,129 @@ -function [b_out, a_out] = invimpinvar (b, a, fs, tol) -//This function converts digital filter with coefficients b and a to analog, conserving impulse response. -//Calling Sequence -//[b, a] = impinvar (b, a) -//[b, a] = impinvar (b, a, fs) -//[b, a] = impinvar (b, a, fs, tol) -//Parameters -//b: real or complex valued scalar or vector -//a: real or complex valued scalar or vector, order should be greater than b -//fs: real or complex value, default value 1Hz -//tol: real or complex value, default value 0.0001 -//Description -//This is an Octave function. -//This function converts digital filter with coefficients b and a to analog, conserving impulse response. -//This function does the inverse of impinvar. -//Examples -//b = 0.0081000 -//a = [2.0000000, 0.56435378, 0.4572792, 0.00705544, 0.091000] -//[ay, by] = invimpinvar(b,a,10) -//ay = -// -1.6940e-16 4.6223e+00 -4.5210e+00 7.2880e+02 -//by = -// Columns 1 through 4: -// 1.0000e+00 3.0900e+01 9.6532e+02 1.2232e+04 -// Column 5: -// 1.1038e+05 -funcprot(0); -rhs = argn(2) -if(rhs<2) -error("Wrong number of input arguments.") -end - - - select(rhs) - case 2 then - [b_out,a_out] = callOctave("invimpinvar",b,a) - case 3 then - [b_out,a_out] = callOctave("invimpinvar",b,a,fs) - case 4 then - [b_out,a_out] = callOctave("invimpinvar",b,a,fs,tol) - end +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +// Impulse invariant conversion from s to z domain +/* + [b_out, a_out] = invimpinvar (b, a, fs, tol) + [b_out, a_out] = invimpinvar (b, a, fs) + [b_out, a_out] = invimpinvar (b, a) + + Converts digital filter with coefficients b and a to analog, conserving impulse response. + + This function does the inverse of impinvar so that the following example should restore the original values of a and b. + + [b, a] = impinvar (b, a); + [b, a] = invimpinvar (b, a); + + If fs is not specified, or is an empty vector, it defaults to 1Hz. + + If tol is not specified, it defaults to 0.0001 (0.1%) + Dependencies + residue + inv_residue + */ +function [b_out, a_out] = invimpinvar (b_in, a_in, fs, tol) + + error("invimpinvar: Missing functionality not implemented in this release .Will be Available soon "); endfunction +// FIXME : fix filter function first . till then drop this fun +/* + if (nargin <2) + error("invimpinvar: Insufficient input arguments"); + end + + if nargin < 3 then fs = 1; end + if nargin < 4 then tol = 0.0001; end + // to be compatible with the matlab implementation where an empty vector can + // be used to get the default + if (isempty(fs)) + ts = 1; + else + ts = 1/fs; // we should be using sampling frequencies to be compatible with Matlab + end + + b_in = [b_in 0]; // so we can calculate in z instead of z^-1 + + [r_in, p_in, k_in] = residue(b_in, a_in); // partial fraction expansion + + // clean r_in for zero values + n = length(r_in); // Number of poles/residues + + if (length(k_in) > 1) // Greater than one means we cannot do impulse invariance + error("Order numerator > order denominator"); + end + + r_out = zeros(1,n); // Residues of H(s) + sm_out = zeros(1,n); // Poles of H(s) + + i=1; + while (i<=n) + m=1; + first_pole = p_in(i); // Pole in the z-domain + while (i<n && abs(first_pole-p_in(i+1))<tol) // Multiple poles at p(i) + i=i+1; // Next residue + m=m+1; // Next multiplicity + end + [r, sm, k]= inv_z_res(r_in(i-m+1:i), first_pole, ts); // Find s-domain residues + k_in = k_in - k; // Just to check, should end up zero for physical system + sm_out(i-m+1:i) = sm; // Copy s-domain pole(s) to output + r_out(i-m+1:i) = r; // Copy s-domain residue(s) to output + + i=i+1; // Next z-domain residue/pole + end + [b_out, a_out] = inv_residue(r_out, sm_out , 0, tol); + a_out = to_real(a_out); // Get rid of spurious imaginary part + b_out = to_real(b_out); + + b_out = polyreduce(b_out); + +endfunction +*/ +// Inverse function of z_res (see impinvar source) + +function [r_out, sm_out, k_out] = inv_z_res (r_in,p_in,ts) + + n = length(r_in); // multiplicity of the pole + r_in = r_in.'; // From column vector to row vector + + j=n; + while (j>1) // Go through residues starting from highest order down + r_out(j) = r_in(j) / ((ts * p_in)^j); // Back to binomial coefficient for highest order (always 1) + r_in(1:j) = r_in(1:j) - r_out(j) * polyrev(h1_z_deriv(j-1,p_in,ts)); // Subtract highest order result, leaving r_in(j) zero + j=j-1; + end + + // Single pole (no multiplicity) + r_out(1) = r_in(1) / ((ts * p_in)); + k_out = r_in(1) / p_in; + sm_out = log(p_in) / ts; +endfunction + +/* +// tests passed +[b_out,a_out]=invimpinvar([1],[1 -0.5],0.01) +[b_out,a_out]=invimpinvar([1],[1 -1 0.25],0.01) +[b_out,a_out]=invimpinvar([1 1],[1 -1 0.25],0.01) +[b_out,a_out]=invimpinvar([1],[1 -1.5 0.75 -0.125],0.01) +[b_out,a_out]=invimpinvar([1 1],[1 -1.5 0.75 -0.125],0.01) + + + +// FIXME : built in filter doesn't support complex parameters +// Because of this thsese test cases are failing +//[b_out,a_out]=invimpinvar([1],[1 0 0.25],0.01) +// [b_out,a_out]=invimpinvar([1 1],[1 0 0.25],0.01) +// [b_out,a_out]=invimpinvar([1],[1 0 0.5 0 0.0625],0.01) +// [b_out,a_out]=invimpinvar([1 1],[1 0 0.5 0 0.0625],0.01) +// [b_out,a_out]=invimpinvar([1 1 1],[1 0 0.5 0 0.0625],0.01 +// [b_out,a_out]=invimpinvar([1 1 1 1],[1 0 0.5 0 0.0625],0.01) + +*/ diff --git a/macros/marcumq.sci b/macros/marcumq.sci index 630ec97..29f4416 100644 --- a/macros/marcumq.sci +++ b/macros/marcumq.sci @@ -1,38 +1,389 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +/* +Calling Sequence + q = marcumq (a, b) + q = marcumq (a, b, m) + q = marcumq (a, b, m, tol) + +Input and Output parameters +a — Noncentrality parameter --- nonnegative scalar | array of nonnegative numbers +b — Argument of Marcum Q-function --- nonnegative scalar | array of nonnegative numbers +m — Order of generalized Marcum Q-function --- positive integer | array of positive integers + +Compute the generalized Marcum Q function of order `m` with noncentrality parameter `a` and argument `b`. +If the order `m` is omitted, it defaults to 1. An optional relative tolerance `tol` may be included, +and the default value is `eps`. + +If the input arguments are commensurate vectors, this function will produce a table of values. + +This function computes Marcum’s Q function using the infinite Bessel series, +which is truncated when the relative error is less than the specified tolerance. +The accuracy is limited by that of the Bessel functions, so reducing the tolerance is probably not useful. + +References: +- Marcum, "Tables of Q Functions", Rand Corporation. +- R.T. Short, "Computation of Noncentral Chi-squared and Rice Random Variables", + www.phaselockedsystems.com/publications +*/ + function q = marcumq (a, b, m, tol) -//This function computes the generalized Marcum Q function of order m with noncentrality parameter a and argument b. -//Calling Sequence -//q = marcumq (a, b) -//q = marcumq (a, b, m) -//q = marcumq (a, b, m, tol) -//Parameters -//a: -//b: -//m: default value 1 -//tol: default value eps -//Description -//This is an Octave function. -//This function computes the generalized Marcum Q function of order m with noncentrality parameter a and argument b. -//The third argument m is the order, which by default is 1. -//The fourth argument tol is the tolerance, which by default is eps. -//If input arguments are vectors which correspond in size and degree, the output is a table of values. -//This function calculates Marcum’s Q function using the infinite Bessel series, which is truncated when the relative error is less than the specified tolerance. -//Examples -//marcumq([1,2,3],4) -//ans = -// 0.0028895 0.0341348 0.1965122 - -funcprot(0); -rhs = argn(2) -if(rhs<2 | rhs>4) -error("Wrong number of input arguments.") -end - select(rhs) - case 2 then - q = callOctave("marcumq",a,b) - case 3 then - q = callOctave("marcumq",a,b,m) - case 4 then - q = callOctave("marcumq",a,b,m,tol) - end + + if ((nargin < 2) || (nargin > 5)) + error(" marcumq : wrong numbers of input arguments "); + end + + if nargin < 3 then m = 1 end + if nargin < 4 then tol = %eps end + if nargin < 5 then max_iter = 100 end + + if (or (a < 0)) + error ("marcumq: A must be a non-negative value"); + end + if (or (b < 0)) + error ("marcumq: B must be a non-negative value"); + end + if (or (m < 1) || or (fix (m) ~= m)) + error ("marcumq: M must be a positive integer"); + end + + [a, b] = tablify (a, b); + q = [] + for i=1:size(a,1) + for j=1:size(a,2) + q(i,j) = mq(a(i,j),b(i,j),m,tol) + end + end + endfunction - + +// Subfunction to compute the actual Marcum Q function. +function q = mq (a, b, m, tol) + + // Special cases. + if (b == 0) + q = 1; + N = 0; + return; + end + if (a == 0) + k = 0:(m - 1); + q = exp (-b^2 / 2) * sum (b.^(2 * k) ./ (2.^k .* factorial (k))); + N = 0; + return; + end + + // The basic iteration. If a<b compute Q_M, otherwise + // compute 1-Q_M. + k = m; + z = a * b; + t = 1; + k = 0; + if (a < b) + s = 1; + c = 0; + x = a / b; + d = x; + S = besseli (0, z, 1); + if (m > 1) + for k = 1:m - 1 + t = (d + 1 / d) * besseli (k, z, 1); + S = S + t; + d = d * x; + end + end + N = k; + k = k + 1 + else + s = -1; + c = 1; + x = b / a; + k = m; + d = x^m; + S = 0; + N = 0; + end + + while (abs (t / S) > tol) + t = d * besseli (abs (k), z, 1); + S = S + t; + d = d * x; + N = k; + k = k + 1 ; + end + q = c + s * exp (-(a - b)^2 / 2) * S; + +endfunction + +// Internal helper function to create a table of like dimensions from arguments. +function [ta , tb] = tablify(a,b) + rows = size(a,1) + cols = size(b,2) + ta=[] + for i=1:cols + ta = [ta a ] + end + tb=[] + for i=1:rows + tb = [ tb ; b] + end +endfunction +/* +test + a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00; + 11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00; + 21.00; 22.00; 23.00; 24.00]; + b = [0.000000, 0.100000, 1.100000, 2.100000, 3.100000, 4.100000]; + Q = [1.000000, 0.995012, 0.546074, 0.110251, 0.008189, 0.000224; + 1.000000, 0.995019, 0.546487, 0.110554, 0.008238, 0.000226; + 1.000000, 0.996971, 0.685377, 0.233113, 0.034727, 0.002092; + 1.000000, 0.999322, 0.898073, 0.561704, 0.185328, 0.027068; + 1.000000, 0.999944, 0.985457, 0.865241, 0.526735, 0.169515; + 1.000000, 0.999998, 0.999136, 0.980933, 0.851679, 0.509876; + 1.000000, 1.000000, 0.999979, 0.998864, 0.978683, 0.844038; + 1.000000, 1.000000, 1.000000, 0.999973, 0.998715, 0.977300; + 1.000000, 1.000000, 1.000000, 1.000000, 0.999969, 0.998618; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; +q = marcumq (a, b); +assert_checkalmostequal (q, Q, %eps,1e-4); + +test + a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00; + 11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00; + 21.00; 22.00; 23.00; 24.00]; + b = [5.100000, 6.100000, 7.100000, 8.100000, 9.100000, 10.10000]; + Q = [0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000049, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.001606, 0.000037, 0.000000, 0.000000, 0.000000, 0.000000; + 0.024285, 0.001420, 0.000032, 0.000000, 0.000000, 0.000000; + 0.161412, 0.022812, 0.001319, 0.000030, 0.000000, 0.000000; + 0.499869, 0.156458, 0.021893, 0.001256, 0.000028, 0.000000; + 0.839108, 0.493229, 0.153110, 0.021264, 0.001212, 0.000027; + 0.976358, 0.835657, 0.488497, 0.150693, 0.020806, 0.001180; + 0.998549, 0.975673, 0.833104, 0.484953, 0.148867, 0.020458; + 0.999965, 0.998498, 0.975152, 0.831138, 0.482198, 0.147437; + 1.000000, 0.999963, 0.998458, 0.974742, 0.829576, 0.479995; + 1.000000, 1.000000, 0.999962, 0.998426, 0.974411, 0.828307; + 1.000000, 1.000000, 1.000000, 0.999961, 0.998400, 0.974138; + 1.000000, 1.000000, 1.000000, 1.000000, 0.999960, 0.998378; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999960; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; + q = marcumq (a, b); + assert_checkalmostequal(q, Q, %eps,1e-4); + + +test + a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00; + 11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00; + 21.00; 22.00; 23.00; 24.00]; + b = [11.10000, 12.10000, 13.10000, 14.10000, 15.10000, 16.10000]; + Q = [0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.000026, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000; + 0.001155, 0.000026, 0.000000, 0.000000, 0.000000, 0.000000; + 0.020183, 0.001136, 0.000025, 0.000000, 0.000000, 0.000000; + 0.146287, 0.019961, 0.001120, 0.000025, 0.000000, 0.000000; + 0.478193, 0.145342, 0.019778, 0.001107, 0.000024, 0.000000; + 0.827253, 0.476692, 0.144551, 0.019625, 0.001096, 0.000024; + 0.973909, 0.826366, 0.475422, 0.143881, 0.019494, 0.001087; + 0.998359, 0.973714, 0.825607, 0.474333, 0.143304, 0.019381; + 0.999959, 0.998343, 0.973546, 0.824952, 0.473389, 0.142803; + 1.000000, 0.999959, 0.998330, 0.973400, 0.824380, 0.472564; + 1.000000, 1.000000, 0.999958, 0.998318, 0.973271, 0.823876; + 1.000000, 1.000000, 1.000000, 0.999958, 0.998307, 0.973158; + 1.000000, 1.000000, 1.000000, 1.000000, 0.999957, 0.998297; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999957; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; + q = marcumq (a, b); + assert_checkalmostequal (q, Q,%eps,1e-4); + + +test + a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00; + 11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00; + 21.00; 22.00; 23.00; 24.00]; + b = [17.10000, 18.10000, 19.10000]; + Q = [0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000000, 0.000000, 0.000000; + 0.000024, 0.000000, 0.000000; + 0.001078, 0.000024, 0.000000; + 0.019283, 0.001071, 0.000023; + 0.142364, 0.019197, 0.001065; + 0.471835, 0.141976, 0.019121; + 0.823429, 0.471188, 0.141630; + 0.973056, 0.823030, 0.470608; + 0.998289, 0.972965, 0.822671; + 0.999957, 0.998281, 0.972883; + 1.000000, 0.999957, 0.998274; + 1.000000, 1.000000, 0.999956; + 1.000000, 1.000000, 1.000000]; + q = marcumq (a, b); + assert_checkalmostequal(q, Q, %eps,1e-4); + +// The tests for M>1 were generating from Marcum's tables by +// using the formula +// Q_M(a,b) = Q(a,b) + exp(-(a-b)^2/2)*sum_{k=1}^{M-1}(b/a)^k*exp(-ab)*I_k(ab) + +test + M = 2; + a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00; + 11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00; + 21.00; 22.00; 23.00; 24.00]; + b = [ 0.00, 0.10, 2.10, 7.10, 12.10, 17.10]; + Q = [1.000000, 0.999987, 0.353353, 0.000000, 0.000000, 0.000000; + 1.000000, 0.999988, 0.353687, 0.000000, 0.000000, 0.000000; + 1.000000, 0.999992, 0.478229, 0.000000, 0.000000, 0.000000; + 1.000000, 0.999999, 0.745094, 0.000001, 0.000000, 0.000000; + 1.000000, 1.000000, 0.934771, 0.000077, 0.000000, 0.000000; + 1.000000, 1.000000, 0.992266, 0.002393, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999607, 0.032264, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999992, 0.192257, 0.000000, 0.000000; + 1.000000, 1.000000, 1.000000, 0.545174, 0.000000, 0.000000; + 1.000000, 1.000000, 1.000000, 0.864230, 0.000040, 0.000000; + 1.000000, 1.000000, 1.000000, 0.981589, 0.001555, 0.000000; + 1.000000, 1.000000, 1.000000, 0.998957, 0.024784, 0.000000; + 1.000000, 1.000000, 1.000000, 0.999976, 0.166055, 0.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 0.509823, 0.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 0.846066, 0.000032; + 1.000000, 1.000000, 1.000000, 1.000000, 0.978062, 0.001335; + 1.000000, 1.000000, 1.000000, 1.000000, 0.998699, 0.022409; + 1.000000, 1.000000, 1.000000, 1.000000, 0.999970, 0.156421; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.495223; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.837820; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.976328; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.998564; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; + q = marcumq (a, b, M); + assert_checkalmostequal (q, Q, %eps,1e-4); + +test + + M = 5; + a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00; + 11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00; + 21.00; 22.00; 23.00; 24.00]; + b = [ 0.00, 0.10, 2.10, 7.10, 12.10, 17.10]; + Q = [1.000000, 1.000000, 0.926962, 0.000000, 0.000000, 0.000000; + 1.000000, 1.000000, 0.927021, 0.000000, 0.000000, 0.000000; + 1.000000, 1.000000, 0.947475, 0.000001, 0.000000, 0.000000; + 1.000000, 1.000000, 0.980857, 0.000033, 0.000000, 0.000000; + 1.000000, 1.000000, 0.996633, 0.000800, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999729, 0.011720, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999990, 0.088999, 0.000000, 0.000000; + 1.000000, 1.000000, 1.000000, 0.341096, 0.000000, 0.000000; + 1.000000, 1.000000, 1.000000, 0.705475, 0.000002, 0.000000; + 1.000000, 1.000000, 1.000000, 0.933009, 0.000134, 0.000000; + 1.000000, 1.000000, 1.000000, 0.993118, 0.003793, 0.000000; + 1.000000, 1.000000, 1.000000, 0.999702, 0.045408, 0.000000; + 1.000000, 1.000000, 1.000000, 0.999995, 0.238953, 0.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 0.607903, 0.000001; + 1.000000, 1.000000, 1.000000, 1.000000, 0.896007, 0.000073; + 1.000000, 1.000000, 1.000000, 1.000000, 0.987642, 0.002480; + 1.000000, 1.000000, 1.000000, 1.000000, 0.999389, 0.034450; + 1.000000, 1.000000, 1.000000, 1.000000, 0.999988, 0.203879; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.565165; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.876284; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.984209; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999165; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999983; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; + q = marcumq (a, b, M); + assert_checkalmostequal (q, Q, %eps,1e-4); + +test passed + M = 10; + a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00; + 11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00; + 21.00; 22.00; 23.00; 24.00]; + b = [ 0.00, 0.10, 2.10, 7.10, 12.10, 17.10]; + Q = [1.000000, 1.000000, 0.999898, 0.000193, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999897, 0.000194, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999931, 0.000416, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999980, 0.002377, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999997, 0.016409, 0.000000, 0.000000; + 1.000000, 1.000000, 0.999999, 0.088005, 0.000000, 0.000000; + 1.000000, 1.000000, 1.000000, 0.302521, 0.000000, 0.000000; + 1.000000, 1.000000, 1.000000, 0.638401, 0.000000, 0.000000; + 1.000000, 1.000000, 1.000000, 0.894322, 0.000022, 0.000000; + 1.000000, 1.000000, 1.000000, 0.984732, 0.000840, 0.000000; + 1.000000, 1.000000, 1.000000, 0.998997, 0.014160, 0.000000; + 1.000000, 1.000000, 1.000000, 0.999972, 0.107999, 0.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 0.391181, 0.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 0.754631, 0.000004; + 1.000000, 1.000000, 1.000000, 1.000000, 0.951354, 0.000266; + 1.000000, 1.000000, 1.000000, 1.000000, 0.995732, 0.006444; + 1.000000, 1.000000, 1.000000, 1.000000, 0.999843, 0.065902; + 1.000000, 1.000000, 1.000000, 1.000000, 0.999998, 0.299616; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.676336; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.925312; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.992390; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999679; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999995; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000; + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; + q = marcumq (a, b, M); + assert_checkalmostequal (q, Q, %eps,1e-4); + +*/
\ No newline at end of file diff --git a/macros/mpoles.sci b/macros/mpoles.sci index bba9997..4556120 100644 --- a/macros/mpoles.sci +++ b/macros/mpoles.sci @@ -1,124 +1,128 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Author:[insert name] +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in +function [multp, idxp] = mpoles (p, tol, reorder) + if (nargin < 1) + error("mpoles: Invalid number of input arguments"); + end -function [multp, indx] = mpoles (p, tol, reorder) - -//Calling sequence: -// [multp, idxp] = mpoles (p) -// [multp, idxp] = mpoles (p, tol) -// [multp, idxp] = mpoles (p, tol, reorder) -// Identify unique poles in p and their associated multiplicity. -// -// The output is ordered from largest pole to smallest pole. -// -// If the relative difference of two poles is less than tol then they are -// considered to be multiples. The default value for tol is 0.001. -// -// If the optional parameter reorder is zero, poles are not sorted. -// -// The output multp is a vector specifying the multiplicity of the poles. -// multp(n) refers to the multiplicity of the Nth pole -// p(idxp(n)). -// -// -// - -// test case: - -// p = [2 3 1 1 2]; -// [m, n] = mpoles (p) -// n = -// 2. -// 5. -// 1. -// 4. -// 3. -// m = -// 1. -// 1. -// 2. -// 1. -// 2. -// - - -[nargout,nargin]=argn(); - - if (nargin < 1 | nargin > 3) - error("wrong number of input arguments"); + if ~( type(p)== 1) + error ("mpoles: P must be a single or double floating point vector"); end - if (nargin < 2 | isempty (tol)) - tol = 0.001; - end + if (nargin < 2 || isempty (tol)) + tol = 0.001; + elseif (~(isscalar (tol) && isreal (tol) && tol > 0)) + error ("mpoles: TOL must be a real scalar greater than 0"); + end - if (nargin < 3 | isempty (reorder)) - reorder = %t; - end + if (nargin < 3 || isempty (reorder)) + reorder = %t; + elseif (~(isscalar (reorder) && isreal (reorder))) + error ("mpoles: REORDER must be a numeric or logical scalar"); + end Np = length (p); - - // Force the poles to be a column vector. - - p = p(:); - - // Sort the poles according to their magnitidues, largest first. + p = p(:); // force poles to be a column vector if (reorder) - // Sort with smallest magnitude first. - [p, ordr] = gsort (p,"r","i"); - // Reverse order, largest maginitude first. - n = Np:-1:1; - p = p(n); - ordr = ordr(n); + //// sort with largest magnitude first + [_, order] = gsort (abs(p)); + p = p(order); else - ordr = 1:Np; + order = (1:Np).'; end - // Find pole multiplicty by comparing the relative differnce in the - // poles. + //// Create vector of tolerances for use in algorithm. + vtol = zeros (Np, 1, typeof(p)); + p_nz = (p ~= 0); // non-zero poles + vtol(~p_nz) = tol; // use absolute tolerance for zero poles - multp = zeros (Np, 1); - indx = []; + //// Find pole multiplicity by comparing relative difference of poles. + multp = zeros (Np, 1, typeof(p)); + idxp = []; n = find (multp == 0, 1); - - - - - - while (n) - dp = abs (p-p(n)); - if (p(n) == 0.0) - if (or (abs (p) > 0 & (abs(p)<%inf))) - p0 = mean (abs (p(abs (p) > 0 & abs(p)<%inf))); - else - p0 = 1; - end - else - p0 = abs (p(n)); + dp = abs (p - p(n)); + vtol(p_nz) = tol * abs (p(n)); + k = find (dp < vtol); + //// Poles can only be members of one multiplicity group. + if (length (idxp)) + k = k(~ismember (k, idxp)); end - k = find (dp < tol * p0)'; - // Poles can only be members of one multiplicity group. -// if (length(indx)) -// mk=members(k,indx); -// k = (~ bool2s(mk~=0)); -// end m = 1:length (k); - multp(k) = m'; - indx = [indx; k]; + multp(k) = m; + // disp("k") + // disp(k) + // disp("idxp") + // disp(idxp) + idxp = [idxp; k(:)]; n = find (multp == 0, 1); end - multp = multp(indx); - indx = ordr(indx); + multp = multp(idxp); + idxp = order(idxp); +endfunction +function y = ismember(a, b) // FIXME : do this in a more efficient . + y = zeros(size(a,1),size(a,2)); + for i = 1:length(a) + for j = 1:length(b) + if a(i) == b(j) + y(i) = 1; + break; + end + end + end endfunction + + +/* + + // test + [mp, ip] = mpoles ([0 0], 0.01); + assert_checkequal (mp, [1; 2]); + + // test + [mp, ip] = mpoles ([-1e4, -0.1, 0]); + assert_checkequal (mp, [1; 1; 1]); + assert_checkequal (ip, [1; 2; 3]); + +// Test single inputs +// test + [mp, ip] = mpoles ([-1e4, -0.1, 0]); + assert_checkequal (mp,[1; 1; 1]); + assert_checkequal (ip, [1; 2; 3]); + +// Test relative tolerance criteria +// test + [mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1); + assert_checkequal (mp, [1; 1; 1]); + [mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1 + %eps); + assert_checkequal (mp, [1; 1; 2]); + +// Test absolute tolerance criteria with a zero pole +// test + [mp, ip] = mpoles ([0, -0.1, 0.3], .1); + assert_checkequal (mp, [1; 1; 1]); + [mp, ip] = mpoles ([0, -0.1, 0.3], .1 + %eps); + assert_checkequal (mp, [1; 1; 2]); + +//// Test input validation +error <Invalid call> mpoles () +error <P must be a single or double floating point vector> mpoles (uint8 (1)) +error <TOL must be a real scalar greater than 0> mpoles (1, [1, 2]) +error <TOL must be a real scalar greater than 0> mpoles (1, 1i) +error <TOL must be a real scalar greater than 0> mpoles (1, 0) +error <REORDER must be a numeric or logical scalar> mpoles (1, 1, [1, 2]) +error <REORDER must be a numeric or logical scalar> mpoles (1, 1, {1}) + +*/
\ No newline at end of file diff --git a/macros/ncauer.sci b/macros/ncauer.sci index d9a0684..dad3025 100644 --- a/macros/ncauer.sci +++ b/macros/ncauer.sci @@ -1,53 +1,85 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Author:Sonu Sharma, RGIT Mumbai +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - function [Zz, Zp, Zg] = ncauer(Rp, Rs, n) -//Analog prototype for Cauer filter (Cauer filter and elliptic filters are same). - -//Calling Sequence -//[Zz, Zp, Zg] = ncauer(Rp, Rs, n) - -//Parameters -//n: Filter Order -//Rp: Peak-to-peak passband ripple in dB -//Rs: Stopband attenuation in dB - -//Description -//It gives an analog prototype for Cauer filter of nth order, with a Peak-to-peak passband ripple of Rp dB and a stopband attenuation of Rs dB. + //Analog prototype for Cauer filter (Cauer filter and elliptic filters are same). + + //Calling Sequence + //[Zz, Zp, Zg] = ncauer(Rp, Rs, n) + + //Parameters + //n: Filter Order + //Rp: Peak-to-peak passband ripple in dB + //Rs: Stopband attenuation in dB + + //Description + //It gives an analog prototype for Cauer filter of nth order, with a Peak-to-peak passband ripple of Rp dB and a stopband attenuation of Rs dB. + + + //Examples + //n = 5; + //Rp = 5; + //Rs = 5; + //[Zz, Zp, Zg] = ncauer(Rp, Rs, n) + + //Zz = + // + // 0.0000 + 2.5546i 0.0000 + 1.6835i -0.0000 - 2.5546i -0.0000 - 1.6835i + // + //Zp = + // + // -0.10199 + 0.64039i -0.03168 + 0.96777i -0.10199 - 0.64039i -0.03168 - 0.96777i -0.14368 + 0.00000i + // + //Zg = 0.0030628 + // Dependencies + // ellipap + + funcprot(0); + lhs = argn(1) + rhs = argn(2) + if (rhs < 3 | rhs > 3) + error("ncauer : Wrong number of input arguments.") + end + + [Zz, Zp, Zg] = ellipap(n, Rp, Rs) ; + // temp fix to permanently fix this change ellipap + Zz = Zz'; + Zp = Zp'; + endfunction -//Examples -//n = 5; -//Rp = 5; -//Rs = 5; -//[Zz, Zp, Zg] = ncauer(Rp, Rs, n) +/* -//Zz = -// -// 0.0000 + 2.5546i 0.0000 + 1.6835i -0.0000 - 2.5546i -0.0000 - 1.6835i -// -//Zp = -// -// -0.10199 + 0.64039i -0.03168 + 0.96777i -0.10199 - 0.64039i -0.03168 - 0.96777i -0.14368 + 0.00000i -// -//Zg = 0.0030628 +// Test Case 1 (ncauer 0.1, 60, 7) +[z, p, g] = ncauer(0.1, 60, 7); +assert_checkalmostequal(z, [2.5574 1.5522 1.3295 -2.5574 -1.5522 -1.3295]*%i, 1e-2); +assert_checkalmostequal(p, [-0.3664+0.5837*%i -0.1802+0.9088*%i -0.0499+1.0285*%i -0.3664-0.5837*%i -0.1802-0.9088*%i -0.0499-1.0285*%i -0.4796], 1e-2); +assert_checkalmostequal(g, 7.4425e-03, 1e-2); +// Test Case 2 (ncauer 1.0, 30, 3) +[z, p, g] = ncauer(1.0, 30, 3); +assert_checkalmostequal(z, [1.9536 -1.9536]*%i, 1e-2); +assert_checkalmostequal(p, [-0.2053+0.9870*%i -0.2053-0.9870*%i -0.5597], 1e-2); +assert_checkalmostequal(g, 0.1490, 1e-2); -funcprot(0); -lhs = argn(1) -rhs = argn(2) -if (rhs < 3 | rhs > 3) -error("ncauer : Wrong number of input arguments.") -end +// Test Case 3 (ncauer 0.25, 50, 6) +[z, p, g] = ncauer(0.25, 50, 6); +assert_checkalmostequal(z, [4.0596 1.6414 1.3142 -4.0596 -1.6414 -1.3142]*%i, 1e-2); +assert_checkalmostequal(p, [-0.4210+0.3665*%i -0.2117+0.8503*%i -0.0550+1.0198*%i -0.4210-0.3665*%i -0.2117-0.8503*%i -0.0550-1.0198*%i], 1e-2); +assert_checkalmostequal(g, 3.1618e-03, 1e-2); -[Zz, Zp, Zg] = ellipap(n, Rp, Rs) ; +// Test Case 4 (ncauer 0.8, 45, 4) +[z, p, g] = ncauer(0.8, 45, 4); +assert_checkalmostequal(z, [4.1768 1.8543 -4.1768 -1.8543]*%i, 1e-2); +assert_checkalmostequal(p, [-0.3861+0.4640*%i -0.1234+1.0000*%i -0.3861-0.4640*%i -0.1234-1.0000*%i], 1e-2); +assert_checkalmostequal(g, 5.6237e-03, 1e-2); -endfunction +*/
\ No newline at end of file diff --git a/macros/parser.sci b/macros/parser.sci new file mode 100644 index 0000000..9611161 --- /dev/null +++ b/macros/parser.sci @@ -0,0 +1,82 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +/// Author : Abinash Singh Under FOSSEE Internship +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +function [ dSided, minH, minD, minW, maxW ] = parser ( varargin ) + // Default values + // This is an helper function for findpeaks + // It parses the input arguments and returns the values of the options + dSided = %f ; + minH = %eps ; + minD = 1 ; + minW = %eps; + maxW = %inf ; + idx = 1 ; + while idx <= length(varargin) + select lower(varargin(idx)) + case 'doublesided' + dSided = %t ; + idx = idx + 1 ; + case 'minpeakheight' + if idx+1 > length(varargin) || ~( isscalar(varargin(idx+1)) && varargin(idx+1) >=0 ) then + error('findpeaks: MinPeakHeight must be a postive scalar') ; + end + minH = varargin(idx+1) ; + idx = idx + 2 ; + case 'minpeakdistance' + if idx+1 > length(varargin) || ~( isscalar(varargin(idx+1)) && varargin(idx+1) >=0 ) then + error('findpeaks: MinPeakDistance must be a postive scalar') ; + end + minD = varargin(idx+1) ; + idx = idx + 2 ; + case 'minpeakwidth' + if idx+1 > length(varargin) || ~( isscalar(varargin(idx+1)) && varargin(idx+1) >=0 ) then + error('findpeaks: MinPeakWidth must be a postive scalar') ; + end + minW = varargin(idx+1) ; + idx = idx + 2 ; + case 'maxpeakwidth' + if idx+1 > length(varargin) || ~( isscalar(varargin(idx+1)) && varargin(idx+1) >=0 ) then + error('findpeaks: MaxPeakWidth must be a postive scalar') ; + end + maxW = varargin(idx+1) ; + idx = idx + 2 ; + else + warning("findpeaks: Ignoring unknown option ") ; + idx = idx + 1 ; + end + + end +endfunction + +function y = lower (y) + // This function converts the input string to lower case + // Returns the input without modification if its not a string + if type(y) == 10 then + y = convstr(y, 'l') ; + else + y = y ; + end + +endfunction + +function out = bsminuseq(A) + // This function returns the difference between the elements of the input array + // This is only useful for the findpeaks function + A = A(:).' ; + atemp = []; + btemp=[] ; + for i=1:length(A) + atemp = [atemp ; A ] + btemp = [btemp A.' ] + end + out = atemp - btemp +endfunction diff --git a/macros/periodgram.sci b/macros/periodgram.sci new file mode 100644 index 0000000..0d406ec --- /dev/null +++ b/macros/periodgram.sci @@ -0,0 +1,232 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +function [pxx, f] = periodogram (x, varargin) + //Calling Sequence: + //[PXX, F] = periodogram (X, WIN, NFFT, FS) + //[PXX, F] = periodogram (..., "RANGE") + + // The possible inputs are: + // + // X + // + // data vector. If X is real-valued a one-sided spectrum is + // estimated. If X is complex-valued, or "RANGE" specifies + // "twosided", the full spectrum is estimated. + // + // WIN + // window weight data. If window is empty or unspecified a + // default rectangular window is used. Otherwise, the window is + // applied to the signal ('X .* WIN') before computing the + // periodogram. The window data must be a vector of the same + // length as X. + // + // NFFT + // number of frequency bins. The default is 256 or the next + // higher power of 2 greater than the length of X ('max (256, + // 2.^nextpow2 (length (x)))'). If NFFT is greater than the + // length of the input then X will be zero-padded to the length + // of NFFT. + // + // FS + // sampling rate. The default is 1. + // + // RANGE + // range of spectrum. "onesided" computes spectrum from + // [0..nfft/2+1]. "twosided" computes spectrum from [0..nfft-1]. + // + // + // Dependencies + // hamming fft1 + [nargout,nargin]=argn(); + // check input arguments + if (nargin < 1 | nargin > 5) + error("wrong no. of input arguments") + end + + nfft = []; + fs = []; + ran = ""; + win = []; + j = 2; + for k = 1:length (varargin) + if (type (varargin(k))==10) + ran = varargin(k); + else + select (j) + case 2 + win = varargin(k); + case 3 + nfft = varargin(k); + case 4 + fs = varargin(k); + end + j=j+1; + end + end + + if (~ isvector (x)) + error ("periodogram: X must be a real or complex vector"); + end + x = x(:); // Use column vectors from now on + + n = size(x,1); + + if (~isempty (win)) + if (~ isvector (win) | length (win) ~= n) + error ("periodogram: WIN must be a vector of the same length as X"); + end + win = win(:); + x =x.* win; + else + win=window("re",length(x)); + win=win(:); + x=x.*win; + + end + + if (isempty (nfft)) + nfft = max (256, 2.^nextpow2 (n)); + elseif (~ isscalar (nfft)) + error ("periodogram: NFFT must be a scalar"); + end + + use_w_freq = isempty (fs); + if (~use_w_freq & ~ isscalar (fs)) + error ("periodogram: FS must be a scalar"); + end + + if (~strcmp (ran, "onesided")) + ran = 1; + elseif (~strcmp (ran, "twosided")) + ran = 2; + elseif (~strcmp (ran, "centered")) + error ('periodogram: centered ran type is not implemented'); + else + ran = 2 - double(isreal (x)); + end + + // compute periodogram + + if (n > nfft) + Pxx = 0; + rr = modulo(length(x), nfft); + if (rr) + x = [x(:); zeros(nfft-rr, 1)]; + end + x = sum (matrix (x, nfft,-1), 2); + end + + if (~ isempty (win)) + n = sum(win.*conj(win)); + end; + + Pxx = (abs (fft1 (x,nfft))) .^ 2 / n; + + if (use_w_freq) + Pxx =Pxx/(2*%pi); + else + Pxx =Pxx/fs; + end + + // generate output arguments + + if (ran == 1) // onesided + if (modulo(nfft,2)==0) // nfft is even + psd_len = (nfft/2)+1; + Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1); 0]; + else // nfft is odd + psd_len = (nfft+1)/2; + Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1)]; + end + end + + //if (nargout() ~= 1) FIXME: fix nargout + if (ran == 1) + f = (0:nfft/2)' / nfft; + elseif (ran == 2) + f = (0:nfft-1)' / nfft; + end + if (use_w_freq) + f =f* 2*%pi; // generate w=2*pi*f + else + f =f* fs; + end + //end + if (nargout() ~= 2) + if (use_w_freq) + plot (f/(2*%pi), 10*log10 (Pxx)); + xlabel ("normalized frequency [x pi rad]"); + ylabel ("Power density [dB/rad/sample]"); + else + plot (f, 10*log10 (Pxx)); + xlabel ("frequency [Hz]"); + ylabel ("Power density [dB/Hz]"); + end + title ("Periodogram Power Spectral Density Estimate"); + end + pxx = Pxx; +endfunction + +/* +pi = %pi; // ezecute on scilab only + +t = 0:0.01:1; +x = sin(2*pi*10*t); +periodogram(x); + +x = complex(0:0.01:1, 0:0.01:1); +periodogram(x); + +x = cos(0:0.01:1); +win = hamming(101); +periodogram(x, win); + + + +x = tan(0:0.01:1); +nfft = 512; +periodogram(x, [], nfft); + + +t = 0:0.01:1; +x = sin(2*pi*10*t); +Fs = 100; +periodogram(x, [], [], Fs) + + + +x = sin(0:0.01:1); +periodogram(x, [], [], [], 'onesided'); + +x = sin(0:0.01:1); +periodogram(x, [], [], [], 'twosided') + + +Fs = 1000; +t = 0:1/Fs:1-1/Fs; +f0 = 100; +x = sin(2*pi*f0*t); +[Pxx, f] = periodogram(x, [], [], Fs); +[_, idx] = max(Pxx); +detected_freq = f(idx); + + +// Test error : invalid window length +x = randn(100,1); +win = hamming(50); +periodogram(x, win) + +// Test invalid nfft (negative) +periodogram(x, [], -256) + +*/ diff --git a/macros/postpad.sci b/macros/postpad.sci index f8042be..3797e08 100644 --- a/macros/postpad.sci +++ b/macros/postpad.sci @@ -1,49 +1,57 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Author:Sonu Sharma, RGIT Mumbai +// Author: Abinash Singh Under FOSSEE Internship +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - -// This is a supporting function - -function y = postpad(x, n, varargin) - - //Calling Sequences - // Y = postpad (X, L) - // Y = postpad (X, L, C) - - // Description : - // Append the scalar value C to the vector X until it is of length L. - // If C is not given, a value of 0 is used. - // - // If 'length (X) > L', elements from the end of X are removed until a - // vector of length L is obtained. - - //Example : - //x = [1 2 3]; - //L = 6; - //y = postpad(x, L) - //Output : - // y = - // - // 1. 2. 3. 0. 0. 0. - - funcprot(0); - if argn(2) > 3 | argn(2) < 2 then - error("postpad : wrong number of input argument ") - elseif argn(2) == 2 +/* +Calling Sequence : + postpad (x, l) + postpad (x, l, c) + postpad (x, l, c, dim) +Append the scalar value c to the vector x until it is of length l. If c is not given, a value of 0 is used. +If length (x) > l, elements from the end of x are removed until a vector of length l is obtained. +If x is a matrix, elements are appended or removed from each row. +If the optional argument dim is given, operate along this dimension. +If dim is larger than the dimensions of x, the result will have dim dimensions. +*/ +function res = postpad(x,l,c,dim) + if nargin < 2 then + error("Usage : postpad(x,l,c(optional),dim(optional))") + end + if nargin < 3 then c = 0 ; - else - c = varargin(1); end - - y = x; - for i = 1:(n-length(x)) - y = [y c]; + if nargin <= 4 then + if size(x,1) == 1 then + dim = 1 ; + elseif size(x,2) == 1 then + dim = 2 ; + else + dim = 2 ; // FIXME dim functionality not implemented + end + end + if l < size(x,dim) then + error("l must be greater then dimension of x") + end + + select dim + case 1 then + res = [x c*ones(size(x,1),l-size(x,2))]; + case 2 then + res = [x;c*ones(l-size(x,1),size(x,2))]; end endfunction + +/* +#test for row vectors +postpad([1 2 3 4],6) //passed +postpad([1 ;2 ;3 ;4],6) // passed +postpad([1 2 3 4;5 6 7 8;9 10 11 12],6) // passed +postpad([1 2 ;3 4;5 6],6,-1) //passed +*/
\ No newline at end of file diff --git a/macros/prepad.sci b/macros/prepad.sci index 4e517f8..d3595bf 100644 --- a/macros/prepad.sci +++ b/macros/prepad.sci @@ -1,83 +1,65 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Author:[insert name] +// Author: Abinash Singh Under FOSSEE Internship +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in +/* +Calling Sequence : + prepad (x, l) + prepad (x, l, c) + prepad (x, l, c, dim) - -function y = prepad (x, l, c, dim) -// Calling sequence: -// y= prepad (x, l) -// y= prepad (x, l, c) -// y= prepad (x, l, c, dim) -// Prepend the scalar value c to the vector x until it is of length -// l. If c is not given, a value of 0 is used. -// -// If length (x) > l, elements from the beginning of x -// are removed until a vector of length l is obtained. -// -// If x is a matrix, elements are prepended or removed from each row. -// -// If the optional argument dim is given, operate along this dimension. -// -// If dim is larger than the dimensions of x, the result will have -// dim dimensions. - -//Test cases: -//prepad ([1,2], 4,0,2) -//Output: [0,0,1,2] - - -[nargout,nargin]=argn(); - if (nargin < 2 | nargin > 4) - error("wrong number of input arguments"); - end - - if (nargin < 3 | isempty (c)) - c = 0; - else - if (~ isscalar (c)) - error ("prepad: pad value C must be empty or a scalar"); +Prepend the scalar value c to the vector x until it is of length l. If c is not given, a value of 0 is used. +If length (x) > l, elements from the beginning of x are removed until a vector of length l is obtained. +If x is a matrix, elements are prepended or removed from each row. +If the optional argument dim is given, operate along this dimension. +If dim is larger than the dimensions of x, the result will have dim dimensions. +*/ +function res = prepad(x,l,c,dim) + if nargin < 2 then + error("Usage : postpad(x,l,c(optional),dim(optional))") end - end -//dim=1; - nd = ndims (x); - sz = size (x); - if (argn(2) < 4) - // Find the first non-singleton dimension. - dim = find (sz > 1, 1) - if (dim == []) - dim = 1 + if nargin < 3 then + c = 0 ; end - else - if (~(isscalar (dim) & dim == fix (dim) & dim >= 1)) - error ("prepad: DIM must be an integer and a valid dimension"); + if nargin < 4 then + dim = find(size(x)>1) + if isempty(dim) then dim = 2 end + if isvector(dim) then dim = dim(1) end end - end - - if (~ isscalar (l) | l < 0) - error ("prepad: length L must be a positive scalar"); - end - if (dim > nd) - sz(nd+1:dim) = 1; - end - - d = sz(dim); + if l < size(x,dim) then + if isvector(x) then + start = length(x) - l+1; + res=x(start:$) + return; + else + error("prepad : l must be greter than dim size for matrices") + end + + end + + if dim == 1 then + res = [ c* ones( l - size(x,1) , size(x,2)) ; x] + elseif dim == 2 then + res = [ c* ones( size(x,1) , l - size(x,2)) x] + else + error("prepad : Invalid value for arg dim 1 or 2 expected") + end +endfunction - if (d >= l) -// idx = repmat ({':'}, nd, 1); -// idx(dim) = (d-l+1) : d; -// y = x(idx(:)); - y=x; - else - sz(dim) = l - d; - y = cat (dim, c(ones (sz(1),sz(2))), x); - end +/* +#test for row vectors +prepad([1 2 3 4],6) +prepad([1 ;2 ;3 ;4],6) +prepad([1 2 3 4;5 6 7 8;9 10 11 12],6) +prepad([1 2 ;3 4;5 6],6,-1) -endfunction +// FIXME : Tests for 2d and high dimesnsional matrices +*/ diff --git a/macros/residue.sci b/macros/residue.sci index 70dc91f..58303d8 100644 --- a/macros/residue.sci +++ b/macros/residue.sci @@ -1,286 +1,318 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Author:[insert name] +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - - +// Dependencies +// prepad deconv mpoles function [r, p, k, e] = residue (b, a, varargin) - -// [r, p, k, e] = residue (b, a) -// [b, a] = residue (r, p, k) -// [b, a] = residue (r, p, k, e) -// The first calling form computes the partial fraction expansion for the -// quotient of the polynomials, b and a. -// -// The quotient is defined as - -// B(s) M r(m) N -// ---- = SUM ------------- + SUM k(i)*s^(N-i) -// A(s) m=1 (s-p(m))^e(m) i=1 - -// where M is the number of poles (the length of the r, p, -// and e), the k vector is a polynomial of order N-1 -// representing the direct contribution, and the e vector specifies the -// multiplicity of the m-th residue's pole. -// -//NOTE that the polynomials 'b' and 'a' should have real coefficients(because of the function 'filter' used in polyval) -// -//Test case -//1. -// b = [1, 1, 1]; -// a = [1, -5, 8, -4]; -// [r, p, k, e] = residue (b, a) -// result r = [-2; 7; 3] -// result p = [2; 2; 1] -// result k = [](0x0) -// result e = [1; 2; 1] -// - -//2. -//[r,p,k,e]=residue([1 2 1],[1 -5 8 -4]) -//OUTPUT -//r = -// -3.0000 -// 9.0000 -// 4.0000 -// -//p = -// 2.0000 -// 2.0000 -// 1.0000 -// -//f = [](0x0) -//e = -// 1 -// 2 -// 1 -// - - - [nargout,nargin]=argn(); - - if (nargin < 2 | nargin > 4) - error ("wrong umber of input arguments"); - end - - toler = .001; - - if (nargin >= 3) - if (nargin >= 4) - e = varargin(2); - else - e = []; + if (nargin < 2 || nargin > 4) + error("residue: Invalid number of arguments"); end - // The inputs are the residue, pole, and direct part. Solve for the - // corresponding numerator and denominator polynomials - [r, p] = rresidue (b, a, varargin(1), toler, e); - return; - end - - // Make sure both polynomials are in reduced form. - - a = polyreduce (a); - b = polyreduce (b); - - b = b / a(1); - a = a / a(1); - - la = length (a); - lb = length (b); - - // Handle special cases here. - - if (la == 0 | lb == 0) - k =[]; - r = []; - p = []; - e = []; - return; - elseif (la == 1) - k = b / a; - r = []; - p = []; - e = []; - return; - end - - // Find the poles. - - p = roots (a); - lp = length (p); - - // Sort poles so that multiplicity loop will work. - - [e, indx] = mpoles (p, toler, 0); - p = p(indx); - - // For each group of pole multiplicity, set the value of each - // pole to the average of the group. This reduces the error in - // the resulting poles. - - p_group = cumsum (e == 1); - for ng = 1:p_group($) - m = find (p_group == ng); - p(m) = mean (p(m)); - end - - // Find the direct term if there is one. - - if (lb >= la) - // Also return the reduced numerator. - [k, b] = deconv (b, a); + + tol = .001; + + if (nargin >= 3) + if (nargin >= 4) + e = varargin(2); + else + e = []; + end + // The inputs are the residue2, pole, and direct part. + // Solve for the corresponding numerator and denominator polynomials. + [r, p] = rresidue (b, a, varargin(1), tol, e); + return; + end + + // Make sure both polynomials are in reduced form, and scaled. + a = polyreduce (a); + b = polyreduce (b); + + b = b / a(1); + a = a/ a(1); + + la = length (a); lb = length (b); - else - k = []; - end - - // Determine if the poles are (effectively) zero. - - small = max (abs (p)); - if (type(a)==1 | type(b)==1) - small = max ([small, 1]) * 1.1921e-07 * 1e4 * (1 + length (p))^2; - else - small = max ([small, 1]) * %eps * 1e4 * (1 + length (p))^2; - end - p(abs (p) < small) = 0; - - // Determine if the poles are (effectively) real, or imaginary. - - index = (abs (imag (p)) < small); - p(index) = real (p(index)); - index = (abs (real (p)) < small); - p(index) = 1*%i * imag (p(index)); - - // The remainder determines the residues. The case of one pole - // is trivial. - - if (lp == 1) - r = polyval (b, p); - return; - end - - // Determine the order of the denominator and remaining numerator. - // With the direct term removed the potential order of the numerator - // is one less than the order of the denominator. - - aorder = length (a) - 1; - border = aorder - 1; - - // Construct a system of equations relating the individual - // contributions from each residue to the complete numerator. - - A = zeros (border+1, border+1); - B = prepad (matrix (b, [length(b), 1]), border+1, 0,2); - for ip = 1:length (p) - ri = zeros (size (p,1),size(p,2)); - ri(ip) = 1; - A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0,2).'; - end - - // Solve for the residues. -if(size(A,1)~=size(B,1)) - if(size(A,1)<size(B,1)) - A=[A;zeros((size(B,1)-size(A,1)),(size(A,2)))]; + + // Handle special cases here. + if (la == 0 || lb == 0) + k = []; + r = []; + p = []; + e = []; + return; + elseif (la == 1) + k = b / a; + r = [] + p = [] + e = []; + return; + end + + // Find the poles. + p = roots (a); + lp = length (p); + + // Sort poles so that multiplicity loop will work. + [e, idx] = mpoles (p, tol, 1); + p = p(idx); + + // For each group of pole multiplicity, set the value of each + // pole to the average of the group. This reduces the error in + // the resulting poles. + p_group = cumsum (e == 1); + for ng = 1:p_group($) + m = find (p_group == ng); + p(m) = mean (p(m)); + end + + // Find the direct term if there is one. + if (lb >= la) + // Also return the reduced numerator. + [k, b] = deconv (b, a); + lb = length (b); else - B=[zeros((size(A,1)-size(B,1)),(size(B,2)));B]; + k = []; end - + + // Determine if the poles are (effectively) zero. + small = max (abs (p)); + if ( (type(a)==1) || (type(b)==1)) + small = max ([small, 1]) * %eps * 1e4 * (1 + length (p))^2; + else + small = max ([small, 1]) * %eps * 1e4 * (1 + length (p))^2; end - r = A \ B; - r=r(:,$); - -endfunction - -function [pnum, pden, e] = rresidue (rm, p, k, toler, e) - - // Reconstitute the numerator and denominator polynomials from the - // residues, poles, and direct term. -[nargout,nargin]=argn(); - if (nargin < 2 | nargin > 5) - error ("wrong number of input arguments"); - end - - if (nargin < 5) - e = []; - end - - if (nargin < 4) - toler = []; - end - - if (nargin < 3) - k = []; - end - - if (length (e)) - indx = 1:length (p); - else - [e, indx] = mpoles (p, toler, 0); - p = p(indx); - rm = rm(indx); - end - - indx = 1:length (p); - - for n = indx - pn = [1, -p(n)]; - if (n == 1) - pden = pn; + p(abs (p) < small) = 0; + + // Determine if the poles are (effectively) real, or imaginary. + idx = (abs (imag (p)) < small); + p(idx) = real (p(idx)); + idx = (abs (real (p)) < small); + p(idx) = 1*%i * imag (p(idx)); + + // The remainder determines the residue2s. The case of one pole is trivial. + if (lp == 1) + r = polyval (b, p); + return; + end + + // Determine the order of the denominator and remaining numerator. + // With the direct term removed, the potential order of the numerator + // is one less than the order of the denominator. + aorder = length (a) - 1; + border = aorder - 1; + + // Construct a system of equations relating the individual + // contributions from each residue2 to the complete numerator. + A = zeros (border+1, border+1); + + B = prepad (matrix (b, [length(b), 1]), border+1, 0); + B = B(:); // incase b have only 1 element + for ip = 1:length (p) + ri = zeros (size (p,1),size(p,2)); + ri(ip) = 1; + // A(:,ip) = prepad (rresidue (ri, p, [], tol), border+1, 0).';// invalid index error + temprr=rresidue (ri, p, [], tol) + ppr=prepad(temprr,border+1,0) + A(:,ip) = ppr + end + + // Solve for the residue2s. + // FIXME: Use a pre-conditioner d to make A \ B work better (bug #53869). + // It would be better to construct A and B so they are not close to + // singular in the first place. + d = max (abs (A),'c'); + + r = (diag (d) \ A) \ (B ./ d); + + endfunction + + // Reconstitute the numerator and denominator polynomials + // from the residue2s, poles, and direct term. + function [pnum, pden, e] = rresidue (r, p, k, tol, e ) + + if nargin < 3 then k = [] ; end + if nargin < 4 then tol = [] ; end + if nargin < 5 then e = [] ; end + if (~isempty (e)) + idx = 1:length (p); else - pden = conv (pden, pn); + [e, idx] = mpoles (p, tol, 0); + p = p(idx); + r = r(idx); end - end - - // D is the order of the denominator - // K is the order of the direct polynomial - // N is the order of the resulting numerator - // pnum(1:(N+1)) is the numerator's polynomial - // pden(1:(D+1)) is the denominator's polynomial - // pm is the multible pole for the nth residue - // pn is the numerator contribution for the nth residue - - D = length (pden) - 1; - K = length (k) - 1; - N = K + D; - pnum = zeros (1, N+1); - for n = indx(abs (rm) > 0) - p1 = [1, -p(n)]; - for m = 1:e(n) - if (m == 1) - pm = p1; + + idx = 1:length (p); + for n = idx + pn = [1, -p(n)]; + if (n == 1) + pden = pn; else - pm = conv (pm, p1); + pden = conv (pden, pn); end end - pn = deconv (real(pden),real(pm)); - pn = rm(n) * pn; - pnum = pnum + prepad (real(pn), N+1, 0, 2); - end - - // Add the direct term. - - if (length (k)) - pnum = pnum + conv (pden, k); - end - - // Check for leading zeros and trim the polynomial coefficients. - if (type(rm)==1 | type(p)==1 | type(k)==1) - small = max ([max(abs(pden)), max(abs(pnum)), 1]) * 1.1921e-07; - else - small = max ([max(abs(pden)), max(abs(pnum)), 1]) *%eps; - end - - pnum(abs (pnum) < small) = 0; - pden(abs (pden) < small) = 0; - - pnum = polyreduce (pnum); - pden = polyreduce (pden); - + + // D is the order of the denominator + // K is the order of the direct polynomial + // N is the order of the resulting numerator + // pnum(1:(N+1)) is the numerator's polynomial + // pden(1:(D+1)) is the denominator's polynomial + // pm is the multiple pole for the nth residue2 + // pn is the numerator contribution for the nth residue2 + + D = length (pden) - 1; + K = length (k) - 1; + N = K + D; + pnum = zeros (1, N+1); + for n = idx(abs (r) > 0) + p1 = [1, -p(n)]; + pn = 1; + for j = 1:n - 1 + pn = conv (pn, [1, -p(j)]); + end + for j = n + 1:length (p) + pn = conv (pn, [1, -p(j)]); + end + for j = 1:e(n) - 1 + pn = deconv (pn, p1); + end + pn = r(n) * pn; + pnum = pnum + prepad (pn, N+1, 0, 2); + end + + // Add the direct term. + if (length (k)) + pnum = pnum + conv (pden, k); + end + + pnum = polyreduce (pnum); + pden = polyreduce (pden); + + endfunction + +function y = polyreduce(p) + // Remove leading zeros from the polynomial + idx = find(p, 1) + if isempty(idx) then + y = 0; + else + y = p(idx(1):$); + end endfunction +/* +//test passed + b = [1, 1, 1]; + a = [1, -5, 8, -4]; + [r, p, k, e] = residue (b, a); + assert_checkalmostequal (r, [-2; 7; 3], 1e-12); + assert_checkalmostequal (p, [2; 2; 1], 1e-12); + assert_checkalmostequal (k,[]); + assert_checkalmostequal (e, [1; 2; 1]); + k = [1 0]; + b = conv (k, a) + prepad (b, length (k) + length (a) - 1, 0); + a = a; + [br, ar] = residue (r, p, k); + assert_checkalmostequal (br, b,1e-12); + assert_checkalmostequal (ar, a,1e-12); + [br, ar] = residue (r, p, k, e); + assert_checkalmostequal (br, b,1e-12); + assert_checkalmostequal (ar, a,1e-12); + +//test passed + b = [1, 0, 1]; + a = [1, 0, 18, 0, 81]; + [r, p, k, e] = residue (b, a); // error filter: Wrong type for input argument #1: Real matrix or polynomial expect + FIXME : builtin filter doesn't support complex paramters + r1 = [-5*%i; 12; +5*%i; 12]/54; + p1 = [+3*%i; +3*%i; -3*%i; -3*%i]; + assert_checkalmostequal (r, r1, 1e-12); + assert_checkalmostequal (p, p1, 1e-12); + assert_checkalmostequal (k,[]); + assert_checkalmostequal (e, [1; 2; 1; 2]); + [br, ar] = residue (r, p, k); + assert_checkalmostequal (br, b, 1e-12); + assert_checkalmostequal (ar, a, 1e-12); + +//test passed + r = [7; 3; -2]; + p = [2; 1; 2]; + k = [1 0]; + e = [2; 1; 1]; + [b, a] = residue (r, p, k, e); + assert_checkalmostequal (b, [1, -5, 9, -3, 1], 1e-12); + assert_checkalmostequal (a, [1, -5, 8, -4], 1e-12); + + + [rr, pr, kr, er] = residue (b, a); + [_, m] = mpoles (rr); + [_, n] = mpoles (r); + assert_checkalmostequal (rr(m), r(n), 1e-12); + assert_checkalmostequal (pr(m), p(n), 1e-12); + assert_checkalmostequal (kr, k, 1e-12); + assert_checkalmostequal (er(m), e(n), 1e-12); + +//test passed + b = [1]; + a = [1, 10, 25]; + [r, p, k, e] = residue (b, a); + r1 = [0; 1]; + p1 = [-5; -5]; + assert_checkalmostequal (r, r1, 1e-12); + assert_checkalmostequal (p, p1, 1e-12); + assert_checkalmostequal (k,[]); + assert_checkalmostequal (e, [1; 2]); + [br, ar] = residue (r, p, k); + assert_checkalmostequal (br, b, 1e-12); + assert_checkalmostequal (ar, a, 1e-12); + +// The following //test is due to Bernard Grung +//test <*34266> passed + z1 = 7.0372976777e6; + p1 = -3.1415926536e9; + p2 = -4.9964813512e8; + r1 = -(1 + z1/p1)/(1 - p1/p2)/p2/p1; + r2 = -(1 + z1/p2)/(1 - p2/p1)/p2/p1; + r3 = (1 + (p2 + p1)/p2/p1*z1)/p2/p1; + r4 = z1/p2/p1; + r = [r1; r2; r3; r4]; + p = [p1; p2; 0; 0]; + k = []; + e = [1; 1; 1; 2]; + b = [1, z1]; + a = [1, -(p1 + p2), p1*p2, 0, 0]; + [br, ar] = residue (r, p, k, e); + assert_checkalmostequal (br, [0,0,b],%eps,1e-8); + assert_checkalmostequal (ar, a, %eps,1e-8); + +//test <*49291> passed + rf = [1e3, 2e3, 1e3, 2e3]; + cf = [316.2e-9, 50e-9, 31.6e-9, 5e-9]; + [num, den] = residue (1./cf,-1./(rf.*cf),0); + assert_checkalmostequal (length (num), 4); + assert_checkalmostequal (length (den), 5); + assert_checkalmostequal (den(1), 1); + +//test <*51148> passed + r = [1.0000e+18, 3.5714e+12, 2.2222e+11, 2.1739e+10]; + pin = [-1.9231e+15, -1.6234e+09, -4.1152e+07, -1.8116e+06]; + k = 0; + [p, q] = residue (r, pin, k); + assert_checkalmostequal (p(4), 4.6828e+42, 1e-5); + +//test <*60384> passed + B = [1315.789473684211]; + A = [1, 1.100000536842105e+04, 1.703789473684211e+03, 0]; + poles1 = roots (A); + [r, p, k, e] = residue (B, A); + [B1, A1] = residue (r, p, k, e); + assert_checkalmostequal (B, B1); + assert_checkalmostequal (A, A1); + +%/ diff --git a/macros/sftrans.sci b/macros/sftrans.sci index 2dfd2bb..f0c092e 100644 --- a/macros/sftrans.sci +++ b/macros/sftrans.sci @@ -1,15 +1,14 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt // Original Source : https://octave.sourceforge.io/signal/ -// Modifieded by:Sonu Sharma, RGIT Mumbai +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - function [Sz, Sp, Sg] = sftrans (Sz, Sp, Sg, W, stop) //Transform band edges of a prototype filter (cutoff at W=1) represented in s-plane zero-pole-gain form (Frequency Transformation in Analog domain). @@ -71,7 +70,8 @@ function [Sz, Sp, Sg] = sftrans (Sz, Sp, Sg, W, stop) // Sz = // // 20. 10. 6.6666667 - + // dependencies + // funcprot(0); [nargout nargin]= argn(); @@ -115,12 +115,12 @@ function [Sz, Sp, Sg] = sftrans (Sz, Sp, Sg, W, stop) Sp = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)]; extend = [sqrt(-Fh*Fl), -sqrt(-Fh*Fl)]; if isempty(Sz) - Sz = [extend(1+rem([1:2*p],2))]; + Sz = [extend(1+modulo([1:2*p],2))]; else b = (C*(Fh-Fl)/2)./Sz; Sz = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)]; if (p > z) - Sz = [Sz, extend(1+rem([1:2*(p-z)],2))]; + Sz = [Sz, extend(1+modulo([1:2*(p-z)],2))]; end end else @@ -178,3 +178,13 @@ function [Sz, Sp, Sg] = sftrans (Sz, Sp, Sg, W, stop) end end endfunction +/** +Note : This function is tested with Octave's outputs as a reference. + +[Sz_new , Sp_new , Sg_new ] = sftrans([0.5;-0.5],[0.3 ; -0.3],2,0.5,0) // passed +[Sz_new , Sp_new , Sg_new ] = sftrans([0.5;-0.5],[0.3 ; -0.3],2,0.5,1) // passed +[Sz_new, Sp_new, Sg_new] = sftrans([0.5], [0.3], 1, [1.0, 2.0], 0) //passed +[Sz_new, Sp_new, Sg_new] = sftrans([0.5], [0.3], 1, [1.0, 2.0], 1) //passed +[Sz_new, Sp_new, Sg_new] = sftrans([], [], 1, 1.0, 0) //error passed +[Sz_new, Sp_new, Sg_new] = sftrans([0], [], 1, 1.0, 0)//error passed +*/
\ No newline at end of file diff --git a/macros/sos2cell.sci b/macros/sos2cell.sci index be40667..0434c6e 100644 --- a/macros/sos2cell.sci +++ b/macros/sos2cell.sci @@ -1,93 +1,69 @@ -function c = sos2cell(s,g) -//Converts a second order section matrix to a cell array -//Calling Sequences -//c=sos2cell(s) -//c=sos2cell(s,g) -//Parameters -//s -//An L-by-6 matrix where L is the number of sections -//g -//The scalar gain -//Description -//c=sos2cell(s) converts an L-by-6 second-order-section matrix s given by: -// s = [B1 A1 -// B2 A2 -// ... -// BL AL] -//to a cell array c = { {B1},{A1}, {B2},{A2}, ... {BL},{AL}} where each -//numerator vector Bi and denominator vector Ai contains the coefficients of a -//linear or quadratic polynomial. If the polynomial is linear, the coefficients -//zero-padded on the right -//c=sos2cell(s,g) adds a leading gain term to the start of the cell array as: -//c={ {[g,1]},{B1},{A1}, {B2},{A2}, ... {BL},{AL}} -//Example -//s=rand(2,6) -// s = -// -// -// column 1 to 5 -// -// 0.0437334 0.2639556 0.2806498 0.7783129 0.1121355 -// 0.4818509 0.4148104 0.1280058 0.2119030 0.6856896 -// -// column 6 -// -// 0.1531217 -// 0.6970851 -// -//sos2cell(s,2) -// ans = -// -// -// -// column 1 to 3 -// -//![2,1] [0.0437334,0.2639556,0.2806498] [0.7783129,0.1121355,0.1531217] ! -// -// column 4 to 5 -// -//![0.4818509,0.4148104,0.1280058] [0.2119030,0.6856896,0.6970851] ! -//Author -//Ankur Mallick - if(argn(2)<2) - g=[]; +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Author: Abinash Singh Under FOSSEE Internship +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +/* +Calling Sequence : + cll = sos2cell(s) + cll = sos2cell(s, g) +Description + sos2cell converts a second-order section matrix to a cell array representation. + The function can handle both unity-gain and non-unity gain filter systems. For non-unity gain systems, the gain factor is stored in the first cell of the output array. +Input Arguments + s - Second-order section matrix (L-by-6 matrix) + Each row represents one second-order section + Must have exactly 6 columns in format: [b0 b1 b2 a0 a1 a2] + Number of rows (L) represents the number of sections + g - Gain factor (optional) + Scalar value representing the system gain + Default value is 1 if not specified +Output Arguments + cll - Cell array containing second-order sections + For unity-gain systems (no gain specified): + Cell array with L elements + Each element contains coefficients: {[b0 b1 b2] [a0 a1 a2]} + For non-unity gain systems: + Cell array with L+1 elements + First element contains gain: {g 1} + Remaining elements contain section coefficients +*/ +function cll = sos2cell(s, g) + if (argn(2) > 2) then + error("sos2cell: Wrong number of input arguments"); + end + gain_inc = 1 ; + if nargin < 2 then + g = 1 ; + gain_inc = 0; + end + [L, n] = size(s); + if n ~= 6 then + error("sos2cell: Input matrix must have 6 columns"); end - if g==1 - g=[]; + if gain_inc then + L = L + 1 ; end - if(~or(type(s)==[1 5 8])|ndims(s)~=2|size(s,2)~=6) - error('Invalid Entry'); + cll = cell(1,L); + start_index=1 + if gain_inc then + cll{1} = { g 1}; + start_index = 2 ; end - L=size(s,1); - if ((L==1)&(~isempty(g))&(s==[1, 0, 0, 1, 0, 0])) - s=g*s; - g=[]; - end - c=cell(1,2*L); - k=0; - if(~isempty(g)) - c=cell(1,2*L+1); - c(1,1).entries=[g, 1]; - k=1; - end - for i=1:2:2*L - j=ceil(i/2); - sa=s(j,1:3); - ma=max(find(sa~=0)); - sb=s(j,4:6); - mb=max(find(sb~=0)); - cs=cell(1,2); - if(~isempty(ma)) - cs(1,1).entries=sa(1:ma); - else - cs(1,1).entries=[]; - end - if(~isempty(mb)) - cs(1,2).entries=sb(1:mb); - else - cs(1,2).entries=[]; - end - c(k+i)=cs(1,1); - c(k+i+1)=cs(1,2); + for i=start_index:L + cll{i}={s(i+1-start_index,1:3) s(i+1-start_index,4:6)} end endfunction +/* +sos = [3. 6. 7. 1. 1. 2. ; 1. 4. 5. 1. 9. 3. ; 2. 7. 1. 1. 7. 8.]; +cll = sos2cell(sos); // passed + +sos = [3. 6. 7. 1. 1. 2. ; 1. 4. 5. 1. 9. 3. ; 2. 7. 1. 1. 7. 8.]; +g = 0.5 ; +cll = sos2cell(sos,g); // passed + +*/
\ No newline at end of file diff --git a/macros/tf2zp.sci b/macros/tf2zp.sci index ca10fd4..f607ba0 100644 --- a/macros/tf2zp.sci +++ b/macros/tf2zp.sci @@ -1,13 +1,14 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -//Original contribution: FOSSEE, IIT Bombay +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - function [z,p,k]=tf2zp(num,den) // Transfer function to zero pole conversion @@ -60,8 +61,9 @@ function [z,p,k]=tf2zp(num,den) if(rd>1) then error("Denominator must be row vector"); - elseif np>cod then - error("Improper transfer function"); + + // elseif np>cod then + // error("Improper transfer function"); it is a bug it should be removed end if (~isempty(den)) then coef=den(1); @@ -98,3 +100,9 @@ function [z,p,k]=tf2zp(num,den) end z=roots(num); endfunction +/* +Note : This function is tested with Octave's outputs as a reference. +[z p k] = tf2zp([9.6500e+02 -1.1773e+05 3.8648e+06 -4.5127e+071.5581e+08],[ 1 -3 2]) //pass +[z p k] = tf2zp([1 2 3 4],[ 1 -3 2]) //pass +[z p k] = tf2zp([4 5 6],[1]) // pass +*/
\ No newline at end of file diff --git a/macros/zp2tf.sci b/macros/zp2tf.sci index da3e8b1..10da276 100644 --- a/macros/zp2tf.sci +++ b/macros/zp2tf.sci @@ -1,14 +1,14 @@ // Copyright (C) 2018 - IIT Bombay - FOSSEE -// // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Author:Sonu Sharma, RGIT Mumbai +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Date of Modification: 3 Feb 2024 // Organization: FOSSEE, IIT Bombay // Email: toolbox@scilab.in - function [num, den] = zp2tf (z, p, k) //Converts zeros / poles to a transfer function. @@ -54,3 +54,7 @@ function [num, den] = zp2tf (z, p, k) den = flipdim(den,2); endfunction +/* +[num,den] = zp2tf([1 3 4 5],[-4 -3 1 4],7) //passed +[num,den] = zp2tf([15 78 6 23],[2 1],965) //passed +*/
\ No newline at end of file |