diff options
Diffstat (limited to 'macros')
| -rw-r--r-- | macros/arch_rnd.sci | 43 | ||||
| -rw-r--r-- | macros/arma_rnd.sci | 52 | ||||
| -rw-r--r-- | macros/cl2bp.sci | 47 | ||||
| -rw-r--r-- | macros/cplxreal.sci | 5 | ||||
| -rw-r--r-- | macros/fft1.sci (renamed from macros/fft.sci) | 4 | ||||
| -rw-r--r-- | macros/fft21.sci (renamed from macros/fft2.sci) | 4 | ||||
| -rw-r--r-- | macros/fftshift1.sci | 8 | ||||
| -rw-r--r-- | macros/fftw1.sci | 64 | ||||
| -rw-r--r-- | macros/ifft1.sci (renamed from macros/ifft.sci) | 4 | ||||
| -rw-r--r-- | macros/invfreqs.sci | 96 | ||||
| -rw-r--r-- | macros/invfreqz.sci | 94 | ||||
| -rw-r--r-- | macros/kaiserord.sci | 70 | ||||
| -rw-r--r-- | macros/lib | bin | 6868 -> 6614 bytes | |||
| -rw-r--r-- | macros/names | 21 | ||||
| -rw-r--r-- | macros/pburg.sci | 93 | ||||
| -rw-r--r-- | macros/pyulear.sci | 85 | ||||
| -rw-r--r-- | macros/qp_kaiser.sci | 37 | ||||
| -rw-r--r-- | macros/sftrans.sci | 74 | ||||
| -rw-r--r-- | macros/tfestimate.sci | 93 | ||||
| -rw-r--r-- | macros/xcorr1.sci | 86 |
20 files changed, 966 insertions, 14 deletions
diff --git a/macros/arch_rnd.sci b/macros/arch_rnd.sci new file mode 100644 index 0000000..50facc4 --- /dev/null +++ b/macros/arch_rnd.sci @@ -0,0 +1,43 @@ +function res = arch_rnd (a, b, t) +//Simulate an ARCH sequence of length t with AR coefficients b and CH coefficients a. +//Calling Sequence +//arch_rnd (a, b, t) +//Parameters +//a: CH coefficients +//b: AR coefficients +//t: Length of ARCH sequence +//Description +//This is an Octave function. +//It Simulates an ARCH sequence of length t with AR coefficients b and CH coefficients a. +//The result y(t) follows the model +// +//y(t) = b(1) + b(2) * y(t-1) + … + b(lb) * y(t-lb+1) + e(t), +//where e(t), given y up to time t-1, is N(0, h(t)), with +// +//h(t) = a(1) + a(2) * e(t-1)^2 + … + a(la) * e(t-la+1)^2 +//Examples +//a = [1 2 3 4 5]; +//b = [7 8 9 10]; +//arch_rnd (a, b, t) +//ans = +// +// 6.1037e+00 +// 5.7294e+01 +// 5.7390e+02 +// 6.3063e+03 +// 6.8695e+04 + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 3 | rhs > 3) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 3 then + res = callOctave("arch_rnd",a, b, t) + + end +endfunction diff --git a/macros/arma_rnd.sci b/macros/arma_rnd.sci new file mode 100644 index 0000000..bc826f4 --- /dev/null +++ b/macros/arma_rnd.sci @@ -0,0 +1,52 @@ +function res = arma_rnd (a, b, v, t, n) +//Return a simulation of the ARMA model. +//Calling Sequence +//arma_rnd (a, b, v, t, n) +//arma_rnd (a, b, v, t) +//Parameters +//a: vector +//b: vector +//v: Variance +//t: Length of output vector +//n: Number of dummy x(i) used for initialization +//Description +//This is an Octave function. +//The ARMA model is defined by +// +//x(n) = a(1) * x(n-1) + … + a(k) * x(n-k) +// + e(n) + b(1) * e(n-1) + … + b(l) * e(n-l) +//in which k is the length of vector a, l is the length of vector b and e is Gaussian white noise with variance v. The function returns a vector of length t. +// +//The optional parameter n gives the number of dummy x(i) used for initialization, i.e., a sequence of length t+n is generated and x(n+1:t+n) is returned. If n is omitted, n = 100 is used. +//Examples +//a = [1 2 3 4 5]; +//b = [7; 8; 9; 10; 11]; +//v = 10; +//t = 5; +//n = 100; +//arma_rnd (a, b, v, t, n) +//ans = +// +// -1.6176e+05 +// -4.1663e+05 +// -1.0732e+06 +// -2.7648e+06 +// -7.1221e+06 + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 5 | rhs > 6) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 5 then + res = callOctave("arma_rnd",a, b, v, t) + + case 6 then + res = callOctave("arma_rnd",a, b, v, t, n) + + end +endfunction diff --git a/macros/cl2bp.sci b/macros/cl2bp.sci new file mode 100644 index 0000000..1069cf0 --- /dev/null +++ b/macros/cl2bp.sci @@ -0,0 +1,47 @@ +function h = cl2bp (m, w1, w2, up, lo, gridsize) +//Constrained L2 bandpass FIR filter design. +//Calling Sequence +//h = cl2bp (m, w1, w2, up, lo, gridsize) +//h = cl2bp (m, w1, w2, up, lo) +//Parameters +//m: degree of cosine polynomial, i.e. the number of output coefficients will be m*2+1 +//w1 and w2: bandpass filter cutoffs in the range 0 <= w1 < w2 <= pi, where pi is the Nyquist frequency +//up: vector of 3 upper bounds for [stopband1, passband, stopband2] +//lo: vector of 3 lower bounds for [stopband1, passband, stopband2] +//gridsize: search grid size; larger values may improve accuracy, but greatly increase calculation time. +//Description +//This is an Octave function. +//Constrained L2 bandpass FIR filter design. Compared to remez, it offers implicit specification of transition bands, a higher likelihood of convergence, and an error criterion combining features of both L2 and Chebyshev approaches. +//Examples +//h = cl2bp(5, 0.3*pi, 0.6*pi, [0.02, 1.02, 0.02], [-0.02, 0.98, -0.02], 2^11) +//h = +// +// 0.038311 +// 0.082289 +// -0.086163 +// -0.226006 +// 0.047851 +// 0.307434 +// 0.047851 +// -0.226006 +// -0.086163 +// 0.082289 +// 0.038311 + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 5 | rhs > 6) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 5 then + res = callOctave("cl2bp", m, w1, w2, up, lo) + + case 6 then + res = callOctave("cl2bp", m, w1, w2, up, lo, gridsize) + + end +endfunction diff --git a/macros/cplxreal.sci b/macros/cplxreal.sci index 8443837..a552998 100644 --- a/macros/cplxreal.sci +++ b/macros/cplxreal.sci @@ -14,6 +14,11 @@ function [zc, zr] = cplxreal (z, thresh) //This is an Octave function. //Every complex element of z is expected to have a complex-conjugate elsewhere in z. From the pair of complex-conjugates, the one with the negative imaginary part is removed. //If the magnitude of the imaginary part of an element is less than the thresh, it is declared as real. +//Examples +//[zc, zr] = cplxreal([1 2 3+i 4 3-i 5]) +//zc = 3 + 1i +//zr = +// 1 2 4 5 funcprot(0); lhs = argn(1) rhs = argn(2) diff --git a/macros/fft.sci b/macros/fft1.sci index 0ea54ab..a0a6438 100644 --- a/macros/fft.sci +++ b/macros/fft1.sci @@ -1,4 +1,4 @@ -function res = fft (x, n, dim) +function res = fft1 (x, n, dim) //Calculates the discrete Fourier transform of a matrix using Fast Fourier Transform algorithm. //Calling Sequence //fft (x, n, dim) @@ -20,7 +20,7 @@ function res = fft (x, n, dim) //x = [1 2 3; 4 5 6; 7 8 9] //n = 3 //dim = 2 -//fft (x, n, dim) +//fft1 (x, n, dim) //ans = // // 6.0000 + 0.0000i -1.5000 + 0.8660i -1.5000 - 0.8660i diff --git a/macros/fft2.sci b/macros/fft21.sci index 58cd70d..2f04200 100644 --- a/macros/fft2.sci +++ b/macros/fft21.sci @@ -1,4 +1,4 @@ -function res = fft2 (A, m, n) +function res = fft21 (A, m, n) //Calculates the two-dimensional discrete Fourier transform of A using a Fast Fourier Transform algorithm. //Calling Sequence //fft2 (A, m, n) @@ -15,7 +15,7 @@ function res = fft2 (A, m, n) //x = [1 2 3; 4 5 6; 7 8 9] //m = 4 //n = 4 -//fft2 (A, m, n) +//fft21 (A, m, n) //ans = // // 45 + 0i -6 - 15i 15 + 0i -6 + 15i diff --git a/macros/fftshift1.sci b/macros/fftshift1.sci index d8721cd..9026025 100644 --- a/macros/fftshift1.sci +++ b/macros/fftshift1.sci @@ -1,5 +1,5 @@ function y= fftshift1(X,DIM) -//Perform a shift of the vector X, for use with the 'fft' and 'ifft' functions, in order the move the frequency 0 to the center of the vector or matrix. +//Perform a shift of the vector X, for use with the 'fft1' and 'ifft1' functions, in order the move the frequency 0 to the center of the vector or matrix. //Calling Sequence // fftshift1 (X) // fftshift1 (X, DIM) @@ -8,9 +8,9 @@ function y= fftshift1(X,DIM) //DIM: The optional DIM argument can be used to limit the dimension along which the permutation occurs //Description //This is an Octave function. -//Perform a shift of the vector X, for use with the 'fft' and 'ifft' functions, in order the move the frequency 0 to the center of the vector or matrix. +//Perform a shift of the vector X, for use with the 'fft1' and 'ifft1' functions, in order the move the frequency 0 to the center of the vector or matrix. // -//If X is a vector of N elements corresponding to N time samples spaced by dt, then 'fftshift1 (fft (X))' corresponds to frequencies +//If X is a vector of N elements corresponding to N time samples spaced by dt, then 'fftshift1 (fft1 (X))' corresponds to frequencies // //f = [ -(ceil((N-1)/2):-1:1)*df 0 (1:floor((N-1)/2))*df ] // @@ -30,4 +30,4 @@ function y= fftshift1(X,DIM) case 2 then y=callOctave("fftshift",X,DIM); end -endfunction
\ No newline at end of file +endfunction diff --git a/macros/fftw1.sci b/macros/fftw1.sci new file mode 100644 index 0000000..a530017 --- /dev/null +++ b/macros/fftw1.sci @@ -0,0 +1,64 @@ +function res = fftw1(a, b) +//Manage FFTW wisdom data. +//Calling Sequence +// method = fftw ("planner") +//fftw ("planner", method) +//wisdom = fftw ("dwisdom") +//fftw ("dwisdom", wisdom) +//fftw ("threads", nthreads) +//nthreads = fftw ("threads") +//Parameters +//Description +//This is an Octave function. +//Wisdom data can be used to significantly accelerate the calculation of the FFTs, but implies an initial cost in its calculation. When the FFTW libraries are initialized, they read a system wide wisdom +//file (typically in /etc/fftw/wisdom), allowing wisdom to be shared between applications other than Octave. Alternatively, the fftw function can be used to import wisdom. For example, +// +//wisdom = fftw ("dwisdom") +//will save the existing wisdom used by Octave to the string wisdom. This string can then be saved to a file and restored using the save and load commands respectively. This existing wisdom can be re +//imported as follows +// +//fftw ("dwisdom", wisdom) +//If wisdom is an empty string, then the wisdom used is cleared. +// +//During the calculation of Fourier transforms further wisdom is generated. The fashion in which this wisdom is generated is also controlled by the fftw function. There are five different manners in which +//the wisdom can be treated: +// +//"estimate" +//Specifies that no run-time measurement of the optimal means of calculating a particular is performed, and a simple heuristic is used to pick a (probably sub-optimal) plan. The advantage of this method +//is that there is little or no overhead in the generation of the plan, which is appropriate for a Fourier transform that will be calculated once. +// +//"measure" +//In this case a range of algorithms to perform the transform is considered and the best is selected based on their execution time. +// +//"patient" +//Similar to "measure", but a wider range of algorithms is considered. +// +//"exhaustive" +//Like "measure", but all possible algorithms that may be used to treat the transform are considered. +// +//"hybrid" +//As run-time measurement of the algorithm can be expensive, this is a compromise where "measure" is used for transforms up to the size of 8192 and beyond that the "estimate" method is used. +// +//The default method is "estimate". The current method can be queried with +// +//method = fftw ("planner") +//or set by using +// +//fftw ("planner", method) + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 1 | rhs > 2) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 1 then + res = callOctave("fftw",a) + + case 2 then + res = callOctave("fftw",a, b) + end +endfunction diff --git a/macros/ifft.sci b/macros/ifft1.sci index 70688e4..3e94934 100644 --- a/macros/ifft.sci +++ b/macros/ifft1.sci @@ -1,4 +1,4 @@ -function res = ifft (x, n, dim) +function res = ifft1 (x, n, dim) //Calculates the inverse discrete Fourier transform of a matrix using Fast Fourier Transform algorithm. //Calling Sequence //ifft (x, n, dim) @@ -22,7 +22,7 @@ function res = ifft (x, n, dim) //x = [1 2 3; 4 5 6; 7 8 9] //n = 3 //dim = 2 -//ifft (x, n, dim) +//ifft1 (x, n, dim) //ans = // // 2.00000 + 0.00000i -0.50000 - 0.28868i -0.50000 + 0.28868i diff --git a/macros/invfreqs.sci b/macros/invfreqs.sci new file mode 100644 index 0000000..32e732b --- /dev/null +++ b/macros/invfreqs.sci @@ -0,0 +1,96 @@ +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 +endfunction + diff --git a/macros/invfreqz.sci b/macros/invfreqz.sci new file mode 100644 index 0000000..d051499 --- /dev/null +++ b/macros/invfreqz.sci @@ -0,0 +1,94 @@ +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 + +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("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 + + 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 + diff --git a/macros/kaiserord.sci b/macros/kaiserord.sci new file mode 100644 index 0000000..47f4105 --- /dev/null +++ b/macros/kaiserord.sci @@ -0,0 +1,70 @@ +function [n, Wn, beta, ftype] = kaiserord (f, m, dev, fs) +//Return the parameters needed to produce a filter of the desired specification from a Kaiser window. +//Calling Sequence +//[n, Wn, beta, ftype] = kaiserord (f, m, dev, fs) +//[…] = kaiserord (f, m, dev, fs) +//[…] = kaiserord (f, m, dev) +//Parameters +//f: Pairs of frequency band edges. +//m: Magnitude response for each band. +//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]. +// +//The returned value n is the required order of the filter (the length of the filter minus 1). The vector Wn contains the band edges of the filter suitable for passing to fir1. The value beta is the +//parameter of the Kaiser window of length n+1 to shape the filter. The string ftype contains the type of filter to specify to fir1. +// +//The Kaiser window parameters n and beta are computed from the relation between ripple (A=-20*log10(dev)) and transition width (dw in radians) discovered empirically by Kaiser: +// +// +// / 0.1102(A-8.7) A > 50 +// beta = | 0.5842(A-21)^0.4 + 0.07886(A-21) 21 <= A <= 50 +// \ 0.0 A < 21 +// +// 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 +//beta = 1.5099 +//ftype = low + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 3 | rhs > 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 +endfunction Binary files differdiff --git a/macros/names b/macros/names index 001dba2..800df1a 100644 --- a/macros/names +++ b/macros/names @@ -3,9 +3,11 @@ ac2rc arParEst ar_psd arburg +arch_rnd arch_test arch_fit arcov +arma_rnd armcov aryule autoreg_matrix @@ -37,6 +39,7 @@ cheby1 cheby2 check chirp +cl2bp clustersegment cmorwavf cohere @@ -66,10 +69,11 @@ ellipord enbw eqtflength falltime -fft -fft2 +fft1 +fft21 fftconv fftfilt +fftw1 fftn fftshift1 fht @@ -109,8 +113,8 @@ icceps idct1 idct2 idst1 -ifft -ifft2 +ifft1 +ifft21 ifftn ifftshift1 ifht @@ -122,6 +126,8 @@ impzlength interp intfilt invfreq +invfreqs +invfreqz invimpinvar is2rc isallpass @@ -131,6 +137,7 @@ ismaxphase isminphase isstable kaiser +kaiserord lar2rc latc2tf latcfilt @@ -155,6 +162,7 @@ ncauer nnls nuttallwin parzenwin +pburg pchip pchips peak2peak @@ -180,7 +188,9 @@ pulseperiod pulsesep pulsewidth pulstran +pyulear pwelch +qp_kaiser rc2ac rc2is rc2lar @@ -204,6 +214,7 @@ sawtooth schtrig schurrc seqperiod +sftrans sgolay sgolayfilt shanwavf @@ -232,6 +243,7 @@ synthesis tf2sos tf2zp tf2zpk +tfestimate tfe transpose trial_iirlp2mb @@ -254,6 +266,7 @@ welchwin window wkeep wrev +xcorr1 xcorr2 xcov1 yulewalker diff --git a/macros/pburg.sci b/macros/pburg.sci new file mode 100644 index 0000000..78198ed --- /dev/null +++ b/macros/pburg.sci @@ -0,0 +1,93 @@ +function [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion) +//Calculate Burg maximum-entropy power spectral density. +//Calling Sequence +//[psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion) +//All but the first two arguments are optional and may be empty. +//Parameters +// x: [vector] sampled data +// poles: [integer scalar] required number of poles of the AR model +// freq: [real vector] frequencies at which power spectral density is calculated [integer scalar] number of uniformly distributed frequency values at which spectral density is calculated. [default=256] +// Fs: [real scalar] sampling frequency (Hertz) [default=1] +// range: 'half', 'onesided' : frequency range of the spectrum is from zero up to but not including sample_f/2. Power from negative frequencies is added to the positive side of the spectrum. 'whole', 'twosided' : frequency range of the spectrum is -sample_f/2 to sample_f/2, with negative frequencies stored in "wrap around" order after the positive frequencies; e.g. frequencies for a 10-point 'twosided' spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1 'shift', 'centerdc' : same as 'whole' but with the first half of the spectrum swapped with second half to put the zero-frequency value in the middle. (See "help fftshift". If "freq" is vector, 'shift' is ignored. If model coefficients "ar_coeffs" are real, the default range is 'half', otherwise default range is 'whole'. +// method: 'fft': use FFT to calculate power spectral density. 'poly': calculate spectral density as a polynomial of 1/z N.B. this argument is ignored if the "freq" argument is a vector. The default is 'poly' unless the "freq" argument is an integer power of 2. +// plot_type: 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db' specifies the type of plot. The default is 'plot', which means linear-linear axes. 'squared' is the same as 'plot'. 'dB' plots "10*log10(psd)". This argument is ignored and a spectrum is not plotted if the caller requires a returned value. +// criterion: [optional string arg] model-selection criterion. Limits the number of poles so that spurious poles are not added when the whitened data has no more information in it (see Kay & Marple, 1981). Recognized values are 'AKICc' -- approximate corrected Kullback information criterion (recommended), 'KIC' -- Kullback information criterion 'AICc' -- corrected Akaike information criterion 'AIC' -- Akaike information criterion 'FPE' -- final prediction error" criterion The default is to NOT use a model-selection criterion. +//Description +//This function is being called from Octave +//This function is a wrapper for arburg and ar_psd. +//The functions "arburg" and "ar_psd" do all the work. +//See "help arburg" and "help ar_psd" for further details. +//Examples +//a = [1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746]; +//[psd,f_out] = pburg(a,2); + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 2 | rhs > 8) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 2 then + if(lhs==1) + psd = callOctave("pburg",x,poless) + elseif(lhs==2) + [psd,f_out] = callOctave("pburg",x,poles) + else + error("Wrong number of output argments.") + end + + case 3 then + if(lhs==1) + psd = callOctave("pburg",x,poles,freq) + elseif(lhs==2) + [psd,f_out] = callOctave("pburg",x,poles,freq) + else + error("Wrong number of output argments.") + end + case 4 then + if(lhs==1) + psd = callOctave("pburg",x,poles,freq,Fs) + elseif(lhs==2) + [psd,f_out] = callOctave("pburg",x,poles,freq,Fs) + else + error("Wrong number of output argments.") + end + case 5 then + if(lhs==1) + psd = callOctave("pburg",x,poles,freq,Fs,range) + elseif(lhs==2) + [psd,f_out] = callOctave("pburg",x,poles,freq,Fs,range) + else + error("Wrong number of output argments.") + end + case 6 then + if(lhs==1) + psd = callOctave("pburg",x,poles,freq,Fs,range,method) + elseif(lhs==2) + [psd,f_out] = callOctave("pburg",x,poles,freq,Fs,range,method) + else + error("Wrong number of output argments.") + end + case 7 then + if(lhs==1) + psd = callOctave("pburg",x,poles,freq,Fs,range,method,plot_type) + elseif(lhs==2) + [psd,f_out] = callOctave("pburg",x,poles,freq,Fs,range,method,plot_type) + else + error("Wrong number of output argments.") + end + case 8 then + if(lhs==1) + psd = callOctave("pburg",x,poles,freq,Fs,range,method,plot_type,criterion) + elseif(lhs==2) + [psd,f_out] = callOctave("pburg",x,poles,freq,Fs,range,method,plot_type,criterion) + else + error("Wrong number of output argments.") + end + end +endfunction + + diff --git a/macros/pyulear.sci b/macros/pyulear.sci new file mode 100644 index 0000000..55278d8 --- /dev/null +++ b/macros/pyulear.sci @@ -0,0 +1,85 @@ +function [psd,f_out] = pyulear(x,poles,freq,Fs,range,method,plot_type) + +//Calculates a Yule-Walker autoregressive (all-pole) model of the data "x" and computes the power spectrum of the model. +//Calling Sequence +//[psd,f_out] = pyulear(x,poles,freq,Fs,range,method,plot_type) +//All but the first two arguments are optional and may be empty. +//Parameters +// x: [vector] sampled data +// poles: [integer scalar] required number of poles of the AR model +// freq: [real vector] frequencies at which power spectral density is calculated [integer scalar] number of uniformly distributed frequency values at which spectral density is calculated. [default=256] +// Fs: [real scalar] sampling frequency (Hertz) [default=1] +// range: 'half', 'onesided' : frequency range of the spectrum is from zero up to but not including sample_f/2. Power from negative frequencies is added to the positive side of the spectrum. 'whole', 'twosided' : frequency range of the spectrum is -sample_f/2 to sample_f/2, with negative frequencies stored in "wrap around" order after the positive frequencies; e.g. frequencies for a 10-point 'twosided' spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1 'shift', 'centerdc' : same as 'whole' but with the first half of the spectrum swapped with second half to put the zero-frequency value in the middle. (See "help fftshift". If "freq" is vector, 'shift' is ignored. If model coefficients "ar_coeffs" are real, the default range is 'half', otherwise default range is 'whole'. +// method: 'fft': use FFT to calculate power spectral density. 'poly': calculate spectral density as a polynomial of 1/z N.B. this argument is ignored if the "freq" argument is a vector. The default is 'poly' unless the "freq" argument is an integer power of 2. +// plot_type: 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db' specifies the type of plot. The default is 'plot', which means linear-linear axes. 'squared' is the same as 'plot'. 'dB' plots "10*log10(psd)". This argument is ignored and a spectrum is not plotted if the caller requires a returned value. + +//Description +//This function is being called from Octave. +//This function is a wrapper for aryule and ar_psd. +//See "help aryule" and "help ar_psd" for further details. +//Examples +//a = [1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746]; +//[psd,f_out] = pyulear(a,2); + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 2 | rhs > 7) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 2 then + if(lhs==1) + psd = callOctave("pyulear",x,poless) + elseif(lhs==2) + [psd,f_out] = callOctave("pyulear",x,poles) + else + error("Wrong number of output argments.") + end + + case 3 then + if(lhs==1) + psd = callOctave("pyulear",x,poles,freq) + elseif(lhs==2) + [psd,f_out] = callOctave("pyulear",x,poles,freq) + else + error("Wrong number of output argments.") + end + case 4 then + if(lhs==1) + psd = callOctave("pyulear",x,poles,freq,Fs) + elseif(lhs==2) + [psd,f_out] = callOctave("pyulear",x,poles,freq,Fs) + else + error("Wrong number of output argments.") + end + case 5 then + if(lhs==1) + psd = callOctave("pyulear",x,poles,freq,Fs,range) + elseif(lhs==2) + [psd,f_out] = callOctave("pyulear",x,poles,freq,Fs,range) + else + error("Wrong number of output argments.") + end + case 6 then + if(lhs==1) + psd = callOctave("pyulear",x,poles,freq,Fs,range,method) + elseif(lhs==2) + [psd,f_out] = callOctave("pyulear",x,poles,freq,Fs,range,method) + else + error("Wrong number of output argments.") + end + case 7 then + if(lhs==1) + psd = callOctave("pyulear",x,poles,freq,Fs,range,method,plot_type) + elseif(lhs==2) + [psd,f_out] = callOctave("pyulear",x,poles,freq,Fs,range,method,plot_type) + else + error("Wrong number of output argments.") + end + + end +endfunction + diff --git a/macros/qp_kaiser.sci b/macros/qp_kaiser.sci new file mode 100644 index 0000000..820fce8 --- /dev/null +++ b/macros/qp_kaiser.sci @@ -0,0 +1,37 @@ +function res = qp_kaiser (nb, at, linear) +//Computes a finite impulse response (FIR) filter for use with a quasi-perfect reconstruction polyphase-network filter bank. +//Calling Sequence +//qp_kaiser (nb, at, linear) +//qp_kaiser (nb, at) +//Parameters +//nb: Number of bands +//at: Attenuation +//linear: When not zero, minimum-phase calculation is omitted. +//Description +//This is an Octave function. +//This version utilizes a Kaiser window to shape the frequency response of the designed filter. Tha number nb of bands and the desired attenuation at in the stop-band are given as parameters. +// +//The Kaiser window is multiplied by the ideal impulse response h(n)=a.sinc(a.n) and converted to its minimum-phase version by means of a Hilbert transform. +//Examples +// qp_kaiser (5, 5, 1) +//ans = +// +// 0.11591 0.25606 0.25606 0.25606 0.11591 + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 2 | rhs > 3) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 2 then + res = callOctave("qp_kaiser", nb, at) + + case 3 then + res = callOctave("qp_kaiser", nb, at, linear) + + end +endfunction diff --git a/macros/sftrans.sci b/macros/sftrans.sci new file mode 100644 index 0000000..916d44d --- /dev/null +++ b/macros/sftrans.sci @@ -0,0 +1,74 @@ +function [Sz, Sp, Sg] = sftrans (Sz, Sp, Sg, W, stop) +//Transform band edges of a generic lowpass filter (cutoff at W=1) represented in splane zero-pole-gain form. +//Calling Sequence +//[Sz, Sp, Sg] = sftrans (Sz, Sp, Sg, W, stop) +//[Sz, Sp] = sftrans (Sz, Sp, Sg, W, stop) +//[Sz] = sftrans (Sz, Sp, Sg, W, stop) +//Parameters +//Sz: Zeros. +//Sp: Poles. +//Sg: Gain. +//W: Edge of target filter. +//stop: True for high pass and band stop filters or false for low pass and band pass filters. +//Description +//This is an Octave function. +//Theory: Given a low pass filter represented by poles and zeros in the splane, you can convert it to a low pass, high pass, band pass or band stop by transforming each of the poles and zeros +//individually. The following table summarizes the transformation: +// +// Transform Zero at x Pole at x +// ---------------- ------------------------- ------------------------ +// Low Pass zero: Fc x/C pole: Fc x/C +// S -> C S/Fc gain: C/Fc gain: Fc/C +// ---------------- ------------------------- ------------------------ +// High Pass zero: Fc C/x pole: Fc C/x +// S -> C Fc/S pole: 0 zero: 0 +// gain: -x gain: -1/x +// ---------------- ------------------------- ------------------------ +// Band Pass zero: b +- sqrt(b^2-FhFl) pole: b +- sqrt(b^2-FhFl) +// S^2+FhFl pole: 0 zero: 0 +// S -> C -------- gain: C/(Fh-Fl) gain: (Fh-Fl)/C +// S(Fh-Fl) b=x/C (Fh-Fl)/2 b=x/C (Fh-Fl)/2 +// ---------------- ------------------------- ------------------------ +// Band Stop zero: b +- sqrt(b^2-FhFl) pole: b +- sqrt(b^2-FhFl) +// S(Fh-Fl) pole: +-sqrt(-FhFl) zero: +-sqrt(-FhFl) +// S -> C -------- gain: -x gain: -1/x +// S^2+FhFl b=C/x (Fh-Fl)/2 b=C/x (Fh-Fl)/2 +// ---------------- ------------------------- ------------------------ +// Bilinear zero: (2+xT)/(2-xT) pole: (2+xT)/(2-xT) +// 2 z-1 pole: -1 zero: -1 +// S -> - --- gain: (2-xT)/T gain: (2-xT)/T +// T z+1 +// ---------------- ------------------------- ------------------------ +// +//where C is the cutoff frequency of the initial lowpass filter, Fc is the edge of the target low/high pass filter and [Fl,Fh] are the edges of the target band pass/stop filter. With abundant tedious +//algebra, you can derive the above formulae yourself by substituting the transform for S into H(S)=S-x for a zero at x or H(S)=1/(S-x) for a pole at x, and converting the result into the form: +// +// H(S)=g prod(S-Xi)/prod(S-Xj) +//Examples +//[Sz, Sp, Sg] = sftrans (5, 10, 15, 20, 30) +//Sz = 4 +//Sp = 2 +//Sg = 7.5000 + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 5 | rhs > 5) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 5 then + if(lhs==1) + Sz = callOctave("sftrans",Sz, Sp, Sg, W, stop) + elseif(lhs==2) + [Sz, Sp] = callOctave("sftrans",Sz, Sp, Sg, W, stop) + elseif(lhs==3) + [Sz, Sp, Sg] = callOctave("sftrans",Sz, Sp, Sg, W, stop) + else + error("Wrong number of output argments.") + end + + end +endfunction diff --git a/macros/tfestimate.sci b/macros/tfestimate.sci new file mode 100644 index 0000000..e52faf5 --- /dev/null +++ b/macros/tfestimate.sci @@ -0,0 +1,93 @@ +function [Pxx, freq] = tfestimate(x, y, window, overlap, Nfft, Fs, range) + +//Estimate transfer function of system with input x and output y. Use the Welch (1967) periodogram/FFT method. +//Calling Sequence +//tfestimate (x, y, window, overlap, Nfft, Fs, range) +//[Pxx, freq] = tfestimate (…) +//Parameters +//x: Input. +//y: Output. +//window: [real vector] of window-function values between 0 and 1; the data segment has the same length as the window. Default window shape is Hamming. [integer scalar] length of each data segment. The default value is window=sqrt(length(x)) rounded up to the nearest integer power of 2; see 'sloppy' argument. +// overlap: [real scalar] segment overlap expressed as a multiple of window or segment length. 0 <= overlap < 1, The default is overlap=0.5 . +// Nfft: [integer scalar] Length of FFT. The default is the length of the "window" vector or has the same value as the scalar "window" argument. If Nfft is larger than the segment length, "seg_len", the data segment is padded with "Nfft-seg_len" zeros. The default is no padding. Nfft values smaller than the length of the data segment (or window) are ignored silently. +// Fs: [real scalar] sampling frequency (Hertz); default=1.0 +// range: 'half', 'onesided' : frequency range of the spectrum is zero up to but not including Fs/2. Power from negative frequencies is added to the positive side of the spectrum, but not at zero or Nyquist (Fs/2) frequencies. This keeps power equal in time and spectral domains. See reference [2]. 'whole', 'twosided' : frequency range of the spectrum is -Fs/2 to Fs/2, with negative frequencies stored in "wrap around" order after the positive frequencies; e.g. frequencies for a 10-point 'twosided' spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1 'shift', 'centerdc' : same as 'whole' but with the first half of the spectrum swapped with second half to put the zero-frequency value in the middle. (See "help fftshift". If data (x and y) are real, the default range is 'half', otherwise default range is 'whole'. +//Description +//This function is being called from Octave. +//Estimate transfer function of system with input x and output y. Use the Welch (1967) periodogram/FFT method. +//Examples +//[Pxx, freq]=tfestimate ([1 2 3], [4 5 6]) +//Pxx = +// +// 1.7500 + 0.0000i +// 1.5947 + 0.3826i +// 1.2824 + 0.0000i +// +//freq = +// +// 0.00000 +// 0.25000 +// 0.50000 + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 2 | rhs > 7) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 2 then + if(lhs==1) + Pxx = callOctave("tfestimate",x,y) + elseif(lhs==2) + [Pxx, freq] = callOctave("tfestimate",x,y) + else + error("Wrong number of output argments.") + end + + case 3 then + if(lhs==1) + Pxx = callOctave("tfestimate",x, y, window) + elseif(lhs==2) + [Pxx, freq] = callOctave("tfestimate",x, y, window) + else + error("Wrong number of output argments.") + end + case 4 then + if(lhs==1) + Pxx = callOctave("tfestimate",x, y, window, overlap) + elseif(lhs==2) + [Pxx, freq] = callOctave("tfestimate",x, y, window, overlap) + else + error("Wrong number of output argments.") + end + case 5 then + if(lhs==1) + Pxx = callOctave("tfestimate",x, y, window, overlap, Nfft) + elseif(lhs==2) + [Pxx, freq] = callOctave("tfestimate",x, y, window, overlap, Nfft) + else + error("Wrong number of output argments.") + end + case 6 then + if(lhs==1) + Pxx = callOctave("tfestimate",x, y, window, overlap, Nfft, Fs) + elseif(lhs==2) + [Pxx, freq] = callOctave("tfestimate",x, y, window, overlap, Nfft, Fs) + else + error("Wrong number of output argments.") + end + case 7 then + if(lhs==1) + Pxx = callOctave("tfestimate",x, y, window, overlap, Nfft, Fs, range) + elseif(lhs==2) + [Pxx, freq] = callOctave("tfestimate",x, y, window, overlap, Nfft, Fs, range) + else + error("Wrong number of output argments.") + end + end +endfunction + + diff --git a/macros/xcorr1.sci b/macros/xcorr1.sci new file mode 100644 index 0000000..6ae884d --- /dev/null +++ b/macros/xcorr1.sci @@ -0,0 +1,86 @@ +function [R, lag] = xcorr1 (X, Y, maxlag, scale) +//Estimates the cross-correlation. +//Calling Sequence +//[R, lag] = xcorr1 (X, Y, maxlag, scale) +//[R, lag] = xcorr1 (X, Y, maxlag) +//[R, lag] = xcorr1 (X, Y) +//Parameters +//X: [non-empty; real or complex; vector or matrix] data. +//Y: [real or complex vector] data. +// If X is a matrix (not a vector), Y must be omitted. Y may be omitted if X is a vector; in this case xcorr estimates the autocorrelation of X. +//maxlag: [integer scalar] maximum correlation lag If omitted, the default value is N-1, where N is the greater of the lengths of X and Y or, if X is a matrix, the number of rows in X. +//scale: [character string] specifies the type of scaling applied to the correlation vector (or matrix). is one of: +// ‘none’ - return the unscaled correlation, R, +//‘biased’ - return the biased average, R/N, +//‘unbiased’ - return the unbiased average, R(k)/(N-|k|), +//‘coeff’ - return the correlation coefficient, R/(rms(x).rms(y)), where "k" is the lag, and "N" is the length of X. If omitted, the default value is "none". If Y is supplied but does not have the same +// length as X, scale must be "none". +//Description +//This is an Octave function. +//Estimate the cross correlation R_xy(k) of vector arguments X and Y or, if Y is omitted, estimate autocorrelation R_xx(k) of vector X, for a range of lags k specified by argument "maxlag". If X is a +//matrix, each column of X is correlated with itself and every other column. +// +//The cross-correlation estimate between vectors "x" and "y" (of length N) for lag "k" is given by +// +// N +// R_xy(k) = sum x_{i+k} conj(y_i), +// i=1 +// +//where data not provided (for example x(-1), y(N+1)) is zero. Note the definition of cross-correlation given above. To compute a cross-correlation consistent with the field of statistics, see xcov. +//Examples +//[R, lag] = xcorr1 ( [5 5], [2 2], 2, 'biased' ) +// +//R = +// +// 0 5 10 5 0 +// +//lag = +// +// -2 -1 0 1 2 + +funcprot(0); +lhs = argn(1) +rhs = argn(2) +if (rhs < 1 | rhs > 4) +error("Wrong number of input arguments.") +end + +select(rhs) + + case 2 then + if(lhs==1) + R = callOctave("xcorr", X) + elseif(lhs==2) + [R, lag] = callOctave("xcorr", X) + else + error("Wrong number of output argments.") + end + + case 2 then + if(lhs==1) + R = callOctave("xcorr", X, Y) + elseif(lhs==2) + [R, lag] = callOctave("xcorr", X, Y) + else + error("Wrong number of output argments.") + end + + case 3 then + if(lhs==1) + R = callOctave("xcorr", X, Y, maxlag) + elseif(lhs==2) + [R, lag] = callOctave("xcorr", X, Y, maxlag) + else + error("Wrong number of output argments.") + end + case 4 then + if(lhs==1) + R = callOctave("xcorr", X, Y, maxlag, scale) + elseif(lhs==2) + [R, lag] = callOctave("xcorr", X, Y, maxlag, scale) + else + error("Wrong number of output argments.") + end + + end +endfunction |