summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--macros/diric.sci12
-rw-r--r--macros/filtic.sci44
-rw-r--r--macros/fractdiff.sci51
-rw-r--r--macros/gauspuls.sci29
-rw-r--r--macros/kaiserord.sci94
-rw-r--r--macros/meyeraux.sci11
-rw-r--r--macros/rceps.sci35
-rw-r--r--macros/rectpuls.sci19
-rw-r--r--macros/sawtooth.sci22
-rw-r--r--macros/sgolay.sci38
-rw-r--r--macros/sos2tf.sci43
-rw-r--r--macros/sos2zp.sci37
-rw-r--r--macros/welchwin.sci35
13 files changed, 323 insertions, 147 deletions
diff --git a/macros/diric.sci b/macros/diric.sci
index 077bb13..f60fee8 100644
--- a/macros/diric.sci
+++ b/macros/diric.sci
@@ -7,7 +7,7 @@ function [y]= diric(x,n)
// x: Real valued vector or matrix
// n: Real positive integer or complex integer
// Description
-// This is an Octave function
+
// y=diric(x,n) returns the dirichlet function values of parameter x.
// Examples
// 1. diric([1 2 3],3)
@@ -17,8 +17,14 @@ function [y]= diric(x,n)
funcprot(0);
rhs=argn(2);
-if (rhs~=2) then
+if (argn(2)~=2) then
error ("Wrong number of input arguments.")
-else y= callOctave("diric",x,n)
+elseif (n <= 0 | ceil(n) ~= n) then
+ error("n must be an positive integer.");
+
+ else
+ y = sin(n.*x./2)./(n.*sin(x./2));
+ y(pmodulo(x,2*%pi)==0) = (-1).^((n-1).*x(pmodulo(x,2*%pi)==0)./(2.*%pi));
end
+
endfunction
diff --git a/macros/filtic.sci b/macros/filtic.sci
index efd48ba..56367b7 100644
--- a/macros/filtic.sci
+++ b/macros/filtic.sci
@@ -1,5 +1,4 @@
-function zf = filtic (b, a, y, x)
-
+function zf = filtic(b,a,y,x)
//This function finds the initial conditions for the delays in the transposed direct-form II filter implementation
//Calling Sequence
//zf = filtic (b, a, y)
@@ -18,19 +17,36 @@ function zf = filtic (b, a, y, x)
// 0.00000 - 22.60000i
// 2.40000 + 0.00000i
// 0.00000 + 0.00000i
-//This function is being called from Octave
-funcprot(0);
-rhs = argn(2)
+ if (argn(2)>4 | argn(2)<3) | (argn(1)>1)
+ error("Wrong number of input agruments.")
+ end
+ if argn(2) < 4, x = []; end
+
+ nz = max(length(a)-1,length(b)-1);
+ zf=zeros(nz,1);
-if(rhs>4 | rhs<3)
- error("Wrong number of input agruments.")
-end
+
+ if length(a)<(nz+1)
+ a(length(a)+1:nz+1)=0;
+ end
+ if length(b)<(nz+1)
+ b(length(b)+1:nz+1)=0;
+ end
+
+ if length(x) < nz
+ x(length(x)+1:nz)=0;
+ end
+ if length(y) < nz
+ y(length(y)+1:nz)=0;
+ end
+
+ for i=nz:-1:1
+ for j=i:nz-1
+ zf(j) = b(j+1)*x(i) - a(j+1)*y(i)+zf(j+1);
+ end
+ zf(nz)=b(nz+1)*x(i)-a(nz+1)*y(i);
+ end
-select(rhs)
-case 3 then
-zf = callOctave("filtic",b,a,y)
-case 4 then
-zf = callOctave("filtic",b,a,y,x)
-end
endfunction
+
diff --git a/macros/fractdiff.sci b/macros/fractdiff.sci
index 979a079..fda26be 100644
--- a/macros/fractdiff.sci
+++ b/macros/fractdiff.sci
@@ -3,15 +3,44 @@ function y= fractdiff(x,d)
//Calling Sequence
// fractdiff (X, D)
//Description
-//This is an Octave function.
//Compute the fractional differences (1-L)^d x where L denotes the lag-operator and d is greater than -1.
- funcprot(0);
- rhs= argn(2);
- if(rhs < 2 | rhs >2)
- error("Wrong number of input arguments");
- end
- select(rhs)
- case 2 then
- y= callOctave("fractdiff",x,d);
- end
-endfunction \ No newline at end of file
+
+if(argn(2)~=2)
+ error("Wrong number of input arguments");
+end
+
+
+N = 100;
+
+ if (~ isvector (x) | ~isscalar(d))
+ error ("X must be a vector and D must be a scalar");
+ end
+
+ if (d<= -1) then
+ error ("D must be > -1");
+ end
+
+ if (d >= 1)
+ for k = 1 : d
+ x = x(2 : length (x)) - x(1 : length (x) - 1);
+ end
+ end
+//disp(x)
+ if (d >-1)
+ d=d-fix(d./1).*1;
+//disp(d)
+ if (d ~= 0)
+ n = (0 : N)';
+ w = real (gamma (-d+n) ./ gamma (-d) ./ gamma (n+1));
+ // disp(w); disp(x)
+ y = fftfilt (w, x);
+ y = y(1 : length (x));
+ else
+ // disp(x)
+ y = x;
+ end
+
+ end
+
+
+endfunction
diff --git a/macros/gauspuls.sci b/macros/gauspuls.sci
index ad853ee..621670b 100644
--- a/macros/gauspuls.sci
+++ b/macros/gauspuls.sci
@@ -1,5 +1,4 @@
function [y]=gauspuls(t,fc,bw)
-
// Generates Gaussian-modulated sinusoidal pulses
// Calling Sequence
// [y]=gauspuls(t,fc,bw)
@@ -10,7 +9,6 @@ function [y]=gauspuls(t,fc,bw)
// fc: Real non negative number or complex number
// bw: Real positive number or complex number
// Description
-// This is an Octave function
// This function returns a Gaussian RF pulse of unity amplitude at the times indicated in array t.
// Examples
// 1. gauspuls(1,2,3)
@@ -18,14 +16,21 @@ function [y]=gauspuls(t,fc,bw)
// 2. gauspuls([1 2 3],1,1)
// ans= 0.0281016 0.0000006 1.093D-14
-funcprot(0);
-rhs=argn(2);
-if ( rhs<1 ) then
- error ("Wrong number of input arguments.")
-elseif (rhs==1)
- y= callOctave("gauspuls",t)
-elseif (rhs==2)
- y= callOctave("gauspuls",t,fc)
-else y= callOctave("gauspuls",t,fc,bw)
-end
+if ( argn(2)<1 | argn(2)>3 ) then
+ error ("Wrong number of input arguments.")
+elseif (argn(2)==1)
+ fc = 1e3; bw = 0.5;
+elseif (argn(2)==2)
+ bw = 0.5;
+ end
+
+
+ if (~isscalar(fc) | ~isreal(fc) | fc<0) then
+ error('fc must be non-negative real scalar')
+ end
+ if (~isscalar(bw) | ~isreal(bw) | bw<=0) then
+ error('bw must be positive real scalar')
+ end
+ y = exp (-t .* t / (2*(1 / (4*%pi^2 * (-(bw^2 * fc^2) / (8 * log (10 ^ (-6/20)))))))) .* cos (2*%pi*fc * t);
+
endfunction
diff --git a/macros/kaiserord.sci b/macros/kaiserord.sci
index 47f4105..5f46193 100644
--- a/macros/kaiserord.sci
+++ b/macros/kaiserord.sci
@@ -10,7 +10,6 @@ function [n, Wn, beta, ftype] = kaiserord (f, m, dev, fs)
//dev: Deviation of the filter.
//fs: Sampling rate.
//Description
-//This is an Octave function.
//The vector f contains pairs of frequency band edges in the range [0,1]. The vector m specifies the magnitude response for each band. The values of m must be zero for all stop bands and must have the
//same magnitude for all pass bands. The deviation of the filter dev can be specified as a scalar or a vector of the same length as m. The optional sampling rate fs can be used to indicate that f is in
//Hz in the range [0,fs/2].
@@ -28,43 +27,70 @@ function [n, Wn, beta, ftype] = kaiserord (f, m, dev, fs)
// n = (A-8)/(2.285 dw)
//Examples
//[n, w, beta, ftype] = kaiserord ([1000, 1200], [1, 0], [0.05, 0.05], 11025)
-//n = 1
-//w = 1100
+//n = 70
+//w = 0.199
//beta = 1.5099
//ftype = low
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 3 | rhs > 4)
+if (argn(2) < 3 | argn(2) > 4)
error("Wrong number of input arguments.")
end
-select(rhs)
-
- case 3 then
- if(lhs==1)
- n = callOctave("kaiserord",f, m, dev)
- elseif(lhs==2)
- [n, Wn] = callOctave("kaiserord",f, m, dev)
- elseif(lhs==3)
- [n, Wn, beta] = callOctave("kaiserord",f, m, dev)
- elseif(lhs==4)
- [n, Wn, beta, ftype] = callOctave("kaiserord",f, m, dev)
- else
- error("Wrong number of output argments.")
- end
- case 4 then
- if(lhs==1)
- n = callOctave("kaiserord",f, m, dev, fs)
- elseif(lhs==2)
- [n, Wn] = callOctave("kaiserord",f, m, dev, fs)
- elseif(lhs==3)
- [n, Wn, beta] = callOctave("kaiserord",f, m, dev, fs)
- elseif(lhs==4)
- [n, Wn, beta, ftype] = callOctave("kaiserord",f, m, dev, fs)
- else
- error("Wrong number of output argments.")
- end
- end
+
+ if argn(2)<4 then
+ fs=2;
+ end
+
+ if length(f)~=(2*length(m)-2) then
+ error("kaiserord must have one magnitude for each frequency band");
+ end
+
+ if or(m(1:length(m)-2)~=m(3:length(m))) then
+ error("kaiserord pass and stop bands must be strictly alternating");
+ end
+ if length(dev)~=length(m) & length(dev)~=1 then
+ error("kaiserord must have one deviation for each frequency band");
+ end
+ dev = min(dev);
+ if dev <= 0 then
+ error("kaiserord must have dev>0");
+ end
+ Wn = (f(1:2:length(f))+f(2:2:length(f)))/fs;
+
+ if length(Wn) == 1 then
+ if m(1)>m(2) then
+ ftype='low';
+ else
+ ftype='high';
+ end
+ elseif length(w) == 2 then
+ if m(1)>m(2) then
+ ftype='stop';
+ else
+ ftype='pass';
+ end
+ else
+ if m(1)>m(2),
+ ftype='DC-1';
+ else
+ ftype='DC-0';
+ end
+ end
+
+ A = -20*log10(dev);
+ if (A > 50)
+ beta = 0.1102*(A-8.7);
+ elseif (A >= 21)
+ beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21);
+ else
+ beta = 0.0;
+ end
+
+ dw = 2*%pi*min(f(2:2:length(f))-f(1:2:length(f)))/fs;
+ n = max(1,ceil((A-8)/(2.285*dw)));
+
+ if ((m(1)>m(2)) == ((length(Wn)-fix(length(Wn)./2).*2)==0)) & (n-fix(n./2).*2)==1 then
+ n = n+1;
+ end
+
endfunction
diff --git a/macros/meyeraux.sci b/macros/meyeraux.sci
index 657e829..7cde1b9 100644
--- a/macros/meyeraux.sci
+++ b/macros/meyeraux.sci
@@ -1,12 +1,10 @@
function [y]=meyeraux(x)
-
// Returns value of Meyer Wavelet Auxiliary function
// Calling Sequence
// [y]=meyeraux(x)
// Parameters
// x: Real or complex valued vector or matrix
// Description
-// This is an Octave function.
// This function returns values of the auxiliary function used for Meyer wavelet generation.
// Examples
// 1. meyeraux([1 2 3])
@@ -14,10 +12,11 @@ function [y]=meyeraux(x)
// 2. meyeraux([1 2 3;4 5 6])
// ans= [1 -208 -10287 ; -118016 -709375 -2940624 ]
-funcprot(0);
-rhs=argn(2);
-if (rhs~=1) then
+if (argn(2)~=1) then
error ("Wrong number of input arguments.")
-else y=callOctave("meyeraux",x)
+else
+ y = 35.*x.^4-84.*x.^5+70.*x.^6-20.*x.^7;
end
endfunction
+
+
diff --git a/macros/rceps.sci b/macros/rceps.sci
index 8b4d89a..f0433a3 100644
--- a/macros/rceps.sci
+++ b/macros/rceps.sci
@@ -12,18 +12,37 @@ function [y, xm]= rceps(x)
// s(1:Fs/f0:length(s)) = 1; # Impulse glottal wave
// x = filter (1, a, s); # Speech signal in x
// [y, xm] = rceps (x .* hanning (1024)); # cepstrum and min phase reconstruction
-funcprot(0)
-lhs= argn(1)
-rhs= argn(2)
-if(rhs <1 | rhs> 1 )
+if(argn(2)~= 1 )
error("Wrong number of Input Arguments");
end
-if(lhs<2 | lhs>2)
+if(argn(1)>2)
error("Wrong number of Output Arguments")
end
-
- [y,xm]= callOctave("rceps",x);
-
+ f = abs(fft1(x));
+ if (or(f == 0))
+ error ("The spectrum of x contains zeros, unable to compute real cepstrum");
+ end
+
+ y = real(ifft1(log(f)));
+
+ if argn(1) == 2 then
+ n=length(x);
+ if size(x,1)==1 then
+ if (n-fix(n./2).*2) ==1 then
+ xm = [y(1), 2*y(2:n/2+1), zeros(1,n/2)];
+ else
+ xm = [y(1), 2*y(2:n/2), y(n/2+1), zeros(1,n/2-1)];
+ end
+ else
+ if (n-fix(n./2).*2)==1
+ xm = [y(1,:); 2*y(2:n/2+1,:); zeros(n/2,size(y,2))];
+ else
+ xm = [y(1,:); 2*y(2:n/2,:); y(n/2+1,:); zeros(n/2-1,size(y,2))];
+ end
+ end
+ xm = real(ifft1(exp(fft1(xm)')));
+ end
+
endfunction
diff --git a/macros/rectpuls.sci b/macros/rectpuls.sci
index 2796ef0..27138d0 100644
--- a/macros/rectpuls.sci
+++ b/macros/rectpuls.sci
@@ -1,5 +1,4 @@
function [y]=rectpuls(t,w)
-
// Generates a Rectangular pulse based on the width and sampling times.
// Calling Sequence
// [y]=rectpuls(t)
@@ -8,7 +7,6 @@ function [y]=rectpuls(t,w)
// t: Real or complex valued vector or matrix
// w: Real or complex valued vector or matrix
// Description
-// This is an Octave function
// y = rectpuls(t) returns a continuous, aperiodic, unity-height rectangular pulse depending upon input t, centered about t=0 and having default width of 1.
// y = rectpuls(t,w) generates a rectangle of width w.
// Examples
@@ -17,12 +15,17 @@ function [y]=rectpuls(t,w)
// 2. rectpuls([1000 1000 100 100])
// ans = 0 0 0 0
-funcprot(0);
-rhs=argn(2);
-if (rhs<1) then
+if (argn(2)<1 | argn(2)>2) then
error ("Wrong number of input arguments.")
-elseif (rhs==1)
- y=callOctave("rectpuls",t)
-else y=callOctave("rectpuls",t,w)
+elseif (argn(2)==1) then
+ w=1;
end
+
+ if (~ isreal (w) | ~ isscalar (w))
+ error ("W must be a real scalar");
+ end
+
+y = zeros (size(t,1), size(t,2));
+idx = find ((t >= -w/2) & (t < w/2));
+y(idx) = 1;
endfunction
diff --git a/macros/sawtooth.sci b/macros/sawtooth.sci
index 50d093d..8c3d689 100644
--- a/macros/sawtooth.sci
+++ b/macros/sawtooth.sci
@@ -8,7 +8,6 @@ function [y]=sawtooth (t,width)
// t: Real valued vector or matrix
// width: Real number between 0 and 1
// Description
-// This is an Octave function
// This function returns a sawtooth wave with period 2*pi with +1/-1 as the maximum and minimum values for elements of t. If width is specified, it determines where the maximum is in the interval [0,2*pi].
// Examples
// 1. sawtooth([1 2 3 4 5],0.5)
@@ -16,12 +15,21 @@ function [y]=sawtooth (t,width)
// 2. sawtooth([1 2; 4 5])
// ans = [-0.68169 -0.36338; 0.27324 0.59155]
-funcprot(0);
-rhs=argn(2);
-if (rhs<1) then
+if (argn(2)<1 | argn(2)>2) then
error ("Wrong number of input arguments.")
-elseif (rhs==1)
- y=callOctave("sawtooth",t)
-else y=callOctave("sawtooth",t,width)
+elseif (argn(2)==1)
+ width=1
end
+ if (width < 0 | width > 1 | ~ isreal (width))
+ error ("width must be a real value between 0 and 1");
+ end
+ t = pmodulo (t / (2 * %pi), 1);
+ y = zeros (size (t,1), size(t,2));
+
+ if (width ~= 0)
+ y (t < width) = 2 * t (t < width) / width - 1;
+ end
+ if (width ~= 1)
+ y( t >= width) = -2 * (t (t >= width) - width) / (1 - width) + 1;
+ end
endfunction
diff --git a/macros/sgolay.sci b/macros/sgolay.sci
index 1156d40..09062a7 100644
--- a/macros/sgolay.sci
+++ b/macros/sgolay.sci
@@ -10,7 +10,6 @@ function F = sgolay (p, n, m, ts)
//m: positive integer less than 2^31 or logical
//ts: real or complex value
//Description
-//This is an Octave function.
//This function computes the filter coefficients for all Savitzsky-Golay smoothing filters of order p for length n (odd).
//m can be used in order to get directly the mth derivative; ts is a scaling factor.
//Examples
@@ -20,21 +19,32 @@ function F = sgolay (p, n, m, ts)
// 0.33333 0.33333 0.33333
// -0.16667 0.33333 0.83333
-funcprot(0);
-rhs = argn(2)
-
-if(rhs<2 | rhs>4)
+if(argn(2)<2 | argn(2)>4) then
error("Wrong number of input arguments.")
+elseif ((n-fix(n./2).*2)~=1) then
+error ("sgolay needs an odd filter length n");
+elseif (p>=n) then
+ error ("sgolay needs filter length n larger than polynomial order p");
+ end
+
+if (argn(2)==2) then
+ m=0; ts=1;
+end
+if (argn(2)==3) then
+ ts=1;
end
- select(rhs)
- case 2 then
- F = callOctave("sgolay",p,n)
- case 3 then
- F = callOctave("sgolay",p,n,m)
- case 4 then
- F = callOctave("sgolay",p,n,m,ts)
- end
-endfunction
+if length(m) > 1, error("weight vector unimplemented"); end
+
+ F = zeros (n, n);
+ k = floor (n/2);
+ for row = 1:k+1
+ C = ( [(1:n)-row]'*ones(1,p+1) ) .^ ( ones(n,1)*[0:p] );
+ A = pinv(C);
+ F(row,:) = A(1+m,:);
+ end
+ F(k+2:n,:) = (-1)^m*F(k:-1:1,n:-1:1);
+ F = F * ( prod(1:m) / (ts^m) );
+endfunction
diff --git a/macros/sos2tf.sci b/macros/sos2tf.sci
index d38a4f7..f674c85 100644
--- a/macros/sos2tf.sci
+++ b/macros/sos2tf.sci
@@ -1,4 +1,4 @@
-function [B,A] = sos2tf(sos, g)
+function [A,B] = sos2tf(sos, g)
//This function converts series second-order sections to direct H(z) = B(z)/A(z) form.
//Calling Sequence
//[B] = sos2tf(sos)
@@ -8,7 +8,6 @@ function [B,A] = sos2tf(sos, g)
//sos: matrix of real or complex numbers
//g: real or complex value, default value is 1
//Description
-//This is an Octave function.
//This function converts series second-order sections to direct H(z) = B(z)/A(z) form.
//The input is the sos matrix and the second parameter is the overall gain, default value of which is 1.
//The output is a vector.
@@ -19,18 +18,40 @@ function [B,A] = sos2tf(sos, g)
// -2 1 2 4 1
//b =
// 1 10 0 -10 -1
-funcprot(0);
-rhs = argn(2)
-if(rhs<1 | rhs>2)
+if(argn(2)<1 | argn(2)>2)
error("Wrong number of input arguments.")
end
+if argn(2)==1 then
+ g=1;
+end
+
+ [N,M] = size(sos);
+
+ if M~=6
+ error('sos matrix must be N by 6');
+ end
+
+ A = 1;
+ B = 1;
+
+ for i=1:N
+ A = conv(A, sos(i,1:3));
+ B = conv(B, sos(i,4:6));
+ end
+
+ nB = length(A);
+
+ while nB & A(nB)==0
+ A=A(1:nB-1);
+ nB=length(A);
+ end
+ nA = length(B);
+ while nA & B(nA)==0
+ B=B(1:nA-1);
+ nA=length(B);
+ end
+ A = A * g;
- select(rhs)
- case 1 then
- [B,A] = callOctave("sos2tf",sos)
- case 2 then
- [B,A] = callOctave("sos2tf",sos,g)
- end
endfunction
diff --git a/macros/sos2zp.sci b/macros/sos2zp.sci
index 98abc36..1f0f335 100644
--- a/macros/sos2zp.sci
+++ b/macros/sos2zp.sci
@@ -11,7 +11,6 @@ function [z,p,k] = sos2zp (sos, g)
//z: column vector
//p: column vector
//Description
-//This is an Octave function.
//This function converts series second-order sections to zeros, poles, and gains (pole residues).
//The input is the sos matrix and the second parameter is the overall gain, default value of which is 1.
//The outputs are z, p, k. z and p are column vectors containing zeros and poles respectively, and k is the overall gain.
@@ -24,18 +23,32 @@ function [z,p,k] = sos2zp (sos, g)
// -0.6250 + 1.0533i
// -0.6250 - 1.0533i
//c = 1
-
-funcprot(0);
-rhs = argn(2)
-if(rhs<1 | rhs>2)
+if(argn(2)<1 | argn(2)>2)
error("Wrong number of input arguments.")
end
+if argn(2)==1 then
+ g=1;
+end
+ gns = sos(:,1);
+ k = prod(gns)*g;
+ if k==0 then
+ error('one or more section gains is zero');
+ end
+ sos(:,1:3) = sos(:,1:3)./ [gns gns gns];
- select(rhs)
- case 1 then
- [z,p,k] = callOctave("sos2zp",sos)
- case 2 then
- [z,p,k] = callOctave("sos2zp",sos,g)
- end
-endfunction
+ [N,m] = size(sos);
+ if m~=6 then
+ error('sos matrix should be N by 6');
+ end
+ z = zeros(2*N,1);
+ p = zeros(2*N,1);
+ for i=1:N
+ ndx = [2*i-1:2*i];
+ zi = roots(sos(i,1:3));
+ z(ndx) = zi;
+ pi = roots(sos(i,4:6));
+ p(ndx) = pi;
+ end
+
+endfunction
diff --git a/macros/welchwin.sci b/macros/welchwin.sci
index a871a59..9b5fac1 100644
--- a/macros/welchwin.sci
+++ b/macros/welchwin.sci
@@ -23,15 +23,36 @@ function w = welchwin (m, opt)
rhs = argn(2)
-if(rhs<1 | rhs>2)
+if(argn(2)<1 | argn(2)>2)
error("Wrong number of input arguments.")
end
- select(rhs)
- case 1 then
- w = callOctave("welchwin",m)
- case 2 then
- w = callOctave("welchwin",m,opt)
- end
+ if ((~isscalar (m) & (m == fix (m)) & (m > 0))) then
+ error (" M must be a positive integer");
+ end
+
+
+ N = (m - 1) / 2;
+ mmin = 3;
+if argn(2)==2 then
+
+ select (opt)
+ case "periodic"
+ N = m / 2;
+ mmin = 2;
+ case "symmetric"
+ N = (m - 1) / 2;
+ else
+ error ("window type must be either periodic or symmetric");
+ end
+ end
+
+
+ if (m < mmin)
+ error ('M must be an integer greater than mmin');
+ end
+
+ n = [0:m-1]';
+ w = 1 - ((n-N)./N).^2;
endfunction
generated by cgit (git 2.34.1) at 2024-12-19 12:58:47 +0000