Winter Internship 2024 work done by Abinash Singh

2 files changed, 296 insertions, 88 deletions

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/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