diff options
| -rw-r--r-- | macros/ellip.sci | 64 | ||||
| -rw-r--r-- | macros/findpeaks.sci | 380 | ||||
| -rw-r--r-- | macros/freqz.sci | 316 | ||||
| -rw-r--r-- | macros/iirlp2mb.sci | 378 | ||||
| -rw-r--r-- | macros/parser.sci | 82 |
5 files changed, 1050 insertions, 170 deletions
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/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/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/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 |