summaryrefslogtreecommitdiff
path: root/macros
diff options
context:
space:
mode:
Diffstat (limited to 'macros')
-rw-r--r--macros/ar_psd.sci301
-rw-r--r--macros/arma_rnd.sci75
-rw-r--r--macros/autoreg_matrix.sci68
-rw-r--r--macros/cceps.sci92
-rw-r--r--macros/clustersegment.sci90
-rw-r--r--macros/detrend1.sci93
-rw-r--r--macros/dst1.sci76
-rw-r--r--macros/durbinlevinson.sci111
-rw-r--r--macros/fftshift1.sci111
-rw-r--r--macros/filter2.sci93
-rw-r--r--macros/fwhm.sci165
-rw-r--r--macros/fwhmjlt.sci171
-rw-r--r--macros/hurst.sci64
-rw-r--r--macros/ifftshift1.sci105
-rw-r--r--macros/ifht.sci90
-rw-r--r--macros/morlet.sci93
-rw-r--r--macros/pburg.sci215
-rw-r--r--macros/pei_tseng_notch.sci111
-rw-r--r--macros/polystab.sci53
-rw-r--r--macros/schtrig.sci153
-rw-r--r--macros/sinetone.sci94
-rw-r--r--macros/sinewave.sci84
-rw-r--r--macros/spencer.sci50
-rw-r--r--macros/stft.sci206
-rw-r--r--macros/synthesis.sci90
-rw-r--r--macros/tf2sos.sci83
-rw-r--r--macros/tfestimate.sci132
-rw-r--r--macros/triang.sci14
-rw-r--r--macros/tripuls.sci101
-rw-r--r--macros/zp2sos.sci287
30 files changed, 2489 insertions, 982 deletions
diff --git a/macros/ar_psd.sci b/macros/ar_psd.sci
index 540dfec..1953043 100644
--- a/macros/ar_psd.sci
+++ b/macros/ar_psd.sci
@@ -1,40 +1,263 @@
-function [P, F]= ar_psd(A, varargin)
-//Calculate the power spectrum of the autoregressive model
-//Calling Sequence
-// [PSD,F_OUT]=ar_psd (A, V)
-// [PSD,F_OUT]=ar_psd (A, V, FREQ)
-// [PSD,F_OUT]=ar_psd (A, V, FREQ, FS)
-// [PSD,F_OUT]=ar_psd (..., RANGE)
-// [PSD,F_OUT]=ar_psd (..., METHOD)
-// [PSD,F_OUT]=ar_psd (..., PLOTTYPE)
-//Parameters
-//A:List of M=(order+1) autoregressive model coefficients. The first element of "ar_coeffs" is the zero-lag coefficient, which always has a value of 1.
-//V:Square of the moving-average coefficient of the AR model.
-//FREQ:Frequencies at which power spectral density is calculated, or a scalar indicating the number of uniformly distributed frequency values at which spectral density is calculated. (default = 256)
-//FS: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 spectrum. 'poly': calculate power spectrum 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.
-//PSD: estimate of power-spectral density
-//F_OUT: frequency values
-//Description
-//If the FREQ argument is a vector (of frequencies) the spectrum is calculated using the polynomial method and the METHOD argument is ignored. For scalar FREQ, an integer power of 2, or METHOD = "FFT", causes the spectrum to be calculated by FFT. Otherwise, the spectrum is calculated as a polynomial. It may be computationally more efficient to use the FFT method if length of the model is not much smaller than the number of frequency values. The spectrum is scaled so that spectral energy (area under spectrum) is the same as the time-domain energy (mean square of the signal).
-//Examples
-//[a,b]= ar_psd([1,2,3],2)
-
- funcprot(0);
- rhs= argn(2);
- if(rhs <2 | rhs>5)
- error("Wrong number of input arguments");
- end
- select(rhs)
- case 2 then
- [P,F]= callOctave("ar_psd", A, varargin(1));
- case 3 then
- [P,F]= callOctave("ar_psd", A, varargin(1), varargin(2));
- case 4 then
- [P,F]= callOctave("ar_psd", A, varargin(1), varargin(2), varargin(3));
- case 5 then
- [P,F]= callOctave("ar_psd", A, varargin(1), varargin(2), varargin(3), varargin(4));
- end
+function varargout = ar_psd(a, v, varargin)
+//Calculate the power spectrum of the autoregressive model.
+
+//Calling Sequence:
+// [psd, f_out] = ar_psd(a, v)
+// [psd, f_out] = ar_psd (a, v, freq)
+// [psd, f_out] = ar_psd (a, v, freq, fs)
+// [psd, f_out] = ar_psd (..., range)
+// [psd, f_out] = ar_psd (..., method)
+// [psd, f_out] = ar_psd (..., plottype)
+
+//Parameters:
+//Every parameter except for the first two is optional.
+//
+//a- List of m=(order + 1) autoregressive model coefficients. The first element of "ar_coeffs" is the zero-lag coefficient, which always has a value of 1.
+//v- Square of the moving-average coefficient of the AR model.
+//freq: Frequencies at which power spectral density is calculated, or a scalar indicating the number of uniformly distributed frequency values at which spectral density is calculated. (default = 256)
+//fs- 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. 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 spectrum.
+//'poly'- calculate power spectrum 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.
+//psd: estimate of power-spectral density.
+//f_out: frequency values.
+
+//Description:
+//If the 'freq' argument is a vector (of frequencies) the spectrum is calculated using the polynomial method and the METHOD argument is ignored. For scalar 'freq', an integer power of 2, or method = "fft", causes the spectrum to be calculated by fft. Otherwise, the spectrum is calculated as a polynomial. It may be computationally more efficient to use the fft methodif length of the model is not much smaller than the number of frequency values. The spectrum is scaled so that spectral energy (area under spectrum) is the same as the time-domain energy (mean square of the signal).
+
+//Examples:
+//[a,b]= ar_psd([1,2,3], 2)
+
+ funcprot(0);
+ // Check fixed arguments
+ if nargin < 2 then
+ error("ar_psd: needs at least 2 args. Use help ar_psd.");
+ elseif ~isvector(a) | length(a) < 2 then
+ error("ar_psd: arg 1 (a) must be vector, length >= 2.");
+ elseif ~isscalar(v) then
+ error("ar_psd: arg 2 (v) must be real scalar >0.");
+ else
+ real_model = isreal(a);
+ // Default values for optional arguments
+ freq = 256;
+ user_freqs = 0; // Boolean: true for user-specified frequencies
+ Fs = 1.0;
+ // FFT padding factor (is also frequency range divisor): 1=whole, 2=half.
+ pad_fact = 1 + real_model;
+ do_shift = 0;
+ force_FFT = 0;
+ force_poly = 0;
+ plot_type = 1;
+ // Decode and check optional arguments
+ end_numeric_args = 0;
+ for iarg = 1:length(varargin)
+ arg = varargin(iarg);
+ end_numeric_args = end_numeric_args | (type(arg) == 10);
+ // Skip empty arguments
+ if isempty(arg) then
+ // Do nothing
+ elseif (type(arg) ~= 10) then
+ if end_numeric_args then
+ error("ar_psd: control arg must be string.");
+ // First optional numeric arg is "freq"
+ elseif iarg == 1 then
+ user_freqs = isvector(arg) & length(arg) > 1;
+ if ~isscalar(arg) & ~user_freqs then
+ error("ar_psd: arg 3 (freq) must be vector or scalar.");
+ elseif ~user_freqs & (~isreal(arg) | fix(arg) ~= arg | arg <= 2 | arg >= 1048576) then
+ error("ar_psd: arg 3 (freq) must be integer >=2, <=1048576");
+ elseif user_freqs & ~isreal(arg) then
+ error("ar_psd: arg 3 (freq) vector must be real.");
+ end
+ freq = arg(:); // -> column vector
+ // Second optional numeric arg is "Fs" - sampling frequency
+ elseif iarg == 2 then
+ if ~isscalar(arg) | ~isreal(arg) | arg <= 0 then
+ error("ar_psd: arg 4 (Fs) must be real positive scalar.");
+ end
+ Fs = arg;
+ else
+ error("ar_psd: control arg must be string.");
+ end
+ // Decode control-string arguments
+ elseif ~strcmp(arg, "plot") | ~strcmp(arg, "squared") then
+ plot_type = 1;
+ elseif ~strcmp(arg, "semilogx") then
+ plot_type = 2;
+ elseif ~strcmp(arg, "semilogy") then
+ plot_type = 3;
+ elseif ~strcmp(arg, "loglog") then
+ plot_type = 4;
+ elseif ~strcmp(arg, "dB") then
+ plot_type = 5;
+ elseif ~strcmp(arg, "fft") then
+ force_FFT = 1;
+ force_poly = 0;
+ elseif ~strcmp(arg, "poly") then
+ force_FFT = 0;
+ force_poly = 1;
+ elseif ~strcmp(arg, "half") | ~strcmp(arg, "onesided") then
+ pad_fact = 2; // FFT zero-padding factor (pad FFT to double length)
+ do_shift = 0;
+ elseif ~strcmp(arg, "whole") | ~strcmp(arg, "twosided") then
+ pad_fact = 1; // FFT zero-padding factor (do not pad)
+ do_shift = 0;
+ elseif ~strcmp(arg, "shift") | ~strcmp(arg, "centerdc") then
+ pad_fact = 1;
+ do_shift = 1;
+ else
+ error("ar_psd: string arg: illegal value: %s", arg);
+ end
+ end
+ // End of decoding and checking args
+ if user_freqs then
+ // User provides (column) vector of frequencies
+ if or(abs(freq) > Fs/2) then
+ error("ar_psd: arg 3 (freq) cannot exceed half sampling frequency.");
+ elseif pad_fact == 2 & or(freq < 0) then
+ error("ar_psd: arg 3 (freq) must be positive in onesided spectrum");
+ end
+ freq_len = length(freq);
+ fft_len = freq_len;
+ use_FFT = 0;
+ do_shift = 0;
+ else
+ // Internally generated frequencies
+ freq_len = freq;
+ freq = (Fs / pad_fact / freq_len) * (0:freq_len - 1)';
+ // Decide which method to use (poly or FFT)
+ is_power_of_2 = modulo(log(freq_len), log(2)) < 10 * %eps;
+ use_FFT = (~force_poly & is_power_of_2) | force_FFT;
+ fft_len = freq_len * pad_fact;
+ end
+ // Calculate denominator of Equation 2.28, Kay and Marple, ref [1] Jr.:
+ len_coeffs = length(a);
+ if use_FFT then
+ // FFT method
+ x = [a(:); zeros(fft_len - len_coeffs, 1)];
+ fft_out = fft(x);
+ else
+ // Polynomial method
+ // Complex data on "half" frequency range needs -ve frequency values
+ if pad_fact == 2 & ~real_model then
+ freq = [freq; -freq(freq_len:-1:1)];
+ fft_len = 2 * freq_len;
+ end
+ fft_out = horner(a($:-1:1), exp((-%i * 2 * %pi / Fs) * freq));
+ end
+ // The power spectrum (PSD) is the scaled squared reciprocal of amplitude
+ // of the FFT/polynomial. This is NOT the reciprocal of the periodogram.
+ // The PSD is a continuous function of frequency. For uniformly
+ // distributed frequency values, the FFT algorithm might be the most
+ // efficient way of calculating it.
+ psd = (v / Fs) ./ (fft_out .* conj(fft_out));
+ // Range='half' or 'onesided',
+ // Add PSD at -ve frequencies to PSD at +ve frequencies
+ // N.B. unlike periodogram, PSD at zero frequency _is_ doubled.
+ if pad_fact == 2 then
+ freq = freq(1:freq_len);
+ if real_model then
+ // Real data, double the psd
+ psd = 2 * psd(1:freq_len);
+ elseif use_FFT then
+ // Complex data, FFT method, internally-generated frequencies
+ psd = psd(1:freq_len) + [psd(1); psd(fft_len:-1:freq_len + 2)];
+ else
+ // Complex data, polynomial method
+ // User-defined and internally-generated frequencies
+ psd = psd(1:freq_len) + psd(fft_len:-1:freq_len + 1);
+ end
+ // Range='shift'
+ // Disabled for user-supplied frequencies
+ // Shift zero-frequency to the middle (pad_fact == 1)
+ elseif do_shift then
+ len2 = fix((fft_len + 1) / 2);
+ psd = [psd(len2 + 1:fft_len); psd(1:len2)];
+ freq = [freq(len2 + 1:fft_len) - Fs; freq(1:len2)];
+ end
+ // Plot the spectrum if there are no return variables.
+ if nargout() >= 2 then
+ varargout(1) = psd;
+ varargout(2) = freq;
+ elseif nargout() == 1 then
+ varargout(1) = psd;
+ else
+ if plot_type == 1 then
+ plot(freq, psd);
+ elseif plot_type == 2 then
+ semilogx(freq, psd);
+ elseif plot_type == 3 then
+ semilogy(freq, psd);
+ elseif plot_type == 4 then
+ loglog(freq, psd);
+ elseif plot_type == 5 then
+ plot(freq, 10 * log10(psd));
+ end
+ end
+ end
endfunction
+
+//tests:
+
+//a = [1, -0.5];
+////v = 1;
+//[psd, freq] = ar_psd(a, v);
+//plot(freq, psd);
+//title('Power Spectral Density of the AR Model');
+//xlabel('Frequency');
+//ylabel('Power/Frequency');
+
+//a = [1, -1.5, 0.7];
+//v = 2;
+//Fs = 2.0;
+//[psd, freq] = ar_psd(a, v, 512, Fs);
+//plot(freq, psd);
+//title('Power Spectral Density with Different Sampling Frequency');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency');
+
+//a = [1, -0.9, 0.4];
+//v = 0.8;
+//Fs = 1.0;
+//ar_psd(a, v, 512, Fs, 'semilogx');
+//title('Power Spectral Density (Semilogx)');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency');
+//
+//figure();
+//ar_psd(a, v, 512, Fs, 'loglog');
+//title('Power Spectral Density (Loglog)');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency');
+
+//a = [1, -0.7, 0.2];
+//v = 1.5;
+//Fs = 1.0;
+//[psd_fft, freq] = ar_psd(a, v, 512, Fs, 'fft');
+//[psd_poly, freq] = ar_psd(a, v, 512, Fs, 'poly');
+//plot(freq, psd_fft, 'r', freq, psd_poly, 'b');
+//title('Power Spectral Density (FFT vs Polynomial)');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency');
+//legend('FFT Method', 'Polynomial Method');
+
+//a = [1, -1.2, 0.5];
+//v = 1;
+//[psd_half, freq_half] = ar_psd(a, v, 512, 1, 'half');
+//[psd_whole, freq_whole] = ar_psd(a, v, 512, 1, 'whole');
+//subplot(2, 1, 1);
+//plot(freq_half, psd_half);
+//title('Power Spectral Density (Half Spectrum)');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency');
+//
+//subplot(2, 1, 2);
+//plot(freq_whole, psd_whole);
+//title('Power Spectral Density (Whole Spectrum)');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency');
diff --git a/macros/arma_rnd.sci b/macros/arma_rnd.sci
index 5e9ca12..fb428cc 100644
--- a/macros/arma_rnd.sci
+++ b/macros/arma_rnd.sci
@@ -1,27 +1,20 @@
+function x = arma_rnd (a, b, v, t, n)
//Return a simulation of the ARMA model.
-
-//Calling Sequence
+//Calling Sequence:
//arma_rnd (a, b, v, t, n)
//arma_rnd (a, b, v, t)
-
-//Parameters
-//a: vector
-//b: vector
+//Parameters:
+//a: Real vector
+//b: Real 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)
+//Description:
+//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
+//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;
@@ -30,36 +23,13 @@
//arma_rnd (a, b, v, t, n)
//Output :
// ans =
-//
// 61400.907
// 158177.11
// 407440.29
// 1049604.
// 2703841.3
-
-//function res = arma_rnd (a, b, v, t, n)
-//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
-
-function x = arma_rnd (a, b, v, t, n)
-
- funcprot();
+ funcprot(0);
[nargout,nargin] = argn() ;
if (nargin == 4)
@@ -69,7 +39,7 @@ function x = arma_rnd (a, b, v, t, n)
error ("arma_rnd: N must be a scalar");
end
else
- error("arma_rnd: invalid input");
+ error("Wrong number of input arguments.");
end
if ((min (size (a)) > 1) | (min (size (b)) > 1))
@@ -85,16 +55,31 @@ function x = arma_rnd (a, b, v, t, n)
a = matrix (a, ar, 1);
b = matrix (b, br, 1);
-
- // Apply our notational convention.
a = [1; -a];
b = [1; b];
n = min (n, ar + br);
-
e = sqrt (v) * rand(t + n, 1);
x = filter (b, a, e);
x = x(n + 1 : t + n);
endfunction
+
+//input validation:
+//a = [1, 2, 3, 4];
+//b = [5, 6, 7];
+//assert_checkerror("arma_rnd()", "Wrong number of input arguments.");
+//assert_checkerror("arma_rnd(1, 2, 3, 4, 5, 6)", "Wrong number of input arguments.");
+//assert_checkerror("arma_rnd(a, b, 5, 2, a);", "arma_rnd: N must be a scalar");
+//assert_checkerror("arma_rnd(a, b, 5, a);", "arma_rnd: T must be a scalar");
+//assert_checkerror("arma_rnd([1 2; 3 4], [5 6; 7 8], 5, 1);", "arma_rnd: A and B must not be matrices");
+
+//tests:
+//NOTE: The output of arma_rnd is supposed to be random, so we cannot expect the same output for equivalent input.
+//a = [1, 2, 3, 4];
+//b = [5, 6, 7, 8];
+//assert_checkequal(size(arma_rnd(a, b, 5, 1)), [1 1]);
+//assert_checkequal(size(arma_rnd(a, b, 5, 2, 100)), [2 1]);
+//assert_checkequal(size(arma_rnd(a', b', 1, 10, 50)), [10 1]);
+//assert_checkequal(size(arma_rnd(a, b', 1, 4, 5)), [4 1]);
diff --git a/macros/autoreg_matrix.sci b/macros/autoreg_matrix.sci
index 2b70ced..95eb356 100644
--- a/macros/autoreg_matrix.sci
+++ b/macros/autoreg_matrix.sci
@@ -1,59 +1,53 @@
+function x = autoreg_matrix (y, k)
// Given a time series (vector) Y, return a matrix with ones in the first column and the first K lagged values of Y in the other columns.
-
-//Calling Sequence
+//Calling Sequence:
//autoreg_matrix(Y, K)
-
-//Parameters
-//Y: Vector
-//K: Scalar or Vector
-
-//Description
+//Parameters:
+//Y: vector
+//K: scalar
+//Description:
// Given a time series (vector) Y, return a matrix with ones in the first column and the first K lagged values of Y in the other columns.
-//
//In other words, for T > K, '[1, Y(T-1), ..., Y(T-K)]' is the t-th row of the result.
-//
//The resulting matrix may be used as a regressor matrix in autoregressions.
-
-//Examples
-//autoreg_matrix([1,2,3],2)
+//Examples:
+//autoreg_matrix([1,2,3], 2)
//ans =
// 1. 0. 0.
// 1. 1. 0.
// 1. 2. 1.
-
-
-//function y = autoreg_matrix(Y, varargin)
-//funcprot(0);
-//rhs = argn(2)
-//if(rhs<2 | rhs>2)
-//error("Wrong number of input arguments.");
-//end
-//
-// select(rhs)
-// case 2 then
-// y = callOctave("autoreg_matrix", Y, varargin(1));
-// end
-//endfunction
-
-function X = autoreg_matrix (y, k)
-
funcprot(0);
- [nargout, nargin] = argn() ;
- if (nargin ~= 2)
- error('autoreg_matrix: invalid input') ;
+ if (argn(2) ~= 2)
+ error("autoreg_matrix: wrong number of input arguments") ;
end
if (~ (isvector (y)))
- error ("autoreg_matrix: Y must be a vector");
+ error ("autoreg_matrix: y must be a vector");
end
T = length (y);
y = matrix(y, T, 1);
- X = ones (T, k+1);
- for j = 1 : k;
- X(:, j+1) = [(zeros (j, 1)); y(1:T-j)];
+ x = ones (T, k+1);
+ for j = 1 : k
+ x(:, j+1) = [(zeros(j, 1)); y(1:T-j)];
end
endfunction
+
+//input validation:
+//assert_checkerror("autoreg_matrix(1)", "autoreg_matrix: wrong number of input arguments");
+//assert_checkerror("autoreg_matrix(1, 2, 3)", "Wrong number of input arguments.");
+//assert_checkerror("autoreg_matrix(1, 2)", "autoreg_matrix: y must be a vector");
+//assert_checkerror("autoreg_matrix([1, 2; 3, 4], 2)", "autoreg_matrix: y must be a vector");
+
+//tests:
+//assert_checkequal(autoreg_matrix([1, 2], -1), []);
+//assert_checkequal(autoreg_matrix([1, 2, 3], 2), [1, 0, 0; 1, 1, 0; 1, 2, 1]);
+//assert_checkequal(autoreg_matrix([1, 2, 3], 2), autoreg_matrix([1; 2; 3], 2));
+//assert_checkequal(autoreg_matrix([1, 2, 3, 4, 5], 0), [1; 1; 1; 1; 1]);
+//assert_checkequal(autoreg_matrix([-1; -3; -5; -7; -9], 5), [1 0 0 0 0 0;1 -1 0 0 0 0;1 -3 -1 0 0 0;1 -5 -3 -1 0 0;1 -7 -5 -3 -1 0])
+//assert_checkequal(autoreg_matrix([1+2*%i, 5+4*%i, -4*%i, -1-6*%i], 1), [1, 0; 1, 1 + 2*%i; 1, 5 + 4*%i; 1, -4*%i]);
+//assert_checkequal(autoreg_matrix([1+2*%i, 5+4*%i, -4*%i, -1-6*%i], 3), autoreg_matrix([1+2*%i; 5+4*%i; -4*%i; -1-6*%i], 3));
+//assert_checkequal(autoreg_matrix([-%i; -3-%i; 5+6*%i; 7+9*%i;], 3), [1 0 0 0;1 -%i 0 0;1 -3-%i -%i 0;1 5+6*%i -3-%i -%i]);
+//assert_checkequal(autoreg_matrix([-%i; -3-%i; 5+6*%i; 7+9*%i;], 0), [1; 1; 1; 1]);
diff --git a/macros/cceps.sci b/macros/cceps.sci
index f026256..b87d4bf 100644
--- a/macros/cceps.sci
+++ b/macros/cceps.sci
@@ -1,63 +1,75 @@
-function y = cceps (x,correct)
-//Return the complex cepstrum of the vector x
-//Calling Sequence
+function y = cceps(x, correct)
+//Return the complex cepstrum of the vector x.
+//Calling Sequence:
//cceps (x)
//cceps(x, correct)
-//Parameters
-//x: vector.
-//correct: if 1, a correction method is applied.
-//Description
-//This function return the complex cepstrum of the vector x. If the optional argument correct has the value 1, a correction method is applied. The default is not to do this.
-//Examples
-//cceps([1,2,3],1)
-//ans = 1.9256506
- // 0.9634573
-// - 1.0973484
-
-if(argn(2)<1 | argn(2)>2)
-error("Wrong number of input arguments.")
-end
- if argn(2)==1 then
- correct=0;
-end
-
- [r, c]=size(x)
- if (c ~= 1)
+//Parameters:
+//x: vector
+//correct: If 1, a correction method is applied.
+//Description:
+//This function return the complex cepstrum of the vector x.
+//If the optional argument correct has the value 1, a correction method is applied. The default is not to do this.
+//Examples:
+//cceps([1,2,3], 1)
+//ans = 1.9256506
+// 0.9634573
+// -1.0973484
+
+ if(argn(2) < 1 | argn(2) > 2)
+ error("Wrong number of input arguments.");
+ end
+
+ if (argn(2) == 1)
+ correct = 0;
+ end
+
+ [r, c] = size(x);
+ if (c ~= 1)
if (r == 1)
- x = x;
+ x = x';
r = c;
else
error ("x must be a vector");
end
end
-
- F = fft1(x);
- if (min (abs (F)) == 0)
- error ('bad signal x, some Fourier coefficients are zero.');
+
+ F = fft(x);
+ if (min(abs(F)) == 0)
+ error('bad signal x, some Fourier coefficients are zero.');
end
- h = fix (r / 2);
+
+ h = fix (r / 2);
cor = 0;
- if (2 * h == r)
- cor = (c & (real (F (h + 1)) < 0));
+ if (2*h == r)
+ cor = (correct & (real(F(h + 1)) < 0));
if (cor)
- F = fft1 (x(1:r-1))
- if (min (abs (F)) == 0)
- error ('bad signal x, some Fourier coefficients are zero.');
+ F = fft(x(1:r-1));
+ if (min(abs(F)) == 0)
+ error('bad signal x, some Fourier coefficients are zero.');
end
end
end
- y = fftshift1 (ifft1 ((log (F))'));
+ y = fftshift(ifft((log (F))));
- //## make result real
- if (c)
- y = real (y);
+ if (correct)
+ y = real(y);
if (cor)
- // ## make cepstrum of same length as input vector
y (r) = 0;
end
end
-
endfunction
+//input validation:
+//assert_checkerror("cceps()", "Wrong number of input arguments.");
+//assert_checkerror("cceps(1, 2, 3)", "Wrong number of input arguments.");
+//assert_checkerror("cceps([1, 2; 3, 4])", "x must be a vector");
+//assert_checkerror("cceps(0)", 'bad signal x, some Fourier coefficients are zero.');
+//assert_checkerror("cceps(zeros(10, 1))", 'bad signal x, some Fourier coefficients are zero.');
+//tests:
+//assert_checkalmostequal(cceps([1, 2, 3]), [1.9257; 0.9635; -1.0973], 0.5*10^-4);
+//assert_checkequal(cceps([1, 2, 3]), cceps([1, 2, 3], 1));
+//assert_checkequal(cceps([-1, -2, -3]), cceps([-1; -2; -3]));
+//assert_checkalmostequal(cceps([1+2*%i; -2-2*%i; -3+2*%i]), [0.4734+1.1174*%i; 1.2144+1.5358*%i; -0.1899+0.0247*%i], 5*10^-3);
+//assert_checkalmostequal(cceps([1+2*%i; -2-2*%i; -3+2*%i], 1), [0.4734; 1.2144; -0.1899], 5*10^-4);
diff --git a/macros/clustersegment.sci b/macros/clustersegment.sci
index 07957ee..811930a 100644
--- a/macros/clustersegment.sci
+++ b/macros/clustersegment.sci
@@ -1,27 +1,83 @@
-function c = clustersegment(s)
+function contRange = clustersegment(xhi)
//This function calculates boundary indexes of clusters of 1’s.
-//Calling Sequence
-//c = clustersegment(s)
+//Calling Sequence:
+//contRange = clustersegment(xhi)
//Parameters
-//s: scalar, vector or matrix of real numbers (clusters of 1s)
-//c: output variable, cell array of size 1 by N, where N is the number of rows in s
-//Description
-//This is an Octave function.
+//xhi: scalar, vector or matrix of real numbers (clusters of 1s)
+//contRange: output variable, cell array of size 1 by Np, where Np is the number of rows in 'xhi'
+//Description:
//This function calculates boundary indexes of clusters of 1’s.
//This function calculates the initial and end indices of the sequences of 1's present in the input argument.
-//The output variable c is a cell array of size 1 by N, where N is the number of rows in s and each element has two rows indicating the initial index and end index of the cluster of 1's respectively. The indexing starts from 1.
-//Examples
-//y = clustersegment ([0,1,0,0,1,1])
+//The output variable 'contRange' is a cell array of size 1 by Np, where Np is the number of rows in 'xhi' and each element has two rows indicating the initial index and end index of the cluster of 1's respectively. The indexing starts from 1.
+//Examples:
+//y = clustersegment([0,1,0,0,1,1])
//y =
// 2. 5.
-// 2. 6.
+// 2. 6.
-funcprot(0);
-rhs = argn(2)
-if(rhs~=1)
-error("Wrong number of input arguments.")
-end
+ funcprot(0);
+ if (argn(2) ~= 1)
+ error("Wrong number of input arguments.");
+ end
+ warning('off');
+
+ bool_discon = diff (xhi, 1, 2);
+ [Np Na] = size (xhi);
+ contRange = cell (1, Np);
-c = callOctave("clustersegment", s)
+ for i = 1:Np
+ idxUp = find (bool_discon(i,:) > 0) + 1;
+ idxDwn = find (bool_discon(i,:) < 0);
+ tLen = length (idxUp) + length (idxDwn);
+
+ if (xhi(i,1) == 1)
+ contRange{i}(1) = 1;
+ contRange{i}(2:2:tLen+1) = idxDwn;
+ contRange{i}(3:2:tLen+1) = idxUp;
+ else
+ contRange{i}(1:2:tLen) = idxUp;
+ contRange{i}(2:2:tLen) = idxDwn;
+ end
+
+ if (xhi(i, $) == 1)
+ contRange{i}($+1) = Na;
+ end
+
+ tLen = length (contRange{i});
+ if (tLen ~= 0)
+ contRange{i} = matrix(contRange{i}, 2, tLen / 2);
+ end
+
+ end
+
+ if (Np == 1)
+ contRange = cell2mat (contRange);
+ end
endfunction
+
+//tests:
+//assert_checkerror("clustersegment()", "Wrong number of input arguments.");
+//assert_checkerror("clustersegment(1, 2)", "Wrong number of input arguments.");
+//
+//assert_checkequal(clustersegment(1), [1; 1]);
+//assert_checkequal(clustersegment(-5), []);
+//assert_checkequal(clustersegment(3*%i), []);
+//assert_checkequal(clustersegment([0 0 1 1 1 0 0 1 0 0 0 1 1]), [3 8 12; 5 8 13]);
+//
+//ranges = clustersegment([-1; 1; 2; 1]);
+//assert_checkequal(ranges{1, 1}, []);
+//assert_checkequal(ranges{1, 2}, [1; 1]);
+//assert_checkequal(ranges{1, 3}, []);
+//assert_checkequal(ranges{1, 4}, [1; 1]);
+//
+//ranges = clustersegment([-1-2*%i; 1; 2*%i; 3+4*%i])
+//assert_checkequal(ranges{1, 1}, []);
+//assert_checkequal(ranges{1, 2}, [1; 1]);
+//assert_checkequal(ranges{1, 3}, []);
+//assert_checkequal(ranges{1, 4}, []);
+//
+//ranges = clustersegment([-1 1 1; 1 -1 1; -1 -1 1])
+//assert_checkequal(ranges{1, 1}, [2; 3]);
+//assert_checkequal(ranges{1, 2}, [1 3; 1 3]);
+//assert_checkequal(ranges{1, 3}, [3; 3]);
diff --git a/macros/detrend1.sci b/macros/detrend1.sci
index 52ac26f..29b0201 100644
--- a/macros/detrend1.sci
+++ b/macros/detrend1.sci
@@ -1,23 +1,72 @@
-function y = detrend1(x, varargin)
-//This function removes the best fit of a polynomial of order P from the data X
-//Calling Sequence
-//detrend1(X,P)
-//Parameters
-//X: Input vecor or matrix.
-//P: The order of polnomial
-//Description
-//If X is a vector, 'detrend1(X, P)' removes the best fit of apolynomial of order P from the data X.If X is a matrix, 'detrend1(X, P)' does the same for each column in X.
-//
-//The second argument P is optional. If it is not specified, a value of 1 is assumed. This corresponds to removing a linear trend.
-//The order of the polynomial can also be given as a string, in which case P must be either "constant" (corresponds to 'P=0') or "linear"(corresponds to 'P=1')
- rhs= argn(2);
- if(rhs<1 | rhs> 2)
- error("Wrong number of input arguments");
- end
- select(rhs)
- case 1 then
- y= callOctave("detrend", x);
- case 2 then
- y= callOctave("detrend", x , varargin(1));
- end
+function y = detrend1 (x, p)
+//Remove the best fit of a polynomial of order p from the data x.
+//Calling Sequence:
+//detrend1(x,p)
+//Parameters:
+//x: Input vecor or matrix
+//p: The order of polnomial
+//Description:
+//If X is a vector, 'detrend1(X, P)' removes the best fit of apolynomial of order P from the data X.
+//If X is a matrix, 'detrend1(X, P)' does the same for each column in X.
+//The second argument p is optional. If it is not specified, a value of 1 is assumed. This corresponds to removing a linear trend.
+//The order of the polynomial can also be given as a string, in which case p must be either "constant" (corresponds to 'P=0') or "linear" (corresponds to 'P=1')
+//Example:
+//detrend1([1, 6, 9])
+//ans = [ -0.3333, 0.6667, -0.3333]
+
+ funcprot(0);
+ rhs = argn(2);
+ if (rhs < 1 | rhs > 2)
+ error("detrend1: wrong number of input arguments");
+ end
+
+ if (rhs == 1)
+ p = 1;
+ end
+
+ if (~or(type(x)==[1 5 8]) | ndims (x) > 2)
+ error ("detrend1: X must be a numeric vector or matrix");
+ end
+
+ if (type(p) == 10 & ~strcmp(p, "constant", 'i'))
+ p = 0;
+ elseif (type(p) == 10 & ~strcmp(p, "linear", 'i'))
+ p = 1;
+ elseif (~isscalar(p) | p < 0 | p ~= fix (p))
+ error ("detrend1: P must be constant, linear, or a positive integer");
+ end
+
+ [m, n] = size (x);
+ if (m == 1)
+ x = x.';
+ end
+
+ r = size(x, 'r');
+ b = ((1 : r).' * ones (1, p + 1)) .^ (ones (r, 1) * (0 : p));
+ y = x - b * (b \ x);
+
+ if (m == 1)
+ y = y.';
+ end
+
endfunction
+//
+//input validation:
+//assert_checkerror("detrend1()", "detrend1: wrong number of input arguments");
+//a = "string";
+//assert_checkerror("detrend1(a)", "detrend1: X must be a numeric vector or matrix");
+//assert_checkerror("detrend1(%T)", "detrend1: X must be a numeric vector or matrix");
+//assert_checkerror("detrend1(1, -1)", "detrend1: P must be constant, linear, or a positive integer");
+//assert_checkerror("detrend1(1, 1.25)", "detrend1: P must be constant, linear, or a positive integer");
+
+//tests:
+//N = 32;
+//x = (0:1:N-1)/N + 2;
+//y = detrend1(x);
+//assert_checktrue(abs (y(:)) < 20*%eps);
+//assert_checkequal(detrend1([2, 5, 8]), detrend1([2. 5, 8], "linear"));
+//assert_checkequal(detrend1([2, 5, 8], 0), detrend1([2. 5, 8], "constant"));
+//assert_checkalmostequal(detrend1([1; 6; 9], "constant"), [-4.33333; 0.66666; 3.66666], 5*10^-5);
+//assert_checkalmostequal(detrend1([5, 12, 14; 8, 16, 14; 5, 10, 12]), [-1, -1.6666, -0.3333; 2, 3.3333, 0.6666; -1, -1.6666, -0.3333], 5*10^-4);
+//assert_checkalmostequal(detrend1([-5-5*%i; 2+%i; -4+3*%i], "linear"), [-2.1667-0.6667*%i; 4.3333+1.3333*%i; -2.1667-0.6667*%i], 5*10^-4);
+//assert_checkalmostequal(detrend1([5*%i, 1+2*%i,; -8, -1-6*%i], 0), [4+2.5*%i, 1+4*%i; -4-2.5*%i, -1-4*%i], 5*10^-4);
diff --git a/macros/dst1.sci b/macros/dst1.sci
index 7f2165f..c2f84ed 100644
--- a/macros/dst1.sci
+++ b/macros/dst1.sci
@@ -1,28 +1,62 @@
-function y = dst1(x, varargin)
+function y = dst1(x, n)
//Computes the type I discrete sine transform of x
-//Calling Sequence
+//Calling Sequence:
//y= dst1(x)
//y= dst1(x,n)
-//Parameters
+//Parameters:
//x: real or complex valued vector
//n= Length to which x is trimmed before transform
-//Description
-//This is an Octave function.
-//Computes the type I discrete sine transform of x. If n is given, then x is padded or trimmed to length n before computing the transform. If x is a matrix, compute the transform along the columns of the the matrix.
-funcprot(0);
-
-lhs= argn(1);
-rhs= argn(2);
-
-if(rhs>2)
-error("Wrong number of input arguments");
-end
-
-select(rhs)
- case 1 then
- y = callOctave("dst", x);
- case 2 then
- y = callOctave("dst", x, varargin(1));
-end
+//Description:
+//Computes the type 1 discrete sine transform of x. If n is given, then x is padded or trimmed to length n before computing the transform.
+//If x is a matrix, compute the transform along the columns of the the matrix.
+//Example:
+//dst1([1 3 6])
+//ans = [7.94974, -5, 1.94974]
+
+ funcprot(0);
+ rhs = argn(2);
+ if (rhs < 1 | rhs > 2)
+ error("ds1: wrong number of input arguments");
+ end
+
+ transpose = (size(x, 'r') == 1);
+ if (transpose)
+ x = x (:);
+ end
+
+ [nr, nc] = size(x);
+ if (rhs == 1)
+ n = nr;
+ elseif (n > nr)
+ x = [ x ; zeros(n-nr, nc) ];
+ elseif (n < nr)
+ x (nr-n+1 : n, :) = [];
+ end
+
+ d = [ zeros(1, nc); x; zeros(1, nc); -flipdim(x, 1) ];
+ y = fft(d, -1, find(size(d) ~= 1, 1))/2*%i;
+ y = y(2:nr+1, :);
+ if (isreal(x))
+ y = real (y);
+ end
+
+ if (transpose)
+ y = y.';
+ end
+
endfunction
+//input validation:
+//assert_checkerror("dst1()", "ds1: wrong number of input arguments");
+//assert_checkerror("dst1(1, 2, 3)", "Wrong number of input arguments.");
+
+//tests:
+//x = log(linspace(0.1,1,32));
+//y = dst1(x);
+//assert_checkalmostequal(y(3), sum(x.*sin(3*%pi*[1:32]/33)), 100*%eps);
+//
+//assert_checkalmostequal(dst1([-1; -3; -6]), [-7.9497; 5.0000; -1.9497], 5*10^-5);
+//assert_checkalmostequal(dst1([-1 2 2], 5), [3.2321, 0.8660, -3.0000], 5*10^-5);
+//assert_checkalmostequal(dst1([1+2*%i, 5+3*%i, 8+2*%i]), [11.3639+5.8284*%i, -7.0000-0*%i, 1.3640-0.1716*%i], 5*10^-4);
+//assert_checkalmostequal(dst1([1+2*%i; 5+3*%i; 8+2*%i], 2), [7.7942+3.4641*%i; -6.0622-0*%i; 0], 5*10^-5);
+//assert_checkalmostequal(dst1([-1-3*%i, 4*%i; -2-7*%i, 3]), [-2.5981-8.6603*%i, 2.5981+3.4641*%i; 0.8660+3.4641*%i, -2.5981+3.4641*%i], 5*10^-5);
diff --git a/macros/durbinlevinson.sci b/macros/durbinlevinson.sci
index 74dbb46..9fb647b 100644
--- a/macros/durbinlevinson.sci
+++ b/macros/durbinlevinson.sci
@@ -1,27 +1,88 @@
-function y= durbinlevinson(C, varargin)
-// Perform one step of the Durbin-Levinson algorithm..
+function [newphi, newv] = durbinlevinson (c, oldphi, oldv)
+// Perform one step of the Durbin-Levinson algorithm.
//Calling Sequence
-// durbinlevinson (C);
-// durbinlevinson (C, OLDPHI);
-// durbinlevinson (C, OLDPHI, OLDV);
-//Parameters
-//C: The vector C specifies the autocovariances '[gamma_0, ..., gamma_t]' from lag 0 to T.
-//OLDPHI: It specifies the coefficients based on C(T-1).
-//OLDV: It specifies the corresponding error.
-//Description
-//This is an Octave function.
+// durbinlevinson(c)
+// durbinlevinson(c, oldphi)
+// durbinlevinson(c, oldphi, oldv)
+//Parameters:
+//c: The vector c specifies the autocovariances '[gamma_0, ..., gamma_t]' from lag 0 to t.
+//oldphi: It specifies the coefficients based on c(t-1).
+//oldv: It specifies the corresponding error.
+//Description:
//Perform one step of the Durbin-Levinson.
-//If OLDPHI and OLDV are omitted, all steps from 1 to T of the algorithm are performed.
- rhs=argn(2);
- if(rhs<1 | rhs>3)
- error("Wrong number of input arguments");
- end
- select(rhs)
- case 1 then
- y=callOctave("durbinlevinson",C);
- case 2 then
- y=callOctave("durbinlevinson",C, varargin(1));
- case 3 then
- y=callOctave("durbinlevinson",C, varargin(1), varargin(2));
- end
-endfunction \ No newline at end of file
+//If 'oldphi' and 'oldv' are omitted, all steps from 1 to t of the algorithm are performed.
+//Example:
+//durbinlevinson([1, 2, 3], 1, 2)
+//ans = [0.5, 0.5]
+
+ funcprot(0);
+ rhs = argn(2);
+
+ if (rhs ~= 1 & rhs ~= 3)
+ error("durbinlevinson: wrong number of input arguments");
+ end
+
+ if (size(c, 'c') > 1)
+ c = c';
+ end
+
+ newphi = 0;
+ newv = 0;
+
+ if (rhs == 3)
+ t = length (oldphi) + 1;
+ if (length (c) < t + 1)
+ error ("durbinlevinson: arg1 (c) is too small");
+ end
+
+ if (oldv == 0)
+ error ("durbinlevinson: arg3 (oldv) must be non-zero");
+ end
+
+ if (size(oldphi, 'r') > 1)
+ oldphi = oldphi';
+ end
+
+ newphi = zeros (1, t);
+ newphi(1) = (c(t+1) - oldphi * c(2:t)) / oldv;
+ for i = 2 : t
+ newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1);
+ end
+ newv = (1 - newphi(1)^2) * oldv;
+
+ else
+ tt = length (c)-1;
+ oldphi = c(2) / c(1);
+ oldv = (1 - oldphi^2) * c(1);
+
+ for t = 2 : tt
+ newphi = zeros (1, t);
+ newphi(1) = (c(t+1) - oldphi * c(2:t)) / oldv;
+ for i = 2 : t
+ newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1);
+ end
+ newv = (1 - newphi(1)^2) * oldv;
+ oldv = newv;
+ oldphi = newphi;
+ end
+
+ end
+endfunction
+
+//input validation:
+//assert_checkerror("durbinlevinson()", "durbinlevinson: wrong number of input arguments");
+//assert_checkerror("durbinlevinson(1, 2, 3, 4)", "Wrong number of input arguments.");
+//assert_checkerror("durbinlevinson([1, 2], [1, 2, 3], 5)", "durbinlevinson: arg1 (c) is too small");
+//assert_checkerror("durbinlevinson([1, 2, 3, 4], [1, 2], 0)", "durbinlevinson: arg3 (oldv) must be non-zero");
+
+//test mx1 input:
+//assert_checkequal(durbinlevinson([2, 4, 3]), durbinlevinson([2, 4, 3], 2, -6));
+//assert_checkequal(durbinlevinson([9, 12, 43, 55], 3, -4), [-1.75, 8.25]);
+//assert_checkequal(durbinlevinson([2, 5, 3, 5, 7], [2, 6], -4), [5.75, -32.5, -5.5]);
+//assert_checkequal(durbinlevinson([6, 5, 8, 4], [1; 2], -1), [17, -33, -15]);
+
+//test 1xm input:
+//assert_checkequal(durbinlevinson([1; 4; 7]), [0.6, 1.6]);
+//assert_checkequal(durbinlevinson([3; 6; 7], -5, -2), [-18.5, -97.5]);
+//assert_checkequal(durbinlevinson([3; 6; 4; 2; 1; 9], [4; 6; 5], 3), [-19, 99, 120, 81]);
+//assert_checkequal(durbinlevinson([6; 5; 8; 4], [1, 2], -1), [17, -33, -15]);
diff --git a/macros/fftshift1.sci b/macros/fftshift1.sci
index cdff99a..db22891 100644
--- a/macros/fftshift1.sci
+++ b/macros/fftshift1.sci
@@ -1,32 +1,81 @@
-function y= fftshift1(X,DIM)
-//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)
-//Parameters
-//X:It is a vector of N elements corresponding to time samples
-//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 '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 (fft1 (X))' corresponds to frequencies
-//
-//f = [ -(ceil((N-1)/2):-1:1)*df 0 (1:floor((N-1)/2))*df ]
-//
-//where df = 1 / dt.
-//
-//If X is a matrix, the same holds for rows and columns. If X is an array, then the same holds along each dimension.
-//
-//The optional DIM argument can be used to limit the dimension along which the permutation occurs.
- rhs= argn(2);
- if(rhs <1 | rhs >2)
- error('Wrong number of Input arguments');
- end
- select(rhs)
- case 1 then
- y=callOctave("fftshift",X);
- case 2 then
- y=callOctave("fftshift",X,DIM);
- end
+function y = fftshift1(x, dim)
+//Performs a shift of the vector x, for use with the fft1 and ifft1 functions, in order to move the frequency 0 to the center of the vector or matrix.
+//Calling Sequence:
+// fftshift1(x)
+// fftshift1(x, dim)
+//Parameters:
+//x- It is a vector of N elements corresponding to time samples
+//dim- The optional DIM argument can be used to limit the dimension along which the permutation occurs
+//Examples:
+//x = [0:6]
+//fftshift1(x)
+//ans =
+//4 5 6 0 1 2 3
+
+ funcprot(0);
+ rhs = argn(2);
+ if (rhs < 1 | rhs > 2) then
+ error ("fftshift1: wrong number of arguments");
+ end
+
+ if (~or(type(x) == [1 5 8 10 4 6]))
+ error ("fftshitft1: arg1 (x) must be a vector or matrix");
+ end
+
+ if (rhs == 1) then
+ dim = 1:max(size(size(x)));
+ else
+ if (~(isscalar(dim) & dim > 0 & dim == fix(dim))) then
+ error ("fftshift1: arg2 (dim) must be a positive integer");
+ end
+ end
+
+ for d = dim
+ sz = size(x);
+ sz2 = ceil(sz(d) / 2);
+ dim_idx = [sz2+1:sz(d), 1:sz2];
+ if (d == 1) then
+ x = x(dim_idx, :);
+ elseif ( d == max(size(size(x))) ) then
+ x = x(:, dim_idx);
+ else
+ idx = repmat({':'}, 1, max(size(size(x))));
+ idx(d) = {dim_idx};
+ x = x(idx{:});
+ end
+ end
+ y = x;
+
endfunction
+
+//input validation:
+//assert_checkerror("fftshift1()", "fftshift1: wrong number of arguments");
+//assert_checkerror("fftshift1(1, 2, 3)", "Wrong number of input arguments.");
+//assert_checkerror("fftshift1(0:2, -1)", "fftshift1: arg2 (dim) must be a positive integer");
+//assert_checkerror("fftshift1(0:2, 0:3)", "fftshift1: arg2 (dim) must be a positive integer");
+
+//test mx1 input:
+//x = [0:7];
+//y = fftshift1 (x);
+//assert_checkequal (y, [4 5 6 7 0 1 2 3]);
+//assert_checkequal (fftshift1(y), x);
+
+//test 1xm input:
+//x = [0:7]';
+//y = fftshift1(x);
+//assert_checkequal(y, [4;5;6;7;0;1;2;3]);
+//assert_checkequal(fftshift1(y), x);
+
+//test mxn input:
+//x = [0:3];
+//x = [x;2*x;3*x+1;4*x+1];
+//y = fftshift1(x);
+//assert_checkequal(y, [[7 10 1 4];[9 13 1 5];[2 3 0 1];[4 6 0 2]]);
+//assert_checkequal(fftshift1(y), x);
+
+//test dim is provided:
+//x = [0:3];
+//x = [x;2*x;3*x+1;4*x+1];
+//y = fftshift1(x, 1);
+//assert_checkequal(y, [[1 4 7 10];[1 5 9 13];[0 1 2 3];[0 2 4 6]]);
+//assert_checkequal(fftshift1(y, 1), x);
diff --git a/macros/filter2.sci b/macros/filter2.sci
index b8b48e4..f565e14 100644
--- a/macros/filter2.sci
+++ b/macros/filter2.sci
@@ -1,33 +1,70 @@
-function Y = filter2 (B, X, SHAPE)
-//Apply the 2-D FIR filter B to X.
-//Calling Sequence
-//Y = filter2(B, X)
-//Y = filter2(B, X, SHAPE)
-//Parameters
-//B, X: Matrix
-// SHAPE:
-// 'full': pad X with zeros on all sides before filtering.
-// 'same': unpadded X (default)
-// 'valid': trim X after filtering so edge effects are no included.
-//Description
-//This function applies the 2-D FIR filter B to X. If the argument SHAPE is specified, return an array of the desired shape.
-//Examples
+function y = filter2 (b, x, shape)
+//Apply the 2-D FIR filter b to x.
+//Calling Sequence:
+//y = filter2(b, x)
+//y = filter2(b, x, shape)
+//Parameters:
+//b, x: vectors
+//If the optional argument 'shape' is specified, return a vector of the desired shape.
+//Possible values are:
+// "full"- pad x with zeros on all sides before filtering.
+// "same"- unpadded x (default)
+// "valid"- trim x after filtering so edge effects are no included.
+//Examples:
//filter2([1,3], [4,5])
-//ans =
-// 19. 5.
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 2 | rhs > 3)
-error("Wrong number of input arguments.")
-end
+//ans =
+//[19, 5]
-select(rhs)
-
- case 2 then
- Y=callOctave("filter2",B,X)
- case 3 then
- Y = callOctave("filter2",B,X,SHAPE)
+ funcprot(0);
+ rhs = argn(2);
+
+ if (rhs < 2 | rhs > 3) then
+ error("filter2: wrong number of input arguments");
+ end
+ if (rhs < 3) then
+ shape = "same";
+ end
+ [nr, nc] = size(b);
+ A = x;
+ B = b(nr:-1:1, nc:-1:1);
+ fullConv = convol2d(A, B);
+ sizeA = size(A);
+ sizeB = size(B);
+ sizeFull = size(fullConv);
+
+ select(shape)
+ case "full"
+ y = fullConv;
+ case "same"
+ startRow = floor(sizeB(1) / 2) + 1;
+ startCol = floor(sizeB(2) / 2) + 1;
+ endRow = startRow + sizeA(1) - 1;
+ endCol = startCol + sizeA(2) - 1;
+ y = fullConv(startRow:endRow, startCol:endCol);
+ case "valid"
+ startRow = sizeB(1);
+ startCol = sizeB(2);
+ endRow = sizeFull(1) - sizeB(1) + 1;
+ endCol = sizeFull(2) - sizeB(2) + 1;
+ if endRow >= startRow & endCol >= startCol then
+ y = fullConv(startRow:endRow, startCol:endCol);
+ else
+ y = [];
+ end
+ else
+ error("Invalid shape parameter.");
end
endfunction
+
+//input validation:
+//assert_checkerror("filter2()", "filter2: wrong number of input arguments");
+//assert_checkerror("filter2(1, 2, 3, 4)", "Wrong number of input arguments.");
+
+//tests:
+//assert_checkequal(filter2([1, 2; 3, 5], [1, 3; 5, 7]), [57, 24; 19, 7]);
+//assert_checkequal(filter2(1, 5), filter2(1, 5, "same"));
+//assert_checkequal(filter2([1, 2], [4; 5; 6], "full"), [8, 4; 10, 5; 12, 6]);
+//assert_checkequal(filter2([3; 1], [5; 4], "valid"), 19);
+//x=1;
+//assert_checkequal(filter2([3 x*2; 3*x+1 9^x], [5 7; 4 8], "valid"), 117);
diff --git a/macros/fwhm.sci b/macros/fwhm.sci
index 732f623..bb9c186 100644
--- a/macros/fwhm.sci
+++ b/macros/fwhm.sci
@@ -1,43 +1,148 @@
-function [f] = fwhm(y, varargin)
+function myfwhm = fwhm(y, varargin)
//This function computes peak full width at half minimum or at another level of peak minimum for vector or matrix data y supplied as input.
-//Calling Sequence
+//Calling Sequence:
//f = fwhm (y)
//f = fwhm (x, y)
//f = fwhm (…, "zero")
//f = fwhm (…, "min")
//f = fwhm (…, "alevel", level)
//f = fwhm (…, "rlevel", level)
-//Parameters
-//y: vector or matrix
-
-//Description
-//This is an Octave function.
-//This function computes peak full width at half minimum or at another level of peak minimum for vector or matrix data y supplied as input.
-//If y is a matrix, fwhm is calculated for each column as a row vector.
+//Parameters:
+//y- vector or matrix. If y is a matrix, fwhm is calculated for each column as a row vector.
//The second argument is by default "zero" which computes the fwhm at half maximum. If it is "min", fwhm is computed at middle curve.
//The option "rlevel" computes full-width at the given relative level of peak profile.
//The option "alevel" computes full-width at the given absolute level of y.
+//Description:
+//This function computes peak full width at half minimum or at another level of peak minimum for vector or matrix data y supplied as input.
//This function returns 0 if FWHM does not exist.
-//Examples
+//Examples:
//fwhm([1,2,3;9,-7,0.6],[4,5,6])
-//ans = 0.
-funcprot(0);
-
-rhs = argn(2)
-if(rhs<1 | rhs>5)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
- case 1 then
- f = callOctave("fwhm",y)
- case 2 then
- f = callOctave("fwhm",y,varargin(1))
- case 3 then
- f = callOctave("fwhm",y,varargin(1),varargin(2))
- case 4 then
- f = callOctave("fwhm",y,varargin(1),varargin(2),varargin(3))
- case 5 then
- f = callOctave("fwhm",y,varargin(1),varargin(2),varargin(3),varargin(4))
- end
+//ans = 0;
+
+ funcprot(0);
+ if nargin < 1 || nargin > 5
+ error("Wrong number of input arguments.");
+ end
+ opt = 'zero';
+ is_alevel = 0;
+ level = 0.5;
+ if nargin == 1
+ x = 1:max(size(y));
+ else
+ if type(varargin(1)) == 10
+ x = 1:max(size(y));
+ k = 1;
+ else
+ x = y;
+ y = varargin(1);
+ k = 2;
+ end
+ while k <= max(size(varargin))
+ if ~strcmp(varargin(k), 'alevel')
+ is_alevel = 1;
+ k = k+1;
+ if k > max(size(varargin))
+ error('option alevel requires an argument');
+ end
+ level = varargin(k);
+ if ~isreal(level) || max(size(level)) > 1
+ error('argument of alevel must be real number');
+ end
+ k = k+1;
+ break
+ end
+ if (~strcmp(varargin(k), 'zero') || ~strcmp(varargin(k), 'min'))
+ opt = varargin(k);
+ k = k+1;
+ end
+ if k > max(size(varargin)) break; end
+ if ~strcmp(varargin(k), 'rlevel')
+ k = k+1;
+ if k > max(size(varargin))
+ error('option rlevel requires an argument');
+ end
+ level = varargin(k);
+ if ~isreal(level) || max(size(level)) > 1 || level(1) < 0 || level(:) > 1
+ error('argument of rlevel must be real number from 0 to 1 (it is 0.5 for fwhm)');
+ end
+ k = k+1;
+ break
+ end
+ break
+ end
+ if k ~= max(size(varargin))+1
+ error('fwhm: extraneous option(s)');
+ end
+ end
+
+ // test the y matrix
+ [nr, nc] = size(y);
+ if (nr == 1 & nc > 1) then
+ y = y';
+ nr = nc;
+ nc = 1;
+ end
+
+ if max(size(x)) ~= nr then
+ error('dimension of input arguments do not match');
+ end
+
+ // Shift matrix columns so that y(+-xfwhm) = 0:
+ if is_alevel then
+ // case: full-width at the given absolute position
+ y = y - level;
+ else
+ if opt == 'zero' then
+ // case: full-width at half maximum
+ y = y - level * repmat(mtlb_max(y), nr, 1);
+ else
+ // case: full-width above background
+ y = y - level * repmat((mtlb_max(y) + mtlb_min(y)), nr, 1);
+ end
+ end
+
+ myfwhm = zeros(1, nc); // default: 0 for fwhm undefined
+ for n = 1:nc
+ yy = y(:, n);
+ ind = find((yy(1:$ - 1) .* yy(2:$)) <= 0);
+ if mtlb_max(size(ind)) >= 2 & yy(ind(1)) > 0 then // must start ascending
+ ind = ind(2:$);
+ end
+ [mx, imax] = mtlb_max(yy); // protection against constant or (almost) monotonous functions
+ if max(size(ind)) >= 2 & imax >= ind(1) & imax <= ind($) then
+ ind1 = ind(1);
+ ind2 = ind1 + 1;
+ xx1 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
+ ind1 = ind($);
+ ind2 = ind1 + 1;
+ xx2 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
+ myfwhm(n) = xx2 - xx1;
+ end
+ end
endfunction
+
+//tests:
+//x = -%pi:0.001:%pi; y = cos(x);
+//assert_checkalmostequal(fwhm(x, y), 2.094395, 0.5*10^-5);
+//assert_checktrue( abs(fwhm(x, y) - 2*%pi/3) < 0.01 );
+
+//assert_checktrue(fwhm(-10:10) == 0 && fwhm(ones(1,50)) == 0);
+
+//x=1:3;
+//y=[-1,3,-1];
+//assert_checktrue(abs(fwhm(x,y)-0.75)<0.001 && abs(fwhm(x,y,'zero')-0.75) < 0.001 && abs(fwhm(x,y,'min')-1.0) < 0.001);
+
+//x=1:3;
+//y=[-1,3,-1];
+//assert_checktrue(abs(fwhm(x,y, 'rlevel', 0.1)-1.35) < 0.001 && abs(fwhm(x,y,'zero', 'rlevel', 0.1)-1.35) < 0.001 && abs(fwhm(x,y,'min', 'rlevel', 0.1)-1.40) < 0.001);
+
+//x=1:3;
+//y=[-1,3,-1];
+//assert_checktrue(abs(fwhm(x,y, 'alevel', 2.5)-0.25) < 0.001 && abs(fwhm(x,y,'alevel', -0.5)-1.75) < 0.001);
+
+//x=-10:10;
+//assert_checktrue(fwhm(x.*x) == 0);
+
+//x=-5:5;
+//y=18-x.*x;
+//assert_checktrue( abs(fwhm(y)-6.0) < 0.001 && abs(fwhm(x,y,'zero')-6.0) < 0.001 && abs(fwhm(x,y,'min')-7.0 ) < 0.001);
diff --git a/macros/fwhmjlt.sci b/macros/fwhmjlt.sci
index 129bffc..0821930 100644
--- a/macros/fwhmjlt.sci
+++ b/macros/fwhmjlt.sci
@@ -1,40 +1,145 @@
-function [f]=fwhmjlt(y,varargin)
-//This function Computes peak full-width at half maximum
-
-//calling sequence
-//f = fwhm (y)
-//f = fwhm (x, y)
-//f = fwhm (…, "zero")
-//f = fwhm (…, "min")
-//f = fwhm (…, "alevel", level)
-//f = fwhm (…, "rlevel", level)
-
-//Description
-//Compute peak full-width at half maximum (FWHM) or at another level of peak maximum for vector or matrix data y, optionally sampled as y(x). If y is a matrix, return FWHM for each column as a row vector.
-//The default option "zero" computes fwhm at half maximum, i.e. 0.5*max(y). The option "min" computes fwhm at the middle curve, i.e. 0.5*(min(y)+max(y)).
-//The option "rlevel" computes full-width at the given relative level of peak profile
+function myfwhm = fwhmjlt(y, varargin)
+//This function computes peak full width at half minimum or at another level of peak minimum for vector or matrix data y supplied as input.
+//Calling Sequence:
+//f = fwhmjlt(y)
+//f = fwhmjlt(x, y)
+//f = fwhmjlt(…, "zero")
+//f = fwhmjlt(…, "min")
+//f = fwhmjlt(…, "alevel", level)
+//f = fwhmjlt(…, "rlevel", level)
+//Parameters:
+//y- vector or matrix. If y is a matrix, fwhm is calculated for each column as a row vector.
+//The second argument is by default "zero" which computes the fwhm at half maximum. If it is "min", fwhm is computed at middle curve.
+//The option "rlevel" computes full-width at the given relative level of peak profile.
//The option "alevel" computes full-width at the given absolute level of y.
-
-//Example
+//Description:
+//This function computes peak full width at half minimum or at another level of peak minimum for vector or matrix data y supplied as input.
+//This function returns 0 if FWHM does not exist.
+//Examples:
//t=-50:0.01:50;
//y=(1/(2*sqrt(2*%pi)))*exp(-(t.^2)/8);
//z=fwhmjlt(y)
//Output: 470.96442
-rhs = argn(2)
-if(rhs<1 | rhs>5)
-error("Wrong number of input arguments.")
-end
- select(rhs)
- case 1 then
- f = callOctave("fwhm",y)
- case 2 then
- f = callOctave("fwhm",y,varargin(1))
- case 3 then
- f = callOctave("fwhm",y,varargin(1),varargin(2))
- case 4 then
- f = callOctave("fwhm",y,varargin(1),varargin(2),varargin(3))
- case 5 then
- f = callOctave("fwhm",y,varargin(1),varargin(2),varargin(3),varargin(4))
- end
+ funcprot(0);
+ if nargin < 1 || nargin > 5
+ error("Wrong number of input arguments.");
+ end
+ opt = 'zero';
+ is_alevel = 0;
+ level = 0.5;
+ if nargin == 1
+ x = 1:max(size(y));
+ else
+ if type(varargin(1)) == 10
+ x = 1:max(size(y));
+ k = 1;
+ else
+ x = y;
+ y = varargin(1);
+ k = 2;
+ end
+ while k <= max(size(varargin))
+ if ~strcmp(varargin(k), 'alevel')
+ is_alevel = 1;
+ k = k+1;
+ if k > max(size(varargin))
+ error('option alevel requires an argument');
+ end
+ level = varargin(k);
+ if ~isreal(level) || max(size(level)) > 1
+ error('argument of alevel must be real number');
+ end
+ k = k+1;
+ break
+ end
+ if (~strcmp(varargin(k), 'zero') || ~strcmp(varargin(k), 'min'))
+ opt = varargin(k);
+ k = k+1;
+ end
+ if k > max(size(varargin)) break; end
+ if ~strcmp(varargin(k), 'rlevel')
+ k = k+1;
+ if k > max(size(varargin))
+ error('option rlevel requires an argument');
+ end
+ level = varargin(k);
+ if ~isreal(level) || max(size(level)) > 1 || level(1) < 0 || level(:) > 1
+ error('argument of rlevel must be real number from 0 to 1 (it is 0.5 for fwhm)');
+ end
+ k = k+1;
+ break
+ end
+ break
+ end
+ if k ~= max(size(varargin))+1
+ error('fwhmjlt: extraneous option(s)');
+ end
+ end
+
+ [nr, nc] = size(y);
+ if (nr == 1 & nc > 1) then
+ y = y';
+ nr = nc;
+ nc = 1;
+ end
+
+ if max(size(x)) ~= nr then
+ error('dimension of input arguments do not match');
+ end
+
+ if is_alevel then
+ y = y - level;
+ else
+ if opt == 'zero' then
+ y = y - level * repmat(mtlb_max(y), nr, 1);
+ else
+ y = y - level * repmat((mtlb_max(y) + mtlb_min(y)), nr, 1);
+ end
+ end
+
+ myfwhm = zeros(1, nc);
+ for n = 1:nc
+ yy = y(:, n);
+ ind = find((yy(1:$ - 1) .* yy(2:$)) <= 0);
+ if mtlb_max(size(ind)) >= 2 & yy(ind(1)) > 0
+ ind = ind(2:$);
+ end
+ [mx, imax] = mtlb_max(yy);
+ if max(size(ind)) >= 2 & imax >= ind(1) & imax <= ind($) then
+ ind1 = ind(1);
+ ind2 = ind1 + 1;
+ xx1 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
+ ind1 = ind($);
+ ind2 = ind1 + 1;
+ xx2 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
+ myfwhm(n) = xx2 - xx1;
+ end
+ end
endfunction
+
+//tests:
+//x = -%pi:0.001:%pi; y = cos(x);
+//assert_checkalmostequal(fwhmjlt(x, y), 2.094395, 0.5*10^-5);
+//assert_checktrue( abs(fwhmjlt(x, y) - 2*%pi/3) < 0.01 );
+
+//assert_checktrue(fwhmjlt(-10:10) == 0 && fwhmjlt(ones(1,50)) == 0);
+
+//x=1:3;
+//y=[-1,3,-1];
+//assert_checktrue(abs(fwhmjlt(x,y)-0.75)<0.001 && abs(fwhmjlt(x,y,'zero')-0.75) < 0.001 && abs(fwhmjlt(x,y,'min')-1.0) < 0.001);
+
+//x=1:3;
+//y=[-1,3,-1];
+//assert_checktrue(abs(fwhmjlt(x,y, 'rlevel', 0.1)-1.35) < 0.001 && abs(fwhmjlt(x,y,'zero', 'rlevel', 0.1)-1.35) < 0.001 && abs(fwhmjlt(x,y,'min', 'rlevel', 0.1)-1.40) < 0.001);
+
+//x=1:3;
+//y=[-1,3,-1];
+//assert_checktrue(abs(fwhmjlt(x,y, 'alevel', 2.5)-0.25) < 0.001 && abs(fwhmjlt(x,y,'alevel', -0.5)-1.75) < 0.001);
+
+//x=-10:10;
+//assert_checktrue(fwhmjlt(x.*x) == 0);
+
+//x=-5:5;
+//y=18-x.*x;
+//assert_checktrue( abs(fwhmjlt(y)-6.0) < 0.001 && abs(fwhmjlt(x,y,'zero')-6.0) < 0.001 && abs(fwhmjlt(x,y,'min')-7.0 ) < 0.001);
diff --git a/macros/hurst.sci b/macros/hurst.sci
index 27507fb..bd48ee7 100644
--- a/macros/hurst.sci
+++ b/macros/hurst.sci
@@ -1,21 +1,47 @@
-function y = hurst(x)
-// Estimate the Hurst parameter of sample X via the rescaled r statistic.
-//Calling Sequence
-//hurst(X)
-//variable=hurst(X)
-//Parameters
-//X: X is a matrix, the parameter of sample X via the rescaled r statistic
-//Description
-//This is an Octave function.
-//This function estimates the Hurst parameter of sample X via the rescaled rstatistic.
-funcprot(0);
-rhs= argn(2);
-if(rhs<1 | rhs>1)
- error("Wrong number of input arguments");
-end
+function H = hurst(x)
+//Estimate the Hurst parameter of sample X via the rescaled range statistic.
+//Calling Sequence:
+//hurst(x)
+//Parameters:
+//x: x is a vector or matrix
+//Description:
+//This function estimates the Hurst parameter of sample x using the rescaled range statistic.
+//If x is a matrix, the parameter is estimated for every column.
+//Examples:
+//hurst([10, 15, 3])
+//ans = 0.045019
+
+ funcprot(0);
+ if (argn(2) ~= 1)
+ error("hurst: wrong number of input arguments");
+ end
+
+ if (isscalar(x))
+ error ("hurst: argument must not be a scalar");
+ end
+ if (isvector(x))
+ x = x(:);
+ end
+
+ [xr, xc] = size (x);
+
+ s = stdev(x, 'r');
+ y = [];
+ for i = 1:xr
+ y(i, :) = x(i, :) - mean(x, 'm');
+ end
+ w = cumsum(y);
+ RS = (max(w, 'r') - min(w, 'r')) ./ s;
+ H = log(RS) / log(xr);
-select(rhs)
- case 1 then
- y= callOctave("hurst", x);
-end
endfunction
+
+//input validation:
+//assert_checkerror("hurst()", "hurst: wrong number of input arguments");
+//assert_checkerror("hurst(1, 2, 3)", "Wrong number of input arguments.");
+//assert_checkerror("hurst(1)", "hurst: argument must not be a scalar");
+
+//tests:
+//assert_checkalmostequal(hurst([1, 5, 7, 14, 6]), 0.31180, 5*10^-5);
+//assert_checkalmostequal(hurst([3; 6; 9; 5]), 0.24271, 5*10^-5);
+//assert_checkalmostequal(hurst([-1, 9, 4; 7, 4, -3; 6, 12, -18]),[ 0.124902, 0.063474, 0.084510], 5*10^-5);
diff --git a/macros/ifftshift1.sci b/macros/ifftshift1.sci
index 7426130..c81a285 100644
--- a/macros/ifftshift1.sci
+++ b/macros/ifftshift1.sci
@@ -1,24 +1,83 @@
-function y= ifftshift1(X,DIM)
+function y = ifftshift1(x, dim)
//Undo the action of the 'fftshift1' function.
-//Calling Sequence
-// ifftshift1 (X)
-// ifftshift1 (X, DIM)
-//Parameters
-//X:It is a vector of N elements corresponding to time samples
-//DIM: The optional DIM argument can be used to limit the dimension along which the permutation occurs
-//Description
-//This is an Octave function.
-//Undo the action of the 'fftshift1' function.
-//
-//For even length X, 'fftshift1' is its own inverse, but odd lengths differ slightly.
- rhs= argn(2);
- if(rhs <1 | rhs >2)
- error('Wrong number of Input arguments');
- end
- select(rhs)
- case 1 then
- y=callOctave("ifftshift",X);
- case 2 then
- y=callOctave("ifftshift",X,DIM);
- end
-endfunction \ No newline at end of file
+//Calling Sequence:
+// ifftshift1(x)
+// ifftshift1(x, dim)
+//Parameters:
+//x- It is a vector of N elements corresponding to time samples
+//dim- The optional DIM argument can be used to limit the dimension along which the permutation occurs
+//Description:
+//Undoes the action of the 'fftshift1' function. For 'x' of even length 'fftshift1' is its own inverse, but odd lengths differ slightly.
+//Examples:
+//x = [1, 2, 3, 4];
+//ifftshift1(fftshift1(x));
+//ans =
+//[1, 2, 3, 4];
+
+ funcprot(0);
+ rhs = argn(2);
+ if (rhs < 1 | rhs > 2) then
+ error ("ifftshift1: wrong number of arguments");
+ end
+
+ if (~or(type(x) == [1 5 8 10 4 6]))
+ error ("ifftshitft1: arg1 (x) must be a vector or matrix");
+ end
+
+ if (rhs == 1) then
+ dim = 1:max(size(size(x)));
+ else
+ if (~(isscalar(dim) & dim > 0 & dim == fix(dim))) then
+ error ("ifftshift1: arg2 (dim) must be a positive integer");
+ end
+ end
+
+ for d = dim
+ sz = size(x);
+ sz2 = floor(sz(d) / 2);
+ dim_idx = [sz2+1:sz(d), 1:sz2];
+ if (d == 1) then
+ x = x(dim_idx, :);
+ elseif ( d == max(size(size(x))) ) then
+ x = x(:, dim_idx);
+ else
+ idx = repmat({':'}, 1, max(size(size(x))));
+ idx(d) = {dim_idx};
+ x = x(idx{:});
+ end
+ end
+ y = x;
+
+endfunction
+
+//input validation:
+//assert_checkerror("ifftshift1()", "ifftshift1: wrong number of arguments");
+//assert_checkerror("ifftshift1(1, 2, 3)", "Wrong number of input arguments.");
+//assert_checkerror("ifftshift1(0:2, -1)", "ifftshift1: arg2 (dim) must be a positive integer");
+//assert_checkerror("ifftshift1(0:2, 0:3)", "ifftshift1: arg2 (dim) must be a positive integer");
+
+//test mx1 input:
+//x = [0:7];
+//y = ifftshift1(x);
+//assert_checkequal(y, [4 5 6 7 0 1 2 3]);
+//assert_checkequal(ifftshift1(y), x);
+
+//test 1xm input:
+//x = [0:6]';
+//y = ifftshift1(x);
+//assert_checkequal(y, [3;4;5;6;0;1;2]);
+//assert_checkequal(ifftshift1(y), [6;0;1;2;3;4;5]);
+
+//test mxn input:
+//x = [0:3];
+//x = [x;2*x;3*x+1;4*x+1];
+//y = ifftshift1(x);
+//assert_checkequal(y, [[7 10 1 4];[9 13 1 5];[2 3 0 1];[4 6 0 2]]);
+//assert_checkequal(ifftshift1(y), x);
+
+//test dim is provided:
+//x = [0:3];
+//x = [x;2*x;3*x+1;4*x+1];
+//y = ifftshift1(x, 2);
+//assert_checkequal(y, [[2 3 0 1];[4 6 0 2];[7 10 1 4];[9 13 1 5]]);
+//assert_checkequal(ifftshift1(y, 2), x);
diff --git a/macros/ifht.sci b/macros/ifht.sci
index 7f865cb..b1a7000 100644
--- a/macros/ifht.sci
+++ b/macros/ifht.sci
@@ -1,34 +1,64 @@
-function m = ifht(d, varargin)
-//Calculate the inverse Fast Hartley Transform of real input D
-//Calling Sequence
-//m= ifht (d)
-//m= ifht (d,n)
-//m= ifht (d,n,dim)
-//Parameters
-//d: real or complex valued scalar or vector
-//n: Similar to the options of FFT function
-//dim: Similar to the options of FFT function
-//Description
-//Calculate the inverse Fast Hartley Transform of real input d. If d is a matrix, the inverse Hartley transform is calculated along the columns by default. The options n and dim are similar to the options of FFT function.
-//
+function m = ifht(d, n, dim)
+//Calculate the inverse Fast Hartley Transform of real input D.
+//Calling Sequence:
+//m: ifht(d)
+//m: ifht(d,n)
+//m: ifht(d,n,dim)
+//Parameters:
+//d: Real or complex scalar or vector.
+//n: Integer. Specifies the number of elements of x to be used.
+//dim: Integer. Specifies the dimention of the matrix along which the ifht is performed.
+//Description:
+//Calculate the inverse Fast Hartley Transform of real input d. If d is a matrix, the inverse Hartley transform is calculated along the columns by default.
+//The options n and dim are similar to the options of FFT function. 'n' is an integer that determines the number of elements of 'd' to use. If 'n' is larger than the dimension along which the ifht is calculated, then 'd' is resized and padded with zeros to match the required size. If n < d, then 'd' is truncated to match the required size.
+//'dim' is an integer that specifies the dimension of the matrix along which the ifht is to be performed. By default, the ifht is computed along the first non-singleton dimension of the array.
//The forward and inverse Hartley transforms are the same (except for a scale factor of 1/N for the inverse hartley transform), but implemented using different functions.
-//
-//The definition of the forward hartley transform for vector d, m[K] = 1/N \sum_{i=0}^{N-1} d[i]*(cos[K*2*pi*i/N] + sin[K*2*pi*i/N]), for 0 <= K < N. m[K] = 1/N \sum_{i=0}^{N-1} d[i]*CAS[K*i], for 0 <= K < N.
-//Examples
+//The definition of the forward hartley transform for vector d, m[K] = 1/N \sum_{i=0}^{N-1} d[i]*(cos[K*2*pi*i/N] + sin[K*2*pi*i/N]),
+//for 0 <= K < N. m[K] = 1/N \sum_{i=0}^{N-1} d[i]*CAS[K*i], for 0 <= K < N.
+//Examples:
//ifht(1 : 4)
-//ifht(1:4, 2)
-funcprot(0);
-rhs= argn(2);
-if(rhs<1 | rhs>3)
-error("Wrong number of Inputs")
-end
+//ans = [2.5, -1, -0.5, 0]
+
+ funcprot(0);
+ rhs = argn(2);
+ if (rhs < 1 | rhs > 3)
+ error("Wrong number of input arguments.");
+ end
+ dimension = size(d);
+ nd = find(dimension ~= 1, 1);
+ if (rhs == 3)
+ if isempty(n)
+ y = fft(d, 1, dim);
+ else
+ dimension(dim)=n;
+ y=fft(resize_matrix(d, dimension), 1, dim);
+ end
+ elseif (rhs == 2)
+ if isempty(n)
+ y = fft(d, 1, nd);
+ else
+ dimension(nd) = n;
+ y=fft(resize_matrix(d, dimension), 1, nd)
+ end
+ else
+ y = fft(d, 1, nd);
+ end
+
+ m = real(y) + imag(y);
-select(rhs)
- case 1 then
- m= callOctave("ifht", d);
- case 2 then
- m= callOctave("ifht", d , varargin(1));
- case 3 then
- m= callOctave("ifht", d , varargin(1),varargin(2) );
-end
endfunction
+
+//input validation:
+//assert_checkerror("ifht()", "Wrong number of input arguments.");
+//assert_checkerror("ifht(1, 2, 3, 4)", "Wrong number of input arguments.");
+
+////tests:
+//assert_checkequal(ifht(1+2*%i), 3);
+//assert_checkequal(ifht((1:4)), [2.5, -1, -0.5, 0]);
+//assert_checkequal(ifht([1:4]', 2), [1.5; -0.5]);
+//assert_checkequal(ifht([1:4]', 2, 2), [0.5 0.5; 1 1; 1.5 1.5; 2 2]);
+//assert_checkequal(ifht([-1 2; 3 -5]), [1 -1.5; -2 3.5]);
+//assert_checkalmostequal(ifht([1:3; -2:-5]), [2, -0.7887, -0.2113], 5*10^-4);
+//assert_checkequal(ifht([1:3; -2:-5], 1, 1), [1:3]);
+//assert_checkequal(ifht([1+2*%i, 3*%i; -4-3*%i, -5*%i]), [-2 -1; 5 4]);
+//assert_checkequal(ifht([1+2*%i, 3*%i; -4-3*%i, -5*%i], 1, 2), [3; -7]);
diff --git a/macros/morlet.sci b/macros/morlet.sci
index 0a160f0..8bd25a1 100644
--- a/macros/morlet.sci
+++ b/macros/morlet.sci
@@ -1,27 +1,68 @@
-function [psi,x] = morlet (lb,ub,n)
-
-// Generates Morlet wavelets
-// Calling Sequence
-// [psi,x]= morlet(lb,ub,n)
-// Parameters
-// lb: Real or complex valued vector or matrix
-// ub: Real or complex valued vector or matrix
-// n: Real strictly positive scalar number
-// Description
-// This is an Octave function
-// This function returns values of the Morlet wavelet in the specified interval for all the sample points.
-// Examples
-// 1. [a,b]=morlet(1,2,3)
-// a = [0.17205 0.11254 -0.11356]
-// b = [1.0000 1.5000 2.0000]
-// 2. [a,b]=morlet([1 2 3],[1 2 3],1)
-// a = [0.1720498; -0.1135560; -0.0084394]
-// b = [1; 2; 3]
-
-funcprot(0);
-rhs=argn(2);
-if (rhs<3) then
- error ("Wrong number of input arguments.")
-else [psi,x] = callOctave("morlet",lb,ub,n)
-end
+function [psi,x] = morlet(lb, ub, n)
+//Compute the Morlet wavelet.
+//Calling sequence:
+//[psi,x]= morlet(lb,ub,n)
+//Parameters:
+//lb: Real or complex valued vector/scalar
+//ub: Real or complex valued vector/scalar
+//n: Real positive scalar number
+//Description:
+//This function returns values of the Morlet wavelet in the specified interval for all the sample points.
+//Example:
+//[a,b] = morlet([1 2 3], [1 2 3], 1)
+//a = [0.1720498; -0.1135560; -0.0084394]
+//b = [1; 2; 3]
+
+ funcprot(0);
+ rhs = argn(2);
+
+ if (rhs ~= 3) then
+ error("morlet: wrong number of input arguments");
+ end
+
+ if (length(lb) ~= length(ub)) then
+ error("morlet: arg1 and arg2 msut have same dimension");
+ end
+
+ if (~isreal(n) | n <= 0) then
+ error("morlet: n must be a strictly positive real number");
+ end
+ x = [];
+ for i=1:length(lb)
+ x(i, :) = linspace(lb(i), ub(i), n);
+ end
+ psi = cos(5.*x) .* exp(-x.^2/2);
+
endfunction
+
+//input validation:
+//assert_checkerror("morlet()", "morlet: wrong number of input arguments");
+//assert_checkerror("morlet(1, 2)", "morlet: wrong number of input arguments");
+//assert_checkerror("morlet(1, 2, 3, 4)", "Wrong number of input arguments.");
+//assert_checkerror("morlet(1, 2, -1)", "morlet: n must be a strictly positive real number");
+//assert_checkerror("morlet(1, 2, 2+3*%i)", "morlet: n must be a strictly positive real number");
+//assert_checkerror(" morlet([5, 2, 7], [1, 3], 3)", "morlet: arg1 and arg2 msut have same dimension");
+
+//test basic input:
+//[a, b] = morlet(1, 2, 3);
+//assert_checkalmostequal(a, [0.17205, 0.11254, -0.11356], 5e-5);
+//assert_checkalmostequal(b, [1, 1.5, 2], %eps);
+
+//test complex input:
+//[a, b] = morlet(3+2*%i, 4+8*%i, 2);
+//assert_checkalmostequal(a, [-495.15886-756.35443*%i, -5.08135E26-3.07588E27*%i], 5e-5);
+//assert_checkalmostequal(b, [3+2*%i, 4+8*%i], %eps);
+
+//test real vector:
+//[a, b] = morlet([1, 2], [3, 4], 2);
+//A = [0.1720498, -0.0084394; -0.113556, 0.0001369];
+//B = [1, 3; 2, 4];
+//assert_checkalmostequal(A, a, 5e-5);
+//assert_checkalmostequal(B, b, %eps);
+
+//test complex vector:
+//[a, b] = morlet([1 + 2*%i, 3*%i], [2*%i, 2 + 3*%i], 1);
+//A = [81377.39587; -19069291.16508 + 5732843.75676*%i];
+//B = [2*%i; 2+3*%i];
+//assert_checkalmostequal(a, A, 5e-5);
+//assert_checkalmostequal(b, B, %eps);
diff --git a/macros/pburg.sci b/macros/pburg.sci
index 78198ed..f677157 100644
--- a/macros/pburg.sci
+++ b/macros/pburg.sci
@@ -1,93 +1,138 @@
-function [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion)
+function [psd,f_out] = pburg(x, poles, varargin)
//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
+
+//Calling Sequence:
+//[psd,f_out] = pburg(x, poles, freq, Fs, range, method, plot_type, criterion)
+
+//Parameters:
+//All but the first two parameters are optional and may be empty.
+//
+//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]
+//
+//CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
+//Control-string arguments can be in any order after the other arguments.
+//
+//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. 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
+//
+// RETURNED VALUES:
+//If return values are not required by the caller, the spectrum is plotted and nothing is returned.
+//psd: [real vector] power-spectral density estimate.
+//f_out: [real vector] frequency values.
+
+//Description:
//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
+//Examples:
+//a = [1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746];
+//[psd, f_out] = pburg(a, 2);
-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
+ funcprot(0);
+ if (nargin < 2)
+ error('pburg: need at least 2 args.');
+ end
+ nvarargin = length(varargin);
+ criterion = [];
+ for iarg = 1:nvarargin
+ arrgh = varargin(iarg);
+ if (type(arrgh) == 10 && ( ~strcmp(arrgh,'AKICc') ||...
+ ~strcmp(arrgh,'KIC') || ~strcmp(arrgh,'AICc') ||...
+ ~strcmp(arrgh,'AIC') || ~strcmp(arrgh,'FPE') ) )
+ criterion = arrgh;
+ if (nvarargin > 1)
+ varargin(iarg) = [];
+ else
+ varargin = list();
+ end
+ end
+ end
+ [ar_coeffs,residual] = arburg(x,poles,criterion);
+ if (nargout == 0)
+ ar_psd(ar_coeffs,residual,varargin(:));
+ elseif (nargout == 1)
+ psd = ar_psd(ar_coeffs,residual,varargin(:));
+ elseif (nargout >= 2)
+ [psd,f_out] = ar_psd(ar_coeffs,residual,varargin(:));
+ 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
+//tests:
+
+//fs = 1000;
+//t = 0:1/fs:1-1/fs;
+//x = cos(2*%pi*100*t);
+//order = 4;
+//[pxx, f] = pburg(x, order, [], fs);
+//figure;
+//plot(f, 10*log10(pxx));
+//title('PSD Estimate using Burg Method - Sinusoidal Signal');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency (dB/Hz)');
+
+//fs = 1000;
+//t = 0:1/fs:1-1/fs;
+//x = cos(2*%pi*100*t);
+//orders = [2, 4, 8, 16];
+//figure;
+//for i = 1:length(orders)
+// order = orders(i);
+// [pxx, f] = pburg(x, order, [], fs);
+// subplot(length(orders), 1, i);
+// plot(f, 10*log10(pxx));
+// title(['PSD Estimate using Burg Method - Order ' string(order)]);
+// xlabel('Frequency (Hz)');
+// ylabel('Power/Frequency (dB/Hz)');
+//end
+
+//fs = 1000;
+//t = 0:1/fs:0.1-1/fs;
+//x = cos(2*%pi*100*t);
+//order = 4;
+//[pxx, f] = pburg(x, order, [], fs);
+//figure;
+//plot(f, 10*log10(pxx));
+//title('PSD Estimate using Burg Method - Short Data Segment');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency (dB/Hz)');
+//fs = 1000;
+//t = 0:1/fs:1-1/fs;
+//x = cos(2*%pi*100*t) + cos(2*%pi*200*t);
+//order = 4;
+//[pxx, f] = pburg(x, order, [], fs);
+//figure;
+//plot(f, 10*log10(pxx));
+//title('PSD Estimate using Burg Method - Multicomponent Signal');
+//xlabel('Frequency (Hz)');
+//ylabel('Power/Frequency (dB/Hz)');
diff --git a/macros/pei_tseng_notch.sci b/macros/pei_tseng_notch.sci
index 1de0fcd..fcb6b9c 100644
--- a/macros/pei_tseng_notch.sci
+++ b/macros/pei_tseng_notch.sci
@@ -1,44 +1,73 @@
-function [b, a] = pei_tseng_notch (frequencies, bandwidths)
-
+function [ b, a ] = pei_tseng_notch ( frequencies, bandwidths )
//Return coefficients for an IIR notch-filter.
-//Calling Sequence
-//[b, a] = pei_tseng_notch (frequencies, bandwidths)
-//b = pei_tseng_notch (frequencies, bandwidths)
-//Parameters
-//frequencies: filter frequencies
-//bandwidths: bandwidths to be used with filter
-//Description
-//This is an Octave function.
-//It return coefficients for an IIR notch-filter with one or more filter frequencies and according bandwidths. The filter is based on a all pass filter that performs phasereversal at filter frequencies.
-//This leads to removal of those frequencies of the original and phase-distorted signal.
-//Examples
+//Calling Sequence:
+//[b, a] = pei_tseng_notch(frequencies, bandwidths)
+//Parameters:
+//frequencies: Real scalar/vector representing filter frequencies.
+//bandwidths: Real scalar scalar/vector representing bandwidths to be used with filter.
+//Description:
+//THis function returns coefficients for an IIR notch-filter with one or more filter frequencies and according bandwidths.
+//The filter is based on a all pass filter that performs phasereversal at filter frequencies. This leads to removal of those frequencies of the original and phase-distorted signal.
+//Examples:
//sf = 800; sf2 = sf/2;
-//data=[[1;zeros(sf-1,1)],sinetone(49,sf,1,1),sinetone(50,sf,1,1),sinetone(51,sf,1,1)];
-//[b,a]=pei_tseng_notch ( 50 / sf2, 2/sf2 )
-//b =
-//
-// 0.99213 -1.83322 0.99213
-//
-//a =
-//
-// 1.00000 -1.83322 0.98426
-
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 2 | rhs > 2)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-
- case 2 then
- if(lhs==1)
- b = callOctave("pei_tseng_notch", frequencies, bandwidths)
- elseif(lhs==2)
- [b, a] = callOctave("pei_tseng_notch", frequencies, bandwidths)
- else
- error("Wrong number of output argments.")
- end
- end
+//data = [[1;zeros(sf-1,1)],sinetone(49,sf,1,1),sinetone(50,sf,1,1),sinetone(51,sf,1,1)];
+//[b,a] = pei_tseng_notch ( 50 / sf2, 2/sf2 )
+//b = 0.99213 -1.83322 0.99213
+//a = 1.00000 -1.83322 0.98426
+
+ if (nargin() ~= 2)
+ error("Wrong number of input arguments.");
+ elseif ( ~isvector(frequencies) | ~isvector(bandwidths) )
+ if ( ~isscalar(frequencies) | ~isscalar(bandwidths) )
+ error("All arguments must be vectors or scalars.");
+ end
+ elseif ( length(frequencies) ~= length (bandwidths) )
+ error("All arguments must be of equal length.");
+ elseif (~and( frequencies > 0 & frequencies < 1 ))
+ error("All frequencies must be in (0, 1).");
+ elseif (~and( bandwidths > 0 & bandwidths < 1 ))
+ error("All bandwidths must be in (0, 1).");
+ end
+
+ frequencies = frequencies (:)';
+ bandwidths = bandwidths (:)';
+
+ frequencies = frequencies*%pi;
+ bandwidths = bandwidths*%pi;
+ M2 = 2 * length(frequencies);
+
+ a = [ frequencies - bandwidths / 2; frequencies ];
+ omega = a(:);
+ factors = ( 1 : 2 : M2 );
+ b = [ -%pi * factors + %pi / 2; -%pi * factors ];
+ phi = b(:);
+ t_beta = tan ( ( phi + M2 * omega ) / 2 );
+
+ Q = zeros (M2, M2);
+
+ for k = 1 : M2
+ Q ( :, k ) = sin ( k .* omega ) - t_beta .* cos ( k .* omega );
+ end
+
+ h_a = ( Q \ t_beta )';
+ denom = [ 1, h_a ];
+ num = [ flipdim(h_a, 2), 1 ];
+
+ a = denom;
+ b = ( num + denom ) / 2;
+
endfunction
+
+//input validation:
+//assert_checkerror("pei_tseng_notch()", "Wrong number of input arguments.");
+//assert_checkerror("pei_tseng_notch([1, 2, 3])", "Wrong number of input arguments.");
+//assert_checkerror("pei_tseng_notch([1, 2; 3, 4], [4; 5; 6])", "All arguments must be vectors or scalars.");
+//assert_checkerror("pei_tseng_notch([1, 2, 3, 4], [5, 6])", "All arguments must be of equal length.");
+//assert_checkerror("pei_tseng_notch([1, 2, 3], [4, 5, 6])", "All frequencies must be in (0, 1).");
+//assert_checkerror("pei_tseng_notch([0.1, 0.2, 0.3], [4, 5, 6])", "All bandwidths must be in (0, 1).");
+
+//tests:
+//assert_checkequal(pei_tseng_notch(0.2, 0.4), [0 0 0]);
+//assert_checkalmostequal(pei_tseng_notch([0.2, 0.5, 0.9], [0.1, 0.3, 0.5]), [0.41549, 0.11803, -0.03227, 0.23606, -0.03227, 0.11803, 0.41549], 5*10^-4);
+//assert_checkequal(pei_tseng_notch([0.2; 0.5; 0.9], [0.1, 0.3, 0.5]), pei_tseng_notch([0.2, 0.5, 0.9], [0.1; 0.3; 0.5]));
+//assert_checkalmostequal(pei_tseng_notch([0.2, 0.7], [0.1, 0.3]), [0.60331, -0.26694, 0.05905, -0.26694, 0.60331], 5*10^-4);
diff --git a/macros/polystab.sci b/macros/polystab.sci
index 84e42d4..c9da6f0 100644
--- a/macros/polystab.sci
+++ b/macros/polystab.sci
@@ -1,21 +1,38 @@
-function b = polystab(a)
-//This function stabilizes the polynomial transfer function.
-//Calling Sequence
-//b = polystab(a)
-//Parameters
-//a:
-//Description
-//This is an Octave function.
-//This function stabilizes the polynomial transfer function by replacing all roots outside the unit circle with their reflection inside the unit circle.
-//Examples
+function y = polystab(x)
+//Stabilize the polynomial transfer function.
+//Calling Sequence:
+//y = polystab(x)
+//Parameters:
+//x: real or complex valued vector
+//Description:
+//This function stabilizes the polynomial transfer function by replacing all
+//roots outside the unit circle with their reflection inside the unit circle.
+//Example:
//polystab([1,3,5])
-//ans =
-// 1. 0.6 0.2
+//ans=
+//[1, 0.6, 0.2]
+
+ funcprot(0);
+ if (argn(2) ~= 1)
+ error("polystab: wrong number of input arguments");
+ end
+ root = roots(x);
+ vals = find(abs(root) > 1);
+ root(vals) = 1 ./ conj(root(vals));
+ b = x(1) * poly(root, 'x');
+ y = flipdim(coeff(b), 2);
+ if (isreal(x))
+ y = real(y);
+ end
-funcprot(0);
-rhs = argn(2)
-if(rhs~=1)
-error("Wrong number of input arguments.")
-end
-b = callOctave("polystab",a)
endfunction
+
+//input validation:
+//assert_checkerror("polystab()", "polystab: wrong number of input arguments");
+//assert_checkerror("polystab(1, 2)", "Wrong number of input arguments.")
+
+//tests:
+//assert_checkalmostequal(polystab([1, 3, 5]), [1, 0.6, 0.2], %eps);
+//assert_checkequal(polystab([1; 3; 5]), polystab([1, 3, 5]));
+//assert_checkalmostequal(polystab([1+3*%i, 4+2*%i, 2+2*%i]), [1 + 3*%i, 1.35664 + 1.84888*%i, 0.88566 + 0.88566*%i], 5*10^-5);
+//assert_checkalmostequal(polystab([-1; 3+%i; -2-6*%i]), [-1, 0.3 + 0.4*%i, -0.05-0.15*%i], 5*10^-5);
diff --git a/macros/schtrig.sci b/macros/schtrig.sci
index 90286e2..f50adca 100644
--- a/macros/schtrig.sci
+++ b/macros/schtrig.sci
@@ -1,30 +1,127 @@
-function v = schtrig (x, lev, rs)
-//This function implements a multisignal Schmitt triggers with lev levels supplied as input.
-//Calling Sequence
-//v = schtrig (x, lev)
-//v = schtrig (x, lev, rs)
-//Parameters
-//x: vector or matrix of real numbers
-//lev: real number
-//rs: default value 1
-//Description
-//This is an Octave function.
-//This function implements a multisignal Schmitt triggers with lev levels supplied as input.
-//The argument 1 is a matrix (or a vector) and this trigger works along its first dimension.
-//Examples
-//schtrig([0.2,-3,5],-4)
-//ans =
-// 0. 0. 1.
-
-funcprot(0);
-rhs = argn(2)
-if(rhs<2 | rhs>3)
-error("Wrong number of input arguments.")
-end
-if(rhs==2)
-v = callOctave("schtrig", x, lev)
-elseif(rhs==3)
-v = callOctave("schtrig",x, lev, rs)
-end
+function [v, rg] = schtrig (x, lvl, rst)
+//This function implements a multisignal Schmitt triggers with 'lvl' levels supplied as input.
+//Calling Sequence:
+//x = schtrig(x, lvl)
+//lvl = schtrig(x, lvl, rst)
+//Parameters:
+//x: Vector or matrix of real numbers
+//lvl: Real number
+//rst: Boolean, default value is 'true'
+//Description:
+//This function implements a multisignal Schmitt triggers with 'lvl' levels supplied as input.
+//The argument 1 is a matrix (or a vector) and this trigger works along its first dimension.
+//Examples:
+//schtrig([0.2,-3,5], -4)
+//ans = [0, 0, 1]
+
+ funcprot(0);
+ warning('off');
+ if (argn(2) < 2 | argn(2) > 3)
+ error("Wrong number of input arguments.");
+ elseif (argn(2) == 2)
+ rst = %T;
+ end
+ if (length (ndims (x)) > 2)
+ error ('The input should be two dimensional.');
+ end
+ if (length (ndims (lvl)) > 2)
+ error ('Only a maximum of two threshold levels accepted.');
+ end
+
+ [nT nc] = size (x);
+
+ global st0;
+ if (rst || isempty (st0))
+ st0 = zeros (1,nc);
+ end
+
+ if (length(lvl) == 1)
+ lvl = abs (lvl) .* [1 -1];
+ else
+ lvl = gsort(lvl, 'g', 'd');
+ end
+
+ v = repmat(%nan, nT, nc);
+ v(1,:) = st0;
+
+ up = x > lvl(1);
+ v(up) = 1;
+
+ dw = x < lvl(2);
+ v(dw) = 0;
+
+ idx = bool2s(isnan(v));
+ xhi = idx';
+ bool_discon = diff (xhi, 1, 2);
+ [Np Na] = size (xhi);
+ ranges = cell (1, Np);
+
+ for i = 1:Np
+ idxUp = find (bool_discon(i,:) > 0) + 1;
+ idxDwn = find (bool_discon(i,:) < 0);
+ tLen = length (idxUp) + length (idxDwn);
+
+ if (xhi(i,1) == 1)
+ ranges{i}(1) = 1;
+ ranges{i}(2:2:tLen+1) = idxDwn;
+ ranges{i}(3:2:tLen+1) = idxUp;
+ else
+ ranges{i}(1:2:tLen) = idxUp;
+ ranges{i}(2:2:tLen) = idxDwn;
+ end
+
+ if (xhi(i, $) == 1)
+ ranges{i}($+1) = Na;
+ end
+
+ tLen = length(ranges{i});
+ if (tLen ~= 0)
+ ranges{i} = matrix(ranges{i}, 2, tLen / 2);
+ end
+
+ end
+
+ if (Np == 1)
+ ranges = cell2mat(ranges);
+ end
+ if (nc == 1)
+ ranges = {ranges};
+ end
+
+ for i=1:nc
+ if (~isempty(ranges{i}))
+ prev = ranges{i}(1,:)-1;
+ prev(prev<1) = 1;
+ st0 = v(prev, i);
+ ini_idx = ranges{i}(1,:);
+ end_idx = ranges{i}(2,:);
+ for j =1:length(ini_idx)
+ v(ini_idx(j):end_idx(j),i) = st0(j);
+ end
+ end
+ end
+
+ st0 = v($, :);
endfunction
+
+//input validation:
+//assert_checkerror("schtrig(1)", "Wrong number of input arguments.");
+//assert_checkerror("schtrig(1, 2, 3, 4)", "Wrong number of input arguments.");
+
+//tests:
+//assert_checkequal(schtrig(ones(128, 1), -1), zeros(128, 1));
+//assert_checkequal(schtrig(ones(128, 1), -1), schtrig(ones(128, 1), -1, %F));
+//
+//assert_checkequal(schtrig([1 5 1], 3), schtrig([1 5 1], 3, %T));
+//assert_checkequal(schtrig([1 5 1], 3), [0 1 0]);
+//assert_checkequal(schtrig([1 5 1], 3, %F), schtrig([1 5 1], 3));
+//
+//assert_checkequal(schtrig([-3; 4; -1], -2), [0; 1; 1]);
+//assert_checkequal(schtrig([-3; 4; -1], -2, %F), schtrig([-3; 4; -1], -2));
+//
+//assert_checkequal(schtrig([1 -3; 2 4; 5 -1], -2), [0 0; 0 1; 1 1]);
+//assert_checkequal(schtrig([1 -3; 2 4; 5 -1], -2, %F), [1 0; 1 1; 1 1]);
+//
+//assert_checkequal(schtrig([1.1 1.5 0; 0.53 1.21 3.57; 5.34 -1.24 0], 0.424), [1 1 0;1 1 1;1 0 1]);
+//assert_checkequal(schtrig([1.1 1.5 0; 0.53 -1.21 3.57; 5.34 -1.24 0], 0.424, %F), [1 1 1;1 0 1;1 0 1]);
diff --git a/macros/sinetone.sci b/macros/sinetone.sci
index 407fa45..f174869 100644
--- a/macros/sinetone.sci
+++ b/macros/sinetone.sci
@@ -1,32 +1,66 @@
-function y= sinetone(x, varargin)
+function y = sinetone(freq, rate, sec, ampl)
//Return a sinetone of the input
-//Calling Sequence
-//y= sinetone(FREQ)
-//y= sinetone(FREQ, RATE)
-//y= sinetone(FREQ, RATE, SEC)
-//y= sinetone(FREQ, RATE, SEC, AMPL)
-//Parameters
-//FREQ: frequency of sinetone
-//RATE: Sampling rate
-//SEC: Length in seconds
-//AMPL: Amplitude
-//Description
-//Return a sinetone of frequency FREQ with a length of SEC seconds atsampling rate RATE and with amplitude AMPL.The arguments FREQ and AMPL may be vectors of common size.The defaults are RATE = 8000, SEC = 1, and AMPL = 64.
-funcprot(0);
-rhs= argn(2);
-if(rhs<1 | rhs>4)
-error("Wrong number of input arguments")
-end
-
-select(rhs)
- case 1 then
- y= callOctave("sinetone", x);
- case 2 then
- y= callOctave("sinetone", x , varargin(1));
- case 3 then
- y= callOctave("sinetone", x , varargin(1), varargin(2));
- case 4 then
- y= callOctave("sinetone", x , varargin(1), varargin(2), varargin(3));
-
-end
+//Calling Sequence:
+//y= sinetone(freq)
+//y= sinetone(freq, rate)
+//y= sinetone(freq, rate, sec)
+//y= sinetone(freq, rate, sec, ampl)
+//Parameters:
+//freq: frequency of sinetone
+//rate: sampling rate
+//sec: length in seconds
+//ampl: amplitude
+//Description:
+//Return a sinetone of frequency 'freq' with a length of 'sec' seconds at sampling rate 'rate' and with amplitude 'ampl'.
+//The arguments freq and ampl may be vectors of common size. The defaults are rate = 8000, sec = 1, and ampl = 64.
+//Examples:
+//sinetone(5, 2, 1, 8)
+//ans = [4.89E-15; -9.79E-15]
+
+ funcprot(0);
+ rhs= argn(2);
+
+ select rhs
+ case 1 then
+ rate = 8000;
+ sec = 1;
+ ampl = 64;
+ case 2 then
+ sec = 1;
+ ampl = 64;
+ case 3 then
+ ampl = 64;
+ case 4 then
+ break;
+ else
+ error("sinetone: wrong number of input arguments");
+ end
+
+ if ( isvector(freq) & isvector(ampl) & ~isequal(size(freq), size(ampl)) ) then
+ error("sinetone: freq and ampl must be vectors of equal size");
+ end
+
+ if ( ~(isscalar(rate) && isscalar(sec)) ) then
+ error("sinetone: rate and sec must be scalars");
+ end
+
+ n = length (freq);
+ ns = round (rate * sec);
+ y = zeros (ns, n);
+
+ for k = 1:n
+ y(:, k) = ampl(k) * sin (2 * %pi * (1:ns) / rate * freq(k))';
+ end
+
endfunction
+
+//input validation:
+//assert_checkerror("sinetone()", "sinetone: wrong number of input arguments");
+//assert_checkerror("sinetone([6, 9, 4, 2], 2, 3, [6, 2, 0])", "sinetone: freq and ampl must be vectors of equal size");
+//assert_checkerror("sinetone(1, [1, 2])", "sinetone: rate and sec must be scalars");
+//assert_checkerror("sinetone(1, 2, [4, 3])", "sinetone: rate and sec must be scalars");
+
+//tests:
+//assert_checkequal(size(sinetone(18e6, 150e6, 19550/150e6, 1)), [19550, 1]);
+//assert_checkequal(sinetone(5), sinetone(5, 8000, 1, 64));
+//assert_checkequal(size(sinetone([1, 2, 3], 4000, 1, [8, 9, 6])), [4000, 3]);
diff --git a/macros/sinewave.sci b/macros/sinewave.sci
index 5263f35..31ee6b0 100644
--- a/macros/sinewave.sci
+++ b/macros/sinewave.sci
@@ -1,29 +1,57 @@
-function y= sinewave(x, varargin)
-//Return an M-element vector with I-th element given by 'sin(2* pi *(I+D-1)/N).'
-//Calling Sequence
-//y= sinewave(M)
-//y= sinewave(M,N)
-//y= sinewave(M,N,D)
-//Parameters
-//M: Input vector
-//N: The default value for N is M
-//D: The default value for D is 0
-//AMPL: Amplitude
-//Description
-//Return an M-element vector with I-th element given by 'sin(2* pi *(I+D-1)/N).'
-funcprot(0);
-rhs= argn(2);
-if(rhs<1 | rhs>3)
-error("Wrong number of input arguments")
-end
-
-select(rhs)
- case 1 then
- y= callOctave("sinewave", x);
- case 2 then
- y= callOctave("sinewave", x , varargin(1));
- case 3 then
- y= callOctave("sinewave", x , varargin(1), varargin(2));
-
-end
+function y = sinewave (m, n, d)
+//Return a m-element vector with the i-th element given by {sin (2 * pi * (i+d-1) / n}.
+//Calling Sequence:
+//y= sinewave(m)
+//y= sinewave(m,n)
+//y= sinewave(m,n,d)
+//Parameters:
+//m: Real positive scalar
+//n: The value of n in the formula {sin (2 * pi * (i+d-1) / n}. The default value for n is m.
+//d: The value of d in the formula {sin (2 * pi * (i+d-1) / n}. The default value for d is 0.
+//Examples:
+//sinewave(1, 4, 1)
+//ans = 1
+
+ funcprot(0);
+ rhs = argn(2);
+ if (rhs < 1 | rhs > 3) then
+ error("sinewave: wrong number of input arguments");
+ end
+
+ if (m < 0) then
+ error("sinewave: arg1 (m) must be non-negative");
+ end
+
+ if (rhs < 2) then
+ n = m;
+ end
+ if (n == 0) then
+ error("sinewave: arg2 (n) must be non-zero");
+ end
+
+ if (~isscalar(m)) then
+ error("sinewave: arg1 (m) must be a real scalar");
+ end
+
+ if (rhs < 3) then
+ d = 0;
+ end
+
+ y = sin(((1 : m) + d - 1) * 2 * %pi / n);
+
endfunction
+
+//input validation:
+//assert_checkerror("sinewave()", "sinewave: wrong number of input arguments");
+//assert_checkerror("sinewave(-1, 2)", "sinewave: arg1 (m) must be non-negative");
+//assert_checkerror("sinewave([-2, 5, -3], 2)", "sinewave: arg1 (m) must be a real scalar");
+//assert_checkerror("sinewave(0, 0, 5)", "sinewave: arg2 (n) must be non-zero");
+
+//tests:
+//assert_checkequal(sinewave(1), 0);
+//assert_checkequal(sinewave(1), sinewave(1, 1, 0));
+//assert_checkequal(sinewave(3, 4), sinewave(3, 4, 0));
+//assert_checkalmostequal(sinewave(1, 12, 1), 1/2, %eps);
+//assert_checkalmostequal(sinewave(1, 20, 1), (sqrt(5)-1)/4, %eps);
+//assert_checkequal(sinewave(3, 4, 1), [1, sin(%pi), -1]);
+//assert_checkalmostequal(sinewave(4, 12, 1), [1/2, sqrt(3)/2, 1, sqrt(3)/2], %eps);
diff --git a/macros/spencer.sci b/macros/spencer.sci
index 63a1b83..1b0f377 100644
--- a/macros/spencer.sci
+++ b/macros/spencer.sci
@@ -1,21 +1,37 @@
-function y= spencer(x)
-//Return Spencer's 15 point moving average of each column of X.
-//Calling Sequence
-//spencer(X)
-//Parameters
-//X: Real scalar or vector
-//Description
-//Return Spencer's 15 point moving average of each column of X.
-funcprot(0);
+function savg = spencer (x)
+//Returns Spencer's 15 point moving average of each column of x.
+//Calling Sequence:
+//spencer(x)
+//Parameters:
+//X: Real vector or matrix
+//Description:
+//Returns Spencer's 15 point moving average of each column of x.
-rhs= argn(2);
+ funcprot(0);
+ if (nargin() ~= 1)
+ error("Wrong number of input arguments.");
+ end
+ [xr, xc] = size (x);
-if(rhs <1 | rhs >1)
-error("Wrong number of input arguments");
-end
+ n = xr;
+ c = xc;
+
+ if (isvector (x))
+ n = length (x);
+ c = 1;
+ x = matrix(x, n, 1);
+ end
+
+ end
+ w = [-3, -6, -5, 3, 21, 46, 67, 74, 67, 46, 21, 3, -5, -6, -3] / 320;
+ savg = filter (w, 1, x);
+ savg = [zeros(7,c); savg(15:n,:); zeros(7,c);];
+ savg = matrix(savg, xr, xc);
-select(rhs)
- case 1 then
- y = callOctave("spencer",x);
-end
endfunction
+
+//tests:
+//assert_checkerror("spencer()", "Wrong number of input arguments.");
+//assert_checkerror("spencer(1, 2)", "Wrong number of input arguments.");
+//assert_checkequal(spencer(linspace(1, 14, 14)'), zeros(14, 1));
+//assert_checkequal(spencer(linspace(-1, -10, 14)), zeros(1, 14));
diff --git a/macros/stft.sci b/macros/stft.sci
index 7c30360..1d4b462 100644
--- a/macros/stft.sci
+++ b/macros/stft.sci
@@ -1,89 +1,139 @@
-function [y,c]= stft(x, varargin)
-//Compute the short-time Fourier transform of the vector X
-//Calling Sequence
-//Y = stft (X)
-//Y = stft (X, WIN_SIZE)
-//Y = stft (X, WIN_SIZE, INC)
-//Y = stft (X, WIN_SIZE, INC, NUM_COEF)
-//Y = stft (X, WIN_SIZE, INC, NUM_COEF, WIN_TYPE)
-//[Y,C] = stft (X)
-//[Y,C] = stft (X, WIN_SIZE)
-//[Y,C] = stft (X, WIN_SIZE, INC)
-//[Y,C] = stft (X, WIN_SIZE, INC, NUM_COEF)
-//[Y,C] = stft (X, WIN_SIZE, INC, NUM_COEF, WIN_TYPE)
-//Parameters
-//X: Real scalar or vector
-//WIN_SIZE: Size of the window used
-//INC: Increment
-//WIN_TYPE: Type of window
-//Description
-//Compute the short-time Fourier transform of the vector X with NUM_COEF coefficients by applying a window of WIN_SIZE data points and an increment of INC points.
-//
+function [y, c] = stft(x, win_size, inc, num_coef, win_type)
+//Compute the short-time Fourier transform of the vector X.
+//Calling Sequence:
+//y = stft(x)
+//y = stft(x, win_size)
+//y = stft(x, win_size, inc)
+//y = stft(x, win_size, inc, num_coef)
+//y = stft(x, win_size, inc, num_coef, win_type)
+//[y, c] = stft(x)
+//[y, c] = stft(x, win_size)
+//[y, c] = stft(x, win_size, inc)
+//[y, c] = stft(x, win_size, inc, num_coef)
+//[y, c] = stft(x, win_size, inc, num_coef, win_type)
+//Parameters:
+//x: Real scalar or vector
+//win_size: Size of the window used
+//inc: Increment
+//num_coef: Coefficients of Fourier transform
+//win_type: Type of window
+//Description:
+//Compute the short-time Fourier transform of the vector X with num_coef coefficients by applying a window of win_size data points and an increment of inc points.
//Before computing the Fourier transform, one of the following windows is applied:
-//
//"hanning" -> win_type = 1
-//
//"hamming" -> win_type = 2
-//
//"rectangle" -> win_type = 3
-//
-//The window names can be passed as strings or by the WIN_TYPE number.
-//
-//The following defaults are used for unspecified arguments:WIN_SIZE= 80, INC = 24, NUM_COEF = 64, and WIN_TYPE = 1.
-//
-//Y = stft (X, ...)' returns the absolute values of the Fourier coefficients according to the NUM_COEF positive frequencies.
-//
-//'[Y, C] = stft (x, ...)' returns the entire STFT-matrix Y and a 3-element vector C containing the window size, increment, and window type, which is needed by the 'synthesis' function.
+//The window names can be passed as strings or by the win_type number.
+//The following defaults are used for unspecified arguments: win_sizse = 80, inc = 24, num_coef = 64, and win_type = 1.
+//y = stft(x, ...)' returns the absolute values of the Fourier coefficients according to the NUM_COEF positive frequencies.
+//'[y, c] = stft (x, ...)' returns the entire stft-matrix y and a 3-element vector c containing the window size, increment, and window type,
+//which is needed by the 'synthesis' function.
+//Examples:
+//[y, c] = stft([1; -2; -5])
+//y = []
+//c = [80, 24, 1]
+
+ funcprot(0);
+ rhs = argn(2);
+ if (rhs < 1 | rhs > 5)
+ error("Wrong number of input arguments.");
+ end
+ if (~isvector (x))
+ error ("stft: X must be a vector");
+ end
+
+ if (rhs == 1)
+ win_size = 80;
+ inc = 24;
+ num_coef = 64;
+ win_type = 1;
+ elseif (rhs == 2)
+ inc = 24;
+ num_coef = 64;
+ win_type = 1;
+ elseif (rhs == 3)
+ num_coef = 64;
+ win_type = 1;
+ elseif (rhs == 4)
+ win_type = 1;
+ end
+
+ if (type(win_type) == 10)
+ switch (convstr(win_type, 'l'))
+ case "hanning" , win_type = 1;
+ case "hamming" , win_type = 2;
+ case "rectangle" , win_type = 3;
+ otherwise
+ error ("stft: unknown window type");
+ end
+ end
-funcprot(0);
-lhs= argn(1);
-rhs= argn(2);
+ x = x(:);
-if(rhs <1 | rhs>5)
- error("Wrong number of input arguments");
-end
+ ncoef = 2 * num_coef;
+ if (win_size > ncoef)
+ win_size = ncoef;
+ printf ("stft: window size adjusted to %f\n", win_size);
+ end
+ num_win = fix ((size(x, 'r') - win_size) / inc);
-if(lhs<1 | lhs>2)
- error("Wrong number of output arguments");
-end
+ switch (win_type)
+ case 1 , win_coef = (window('hn', win_size))';
+ case 2 , win_coef = (window('hm', win_size))';
+ case 3 , win_coef = ones (win_size, 1);
+ end
+
+ z = zeros (ncoef, num_win + 1);
+ start = 1;
+ for i = 0:num_win
+ z(1:win_size, i+1) = x(start:start+win_size-1) .* win_coef;
+ start = start + inc;
+ end
+
+ y = fft(z, -1, find(size(z) ~= 1, 1));
+
+ if (nargout() == 1)
+ y = abs(y(1:num_coef, :));
+ else
+ c = [win_size, inc, win_type];
+ end
-select(rhs)
- case 1 then
- select(lhs)
- case 1 then
- y= callOctave("stft", x);
- case 2 then
- [y,c]= callOctave("stft", x);
- end
- case 2 then
- select(lhs)
- case 1 then
- y= callOctave("stft", x,varargin(1));
- case 2 then
- [y,c]= callOctave("stft", x, varargin(1));
- end
- case 3 then
- select(lhs)
- case 1 then
- y= callOctave("stft", x,varargin(1), varargin(2));
- case 2 then
- [y,c]= callOctave("stft", x,varargin(1), varargin(2));
- end
- case 4 then
- select(lhs)
- case 1 then
- y= callOctave("stft", x,varargin(1), varargin(2), varargin(3));
- case 2 then
- [y,c]= callOctave("stft", x,varargin(1), varargin(2), varargin(3));
- end
- case 5 then
- select(lhs)
- case 1 then
- y= callOctave("stft", x,varargin(1), varargin(2), varargin(3), varargin(4));
- case 2 then
- [y,c]= callOctave("stft", x,varargin(1), varargin(2), varargin(3), varargin(4));
- end
-end
endfunction
+//input validation:
+//assert_checkerror("stft()", "Wrong number of input arguments.");
+//assert_checkerror("stft(1, 2, 3, 4, 5, 6)", "Wrong number of input arguments.");
+//assert_checkerror("stft([1 2; 3 4])", "stft: X must be a vector");
+//s = "square";
+//assert_checkerror("stft([1 2 3], 8, 4, 6, s)", "stft: unknown window type");
+
+//tests:
+//assert_checkequal(stft([1 2 3]), stft([1 21 3], 80, 24, 64, 1));
+//assert_checkequal(stft([1 2 3], 80), stft([1 21 3], 80, 24, 64, 1));
+//assert_checkequal(stft([1 2 3], 80, 24), stft([1 21 3], 80, 24, 64, 1));
+//assert_checkequal(stft([1 2 3], 80, 24, 64), stft([1 21 3], 80, 24, 64, 1));
+
+//[y, c] = stft([1; -2; -5]);
+//assert_checkequal(y, []);
+//assert_checkequal(c, [80, 24, 1]);
+
+//y = stft([1 21 3], 3, 2, 2);
+//assert_checkequal(y, [21; 21]);
+
+//[y, c] = stft([1 21 3], 3, 2, 2);
+//assert_checkequal(y, [21; -21*%i; -21; 21*%i]);
+//assert_checkequal(c, [3 2 1]);
+
+//y = stft([1, 3-2*%i, -6-7*%i], 3, 2, 3);
+//assert_checkalmostequal(y, [3.60555; 3.60555; 3.60555], 0.5*10^-5);
+
+//[y c] = stft([1; 3-2*%i; -6-7*%i], 3, 2, 3);
+//assert_checkalmostequal(y, [3-2*%i; -0.2321-3.5981*%i; -3.2321-1.5981*%i; -3+2*%i; 0.2321+3.5981*%i; 3.2321+1.5981*%i], 5*10^-4);
+//assert_checkequal(c, [3 2 1]);
+
+//[y c] = stft([1; 3-2*%i; -1-2*%i], 3, 2, 3, 3);
+//assert_checkalmostequal(y, [3-4*%i; -0.4641-1.7321*%i; 0-1.4641*%i; -3; 5.4641*%i; 6.4641+1.7321*%i], 5*10^-4);
+//assert_checkequal(c, [3 2 3]);
+//y = stft([3*%i; 4*%i; -5*%i], 3, 2, 3, 1);
+//assert_checkequal(y, [4; 4; 4]);
diff --git a/macros/synthesis.sci b/macros/synthesis.sci
index 7224686..ede8601 100644
--- a/macros/synthesis.sci
+++ b/macros/synthesis.sci
@@ -1,25 +1,67 @@
-function x= synthesis(Y,C)
-//Compute a signal from its short-time Fourier transform
-//Calling Sequence
-//X= synthesis(Y,C)
-//Parameters
-//Y: Shirt-time fourier transform
-//C: 3-element vector C specifying window size, increment, window type.
-//Description
-//Compute a signal from its short-time Fourier transform Y and a 3-element vector C specifying window size, increment, and window type.
-//The values Y and C can be derived by
-//[Y, C] = stft (X , ...)
-funcprot(0);
-lhs= argn(1);
-rhs= argn(2);
-
-if(rhs<2 | rhs >2)
- error("Wrong number of input arguments");
-end
-
-select(rhs)
- case 2 then
- x= callOctave("synthesis", Y,C);
-
-end
+function x = synthesis (y, c)
+//Compute a signal from its short-time Fourier transform.
+//Calling Sequence:
+//synthesis(y, c)
+//Parameters:
+//y: Real or complex matrix representing a signal's short-time fourier transform.
+//c: 3-element vector C specifying window size, increment and window type.
+//Description:
+//Compute a signal from its short-time Fourier transform 'y' and a 3-element vector 'c' specifying window size, increment, and window type.
+//A window type of 1 represents a hanning window, 2 represents a hamming window and 3 represents a rectangular window.
+//The values 'y' and 'c' can be derived by [y, c] = stft (x, ...)
+
+ funcprot(0);
+ rhs = argn(2);
+ if (rhs ~= 2)
+ error("synthesis: wrong number of input arguments");
+ end
+
+ if (length(c) ~= 3)
+ error ("synthesis: C must contain exactly 3 elements");
+ end
+
+ w_size = c(1);
+ inc = c(2);
+ w_type = c(3);
+
+ if (w_type == 1)
+ w_coeff = window('hn', w_size);
+ elseif (w_type == 2)
+ w_coeff = window('hm', w_size);
+ elseif (w_type == 3)
+ w_coeff = ones (w_size, 1);
+ else
+ error ("synthesis: window_type must be 1, 2, or 3");
+ end
+
+ if (isnan(w_coeff))
+ w_coeff = 1;
+ end
+ z = real(fft(y, 1, find(size(y) ~= 1, 1)));
+ st = fix((w_size-inc) / 2);
+ z = z(st+1:st+inc, :);
+ w_coeff = w_coeff(st+1:st+inc);
+
+ nc = size(z, 'c');
+ for i = 1:nc
+ z(:, i) = z(:, i) ./ w_coeff;
+ end
+
+ x = matrix(z, inc * nc, 1);
+
endfunction
+
+//input validation:
+//assert_checkerror("synthesis(1)", "synthesis: wrong number of input arguments");
+//assert_checkerror("synthesis(1, [1 2])", "synthesis: C must contain exactly 3 elements");
+//assert_checkerror("synthesis(1, [1 2 3 4])", "synthesis: C must contain exactly 3 elements");
+//assert_checkerror("synthesis(1, [1 2 4])", "synthesis: window_type must be 1, 2, or 3");
+
+//tests:
+//assert_checkequal(synthesis(1, [1 1 1]), 1);
+//assert_checkequal(synthesis(4, [2 1 3]), 4);
+//assert_checkalmostequal(synthesis(-6, [2 1 2]), -75, %eps);
+//assert_checkalmostequal(synthesis(-9, [2; 1; 2]), -112.50, 10*%eps);
+//assert_checkequal(synthesis(5*%i, [2; 1; 3]), 0);
+//assert_checkequal(synthesis([1 2; -3, -6], [1; 1; 2]), [-1; -2]);
+//assert_checkalmostequal(synthesis([1 1+2*%i; -3-4*%i, -5], [2; 1; 2]), [-12.5; -25], %eps);
diff --git a/macros/tf2sos.sci b/macros/tf2sos.sci
index 7570707..ad7f10d 100644
--- a/macros/tf2sos.sci
+++ b/macros/tf2sos.sci
@@ -1,35 +1,70 @@
-function [sos,varargout] = tf2sos (b, a)
+function [sos,g] = tf2sos (B, A)
//This function converts direct-form filter coefficients to series second-order sections.
-//Calling Sequence
+//Calling Sequence:
//[sos] = tf2sos (b, a)
//[sos, g] = tf2sos (b, a)
-//Parameters
+//Parameters:
//b: matrix of real numbers
//a: matrix of real numbers
-//Description
-//This is an Octave function.
+//Description:
//This function converts direct-form filter coefficients to series second-order sections.
//The input parameters b and a are vectors specifying the digital filter H(z) = B(z)/A(z).
//The output is the sos matrix and the overall gain.
//If there is only one output argument, the overall filter gain is applied to the first second-order section in the sos matrix.
-//Examples
+//Examples:
//tf2sos([1,2,3,4,5,6],2)
//ans =
-// 0.50000 0.80579 1.07239 1.00000 0.00000 0.00000
-// 1.00000 -1.10337 1.87524 1.00000 0.00000 0.00000
-// 1.00000 1.49180 -0.00000 1.00000 0.00000 0.00000
-
-funcprot(0);
-rhs = argn(2)
-lhs = argn(1)
-
-if(rhs~=2)
-error("Wrong number of input arguments.")
-end
- select(lhs)
- case 1 then
- sos = callOctave("tf2sos",b,a)
- case 2 then
- [sos,g] = callOctave("tf2sos",b,a)
- end
-endfunction
+// 0.50000 0.80579 1.07239 0.00000 0.00000 1.00000
+// 1.00000 -1.10337 1.87524 0.00000 0.00000 1.00000
+// 1.00000 1.49180 -0.00000 0.00000 1.00000 0.00000
+
+ funcprot(0);
+ if (nargin() ~= 2) then
+ error("Wrong number of input arguments.");
+ end
+ S = syslin([], inv_coeff(flipdim(B(:)', 2)), inv_coeff(flipdim(A(:)', 2)));
+ [z,p,k] = tf2zp(S);
+ if (nargout() < 2)
+ sos = clean(zp2sos(z,p,k));
+ else
+ [sos,g] = clean(zp2sos(z,p,k));
+ end
+
+endfunction
+
+//tests:
+//assert_checkerror("tf2sos(1)", "Wrong number of input arguments.");
+//assert_checkerror("tf2sos(1, 2, 3)", "Wrong number of input arguments.");
+
+//A = [-7 -12];
+//B = [-1 -7];
+//assert_checkequal(tf2sos(A, B), [7 12 0 1 7 0]);
+
+//A = [1 2; 3 4];
+//B = [3 2; 1 4];
+//assert_checkalmostequal(tf2sos(A, B), [0.3333 0.0679 0.4768 1 -0.6667 1.3333; 1 2.7963 0 1 1 0], 5*10^-4);
+
+//A = [1 1 1 1];
+//B = [-1 -1; -1 -1];
+//assert_checkalmostequal(tf2sos(A, B), [-1 0 -1 1 0 1; 1 1 0 1 1 0], 100*%eps);
+
+//b = [1, 1];
+//a = [1, 2];
+//expected_sos = [1, 1, 0, 1, 2, 0];
+//sos = tf2sos(b, a);
+//assert_checkalmostequal(sos, expected_sos, %eps);
+
+//b = [1, 2, 2];
+//a = [1, -1, 1];
+//expected_sos = [1, 2, 2, 1, -1, 1];
+//sos = tf2sos(b, a);
+//assert_checkalmostequal(sos, expected_sos, 100*%eps);
+
+//b = [1, 1, 1, 1];
+//a = [1, 1, 1, 1];
+//expected_sos = [
+// 1, 0, 1, 1, 0, 1;
+// 1, 1, 0, 1, 1, 0
+//];
+//sos = tf2sos(b, a);
+//assert_checkalmostequal(sos, expected_sos, 100*%eps);
diff --git a/macros/tfestimate.sci b/macros/tfestimate.sci
index e52faf5..f437f52 100644
--- a/macros/tfestimate.sci
+++ b/macros/tfestimate.sci
@@ -1,93 +1,73 @@
-function [Pxx, freq] = tfestimate(x, y, window, overlap, Nfft, Fs, range)
-
+function varargout = tfestimate(varargin)
//Estimate transfer function of system with input x and output y. Use the Welch (1967) periodogram/FFT method.
-//Calling Sequence
+//This function depends on pwelch (and hence, fft1).
+//Calling Sequence:
+//tfestimate (x, y)
+//tfestimate (x, y, window)
+//tfestimate (x, y, window, overlap)
+//tfestimate (x, y, window, overlap, Nfft)
+//tfestimate (x, y, window, overlap, Nfft, Fs)
//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.
+//Parameters:
+//x: system-input time-series data. Real pr complex vector.
+//y: system-output time-series data. Real or complex vector, same dimension as x.
+//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.
+// 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. '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. If data (x and y) are real, the default range is 'half', otherwise default range is 'whole'.
+//Description:
//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])
+//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
+ funcprot(0);
+ if (nargin() < 2 || nargin() > 7)
+ error("Wrong number of input arguments.");
+ end
+ nvarargin = length(varargin);
+ for iarg = 1:nvarargin
+ arg = varargin(iarg);
+ if ( ~isempty(arg) && (type(arg) == 10) && ( ~strcmp(arg,'power') || ~strcmp(arg,'cross') || ~strcmp(arg,'trans') || ~strcmp(arg,'coher') || ~strcmp(arg,'ypower') ))
+ varargin(iarg) = [];
+ end
+ end
+ varargin(nvarargin + 1) = 'trans';
-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
+ if (nargout() == 0)
+ pwelch(varargin(:));
+ elseif (nargout() == 1)
+ Pxx = pwelch(varargin(:));
+ varargout(1) = Pxx;
+ elseif (nargout() >= 2)
+ [Pxx,f] = pwelch(varargin(:));
+ varargout(1) = Pxx;
+ varargout(2) = f;
+ 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
-
+//tests:
+//assert_checkerror("tfestimate(1)", "Wrong number of input arguments.");
+//assert_checkerror("tfestimate(1, 2, 3, 4, 5, 6, 7, 8)", "Wrong number of input arguments.");
+//tfestimate([1 2 3], [4 5 6]);
+//tfestimate([-1 -2 -3], [4 5 6], 'power');
+//tfestimate([-1; -2; -3], [-4 -5 -6], 'power', 'cross');
+//y = tfestimate([1; 2; 3], [4; 5; 6], 'coher', 'trans');
+//assert_checkalmostequal(y, [1.75; 1.59472+0.38256*%i; 1.28235], 5*10^-5);
+//[a b] = tfestimate([1 -2 -3], [4; 5; 6], 'ypower', 'trans');
+//assert_checkalmostequal(a, [-1.1667; -1.1029-0.1912*%i; -1.5797], 5*10^-4);
+//assert_checkalmostequal(b, [0; 0.25; 0.5], %eps);
+//tfestimate([5-1 4+2*%i 4 -1-3*%i], [1+%i -5*%i 6 -6]);
diff --git a/macros/triang.sci b/macros/triang.sci
index f450d97..596bfab 100644
--- a/macros/triang.sci
+++ b/macros/triang.sci
@@ -6,7 +6,6 @@ function w = triang (m)
//m: positive integer value
//w: output variable, vector of real numbers
//Description
-//This is an Octave function.
//This function returns the filter coefficients of a triangular window of length m supplied as input, to the output vector y.
//Examples
//triang(5)
@@ -29,3 +28,16 @@ rhs = argn(2)
w = 1 - abs ([-(m-1):2:(m-1)]' / (m+modulo(m,2)));
endfunction
+
+//test input validation:
+//assert_checkerror("triang()", "Wrong number of input arguments.");
+//assert_checkerror("triang(1, 2)", "Wrong number of input arguments.");
+//assert_checkerror("triang(0.5)", "parzenwin: M must be a positive integer");
+//assert_checkerror("triang(-1)", "parzenwin: M must be a positive integer");
+//assert_checkerror("triang(zeros (2, 5))", "parzenwin: M must be a positive integer");
+
+//tests:
+//assert_checkequal(triang(1), 1);
+//assert_checkequal(triang(2), [1; 1]/2);
+//assert_checkequal(triang(3), [1; 2; 1]/2);
+//assert_checkequal(triang(4), [1; 3; 3; 1]/4);
diff --git a/macros/tripuls.sci b/macros/tripuls.sci
index 498160b..b4fc2fb 100644
--- a/macros/tripuls.sci
+++ b/macros/tripuls.sci
@@ -1,44 +1,69 @@
-function [y] = tripuls(t,w,skew)
-
-//This function generates a triangular pulse over the interval [-w/2,w/2] sampled at times t
-//Calling Sequence
-//[y] = tripuls(t)
-//[y] = tripuls(t,w)
-//[y] = tripuls(t,w,skew)
-//Parameters
-//t: vector of real or complex numbers
-//w: real or complex number
-//skew: real number, -1 <= s <= 1
-//Description
-//This function generates a triangular pulse which is sampled at times t over the interval [-w/2,w/2]. The value of skew lies between -1
-//and 1.
-//The value of skew represents the relative placement of the peak in the given width.
-//Examples
+function y = tripuls(t, w, skew)
+//This function generates a triangular pulse over the interval [-w/2, w/2] sampled at times t.
+//Calling Sequence:
+//y = tripuls(t)
+//y = tripuls(t, w)
+//y = tripuls(t, w, skew)
+//Parameters:
+//t: Vector of real numbers
+//w: Real number. The default value is 1.
+//skew: Real number in the interval [-1, 1], The default value is 0.
+//Description:
+//This function generates a triangular pulse which is sampled at times t over the interval [-w/2,w/2].
+//The value of skew lies between -1 and 1. The value of skew represents the relative placement of the peak in the given width.
+//Examples:
//tripuls([0, .5, .6, 1], 0.9)
//ans =
// 1 0 0 0
-//This function is being called from Octave
-
-funcprot(0);
-
-rhs = argn(2)
-
-if(rhs<1 | rhs>3)
-error("Wrong number of input arguments.")
-end
- select(rhs)
- case 1 then
- y = callOctave("tripuls",t)
- case 2 then
- y = callOctave("tripuls",t,w)
- case 3 then
- //tip = type(skew)
- //if(tip==1)
- y = callOctave("tripuls",t,w,skew)
- //else
- //error("Wrong arguments.")
- //end
- end
+
+ funcprot(0);
+ rhs= argn(2);
+ if (rhs < 1 | rhs > 3)
+ error("tripuls: wrong number of input arguments");
+ end
+
+ if (rhs == 1)
+ w = 1;
+ skew = 0;
+ else if (rhs == 2)
+ skew = 0;
+ end
+ end
+
+ if (~isreal(w) | ~isscalar(w))
+ error ("tripuls: arg2 (W) must be a real scalar");
+ end
+
+ if (~isreal(skew) | ~isscalar(skew) | skew < -1 | skew > 1)
+ error ("tripuls: arg3 (SKEW) must be a real scalar in the range [-1, 1]");
+ end
+
+ y = zeros(size(t, 'r'), size(t, 'c'));
+ peak = skew * w/2;
+
+ idx = find((t >= -w/2) & (t <= peak));
+ if (idx)
+ y(idx) = (t(idx) + w/2) / (peak + w/2);
+ end
+
+ idx = find((t > peak) & (t < w/2));
+ if (idx)
+ y(idx) = (t(idx) - w/2) / (peak - w/2);
+ end
+
endfunction
+//input validation:
+//assert_checkerror("tripuls()", "tripuls: wrong number of input arguments");
+//assert_checkerror("tripuls(1, 2, 3, 4)", "Wrong number of input arguments.");
+//assert_checkerror("tripuls(1, 2*%i)", "tripuls: arg2 (W) must be a real scalar");
+//assert_checkerror("tripuls(1, 2, 2)", "tripuls: arg3 (SKEW) must be a real scalar in the range [-1, 1]");
+//assert_checkerror("tripuls(1, 2, -2)", "tripuls: arg3 (SKEW) must be a real scalar in the range [-1, 1]");
+//tests:
+//assert_checkequal(tripuls([]), []);
+//assert_checkequal(tripuls([], 0.1), []);
+//assert_checkequal(tripuls(zeros (10, 1)), ones (10, 1));
+//assert_checkequal(tripuls(-1:1), [0, 1, 0]);
+//assert_checkequal(tripuls(-5:5, 9), [0, 1, 3, 5, 7, 9, 7, 5, 3, 1, 0] / 9);
+//assert_checkequal(tripuls(0:1/100:0.3, 0.1), tripuls([0:1/100:0.3]', 0.1)');
diff --git a/macros/zp2sos.sci b/macros/zp2sos.sci
index 6bb9426..f299099 100644
--- a/macros/zp2sos.sci
+++ b/macros/zp2sos.sci
@@ -1,4 +1,135 @@
-function [sos,g] = zp2sos(z,p,k)
+function B = ipermute(A, perm)
+ funcprot(0);
+ // ipermute : Inverse permute the dimensions of a matrix A.
+ // B = ipermute(A, perm) returns the array A with dimensions inverted
+ // according to the permutation vector `perm`.
+ // Validate the permutation vector
+ if max(size(perm)) ~= ndims(A) || or(gsort(perm, "g", "i") ~= 1:ndims(A))
+ error('Permutation vector must contain unique integers from 1 to ndims(A).');
+ end
+ // Compute the inverse permutation vector
+ invPerm = zeros(size(perm,1),size(perm , 2));
+ for i = 1:max(size(perm))
+ invPerm(perm(i)) = i;
+ end
+ // Use the permute function with the inverse permutation
+ B = permute(A, invPerm);
+endfunction
+
+function zsort = cplxpair (z, tol, dim)
+ funcprot(0);
+ if (nargin < 1)
+ error("Invalid inputs");
+ end
+ // default double
+ realmin = 2.2251e-308
+ if (isempty (z))
+ zsort = zeros (size (z,1) , size (z,2));
+ return;
+ end
+ if (nargin < 2 || isempty (tol))
+ tol = 100* %eps;
+ elseif (~ isscalar (tol) || tol < 0 || tol >= 1)
+ error ("cplxpair: TOL must be a scalar number in the range 0 <= TOL < 1");
+ end
+ nd = ndims (z);
+ if (nargin < 3)
+ // Find the first singleton dimension.
+ sz = size (z);
+ dim = find (sz > 1, 1);
+ if isempty(dim)
+ dim = 1;
+ end
+ else
+ dim = floor (dim);
+ if (dim < 1 || dim > nd)
+ error ("cplxpair: invalid dimension DIM");
+ end
+ end
+ // Move dimension to analyze to first position, and convert to a 2-D matrix.
+ perm = [dim:nd, 1:dim-1];
+ z = permute (z, perm);
+ sz = size (z);
+ n = sz(1);
+ m = prod (sz) / n;
+ z = matrix (z, n, m);
+ // Sort the sequence in terms of increasing real values.
+ [temp, idx] = gsort (real (z), 1 , "i");
+ z = z(idx + n * ones (n, 1) * [0:m-1]);
+ // Put the purely real values at the end of the returned list.
+ [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin) <= tol);
+ // Force values detected to be real within tolerance to actually be real.
+ z(idxi + n*(idxj-1)) = real (z(idxi + n*(idxj-1)));
+ //if idxi and idxj are not scalers
+ if ~isscalar(idxi) then
+ v = ones(size(idxi,1),size(idxi,2));
+ else
+ v = 1 ;
+ end
+ q = sparse ([idxi' idxj'], v, [n m]);
+ nr = sum (q, 1);
+ [temp, idx] = gsort (q, 'r','i');
+ midx = idx + size (idx,1) * ones (size (idx,1), 1) * [0:size(idx,2)-1];
+ z = z(midx);
+ zsort = z;
+ // For each remaining z, place the value and its conjugate at the start of
+ // the returned list, and remove them from further consideration.
+ for j = 1:m
+ p = n - nr(j);
+ for i = 1:2:p
+ if (i+1 > p)
+ error ("cplxpair: could not pair all complex numbers");
+ end
+ [v, idx] = min (abs (z(i+1:p,j) - conj (z(i,j))));
+ if (v >= tol * abs (z(i,j)))
+ error ("cplxpair: could not pair all complex numbers");
+ end
+ // For pairs, select the one with positive imaginary part and use it and
+ // it's conjugate, but list the negative imaginary pair first.
+ if (imag (z(i,j)) > 0)
+ zsort([i, i+1],j) = [conj(z(i,j)), z(i,j)];
+ else
+ zsort([i, i+1],j) = [conj(z(idx+i,j)), z(idx+i,j)];
+ end
+ z(idx+i,j) = z(i+1,j);
+ end
+ end
+ // Reshape the output matrix.
+ zsort = ipermute (matrix (zsort, sz), perm);
+endfunction
+
+function [zc, zr] = cplxreal (z, tol, dim)
+ funcprot(0);
+ if (nargin < 1 || nargin > 3)
+ error("invalid inputs");
+ end
+ if (isempty (z))
+ zc = zeros (size (z,1),size(z,2));
+ zr = zeros (size (z,1),size(z,2));
+ return;
+ end
+ if (nargin < 2 || isempty (tol))
+ tol = 100 * %eps ;
+ end
+ if (nargin >= 3)
+ zcp = cplxpair(z,tol,dim);
+ else
+ zcp = cplxpair (z , tol);
+ end
+ nz = max(size (z) );
+ idx = nz;
+ while ((idx > 0) && (zcp(idx) == 0 || (abs (imag (zcp(idx))) ./ abs (zcp(idx))) <= tol))
+ zcp(idx) = real (zcp(idx));
+ idx = idx - 1;
+ end
+ if (pmodulo (idx, 2) ~= 0)
+ error ("cplxreal: odd number of complex values was returned from cplxpair");
+ end
+ zc = zcp(2:2:idx);
+ zr = zcp(idx+1:nz);
+endfunction
+
+function [SOS, G] = zp2sos(z, p, k, DoNotCombineReal)
//This function converts filter poles and zeros to second-order sections.
//Calling Sequence
//[sos] = zp2sos(z)
@@ -10,7 +141,6 @@ function [sos,g] = zp2sos(z,p,k)
//p: column vector
//k: real or complex value, default value is 1
//Description
-//This is an Octave function.
//This function converts filter poles and zeros to second-order sections.
//The first and second parameters are column vectors containing zeros and poles. The third parameter is the overall filter gain, the default value of which is 1.
//The output is the sos matrix and the overall gain.
@@ -20,32 +150,133 @@ function [sos,g] = zp2sos(z,p,k)
//ans =
// 6 -18 12 1 -2 0
// 1 -3 0 1 0 0
+
+ if argn(2) < 3 then
+ k = 1;
+ end
+ if argn(2) < 2 then
+ p = [];
+ end
+
+ DoNotCombineReal = 0;
+
+ [zc, zr] = cplxreal(z(:));
+ [pc, pr] = cplxreal(p(:));
+
+ nzc = length(zc);
+ npc = length(pc);
+
+ nzr = length(zr);
+ npr = length(pr);
+
+ if DoNotCombineReal then
+
+ // Handling complex conjugate poles
+ for count = 1:npc
+ SOS(count, 4:6) = [1, -2 * real(pc(count)), abs(pc(count))^2];
+ end
+
+ // Handling real poles
+ for count = 1:npr
+ SOS(count + npc, 4:6) = [0, 1, -pr(count)];
+ end
+ // Handling complex conjugate zeros
+ for count = 1:nzc
+ SOS(count, 1:3) = [1, -2 * real(zc(count)), abs(zc(count))^2];
+ end
-funcprot(0);
-rhs = argn(2)
-lhs = argn(1)
-if(rhs<1 | rhs>3)
-error("Wrong number of input arguments.")
-end
- select(rhs)
- case 1 then
- if (lhs<2)
- sos=callOctave("zp2sos",z)
- else
- [sos,g]=callOctave("zp2sos",z)
- end
- case 2 then
- if(lhs<2)
- [sos]=callOctave("zp2sos",z,p)
- else
- [sos,g]=callOctave("zp2sos",z,p)
- end
- case 3 then
- if(lhs<2)
- sos=callOctave("zp2sos",z,p,k)
- else
- [sos,g]=callOctave("zp2sos",z,p,k)
- end
- end
+ // Handling real zeros
+ for count = 1:nzr
+ SOS(count + nzc, 1:3) = [0, 1, -zr(count)];
+ end
+
+ // Completing SOS if needed (sections without pole or zero)
+ if npc + npr > nzc + nzr then
+ for count = nzc + nzr + 1 : npc + npr // sections without zero
+ SOS(count, 1:3) = [0, 0, 1];
+ end
+ else
+ for count = npc + npr + 1 : nzc + nzr // sections without pole
+ SOS(count, 4:6) = [0, 0, 1];
+ end
+ end
+
+ else
+
+ // Handling complex conjugate poles
+ for count = 1:npc
+ SOS(count, 4:6) = [1, -2 * real(pc(count)), abs(pc(count))^2];
+ end
+
+ // Handling pair of real poles
+ for count = 1:floor(npr / 2)
+ SOS(count + npc, 4:6) = [1, -pr(2 * count - 1) - pr(2 * count), pr(2 * count - 1) * pr(2 * count)];
+ end
+
+ // Handling last real pole (if any)
+ if pmodulo(npr, 2) == 1 then
+ SOS(npc + floor(npr / 2) + 1, 4:6) = [0, 1, -pr($)];
+ end
+
+ // Handling complex conjugate zeros
+ for count = 1:nzc
+ SOS(count, 1:3) = [1, -2 * real(zc(count)), abs(zc(count))^2];
+ end
+
+ // Handling pair of real zeros
+ for count = 1:floor(nzr / 2)
+ SOS(count + nzc, 1:3) = [1, -zr(2 * count - 1) - zr(2 * count), zr(2 * count - 1) * zr(2 * count)];
+ end
+
+ // Handling last real zero (if any)
+ if pmodulo(nzr, 2) == 1 then
+ SOS(nzc + floor(nzr / 2) + 1, 1:3) = [0, 1, -zr($)];
+ end
+
+ // Completing SOS if needed (sections without pole or zero)
+ if npc + ceil(npr / 2) > nzc + ceil(nzr / 2) then
+ for count = nzc + ceil(nzr / 2) + 1 : npc + ceil(npr / 2) // sections without zero
+ SOS(count, 1:3) = [0, 0, 1];
+ end
+ else
+ for count = npc + ceil(npr / 2) + 1 : nzc + ceil(nzr / 2) // sections without pole
+ SOS(count, 4:6) = [0, 0, 1];
+ end
+ end
+ end
+
+ if ~exists('SOS') then
+ SOS = [0, 0, 1, 0, 0, 1]; // leading zeros will be removed
+ end
+
+ // Removing leading zeros if present in numerator and denominator
+ for count = 1:size(SOS, 1)
+ B = SOS(count, 1:3);
+ A = SOS(count, 4:6);
+ while B(1) == 0 & A(1) == 0 do
+ A(1) = [];
+ A($ + 1) = 0;
+ B(1) = [];
+ B($ + 1) = 0;
+ end
+ SOS(count, :) = [B, A];
+ end
+
+ // If no output argument for the overall gain, combine it into the first section.
+ if argn(1) < 2 then
+ SOS(1, 1:3) = k * SOS(1, 1:3);
+ else
+ G = k;
+ end
endfunction
+
+//tests
+//sos = zp2sos ([]);
+//sos = zp2sos ([], []);
+//sos = zp2sos ([], [], 2);
+//[sos, g] = zp2sos ([], [], 2);
+//sos = zp2sos([], [0], 1);
+//sos = zp2sos([0], [], 1);
+//sos = zp2sos([1,2,3,4,5,6], 2);
+//sos = zp2sos([-1-%i, -1+%i], [-1-2*%i, -1+2*%i], 10);
generated by cgit (git 2.34.1) at 2024-12-19 17:04:28 +0000