1 # Created by Octave 3.6.1, Thu May 17 20:38:50 2012 UTC <root@brouzouf>
13 # name: <cell-element>
18 [psd,f_out] = ar_psd(a,v,freq,Fs,range,method,plot_type)
20 Calculate the power spectrum of the autoregressive model
23 x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
25 where x(n) is the output of the model and e(n) is white noise.
26 This function is intended for use with
27 [a,v,k] = arburg(x,poles,criterion)
28 which use the Burg (1968) method to calculate a "maximum entropy"
29 autoregressive model of "x". This function runs on octave and matlab.
31 If the "freq" argument is a vector (of frequencies) the spectrum is
32 calculated using the polynomial method and the "method" argument is
33 ignored. For scalar "freq", an integer power of 2, or "method='FFT'",
34 causes the spectrum to be calculated by FFT. Otherwise, the spectrum
35 is calculated as a polynomial. It may be computationally more
36 efficient to use the FFT method if length of the model is not much
37 smaller than the number of frequency values. The spectrum is scaled so
38 that spectral energy (area under spectrum) is the same as the
39 time-domain energy (mean square of the signal).
42 All but the first two arguments are optional and may be empty.
44 a %% [vector] list of M=(order+1) autoregressive model
45 %% coefficients. The first element of "ar_coeffs" is the
46 %% zero-lag coefficient, which always has a value of 1.
48 v %% [real scalar] square of the moving-average coefficient of
51 freq %% [real vector] frequencies at which power spectral density
53 %% [integer scalar] number of uniformly distributed frequency
54 %% values at which spectral density is calculated.
57 Fs %% [real scalar] sampling frequency (Hertz) [default=1]
59 CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
60 Control-string arguments can be in any order after the other arguments.
62 range %% 'half', 'onesided' : frequency range of the spectrum is
63 %% from zero up to but not including sample_f/2. Power
64 %% from negative frequencies is added to the positive
65 %% side of the spectrum.
66 %% 'whole', 'twosided' : frequency range of the spectrum is
67 %% -sample_f/2 to sample_f/2, with negative frequencies
68 %% stored in "wrap around" order after the positive
69 %% frequencies; e.g. frequencies for a 10-point 'twosided'
70 %% spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
71 %% 'shift', 'centerdc' : same as 'whole' but with the first half
72 %% of the spectrum swapped with second half to put the
73 %% zero-frequency value in the middle. (See "help
74 %% fftshift". If "freq" is vector, 'shift' is ignored.
75 %% If model coefficients "ar_coeffs" are real, the default
76 %% range is 'half', otherwise default range is 'whole'.
78 method %% 'fft': use FFT to calculate power spectrum.
79 %% 'poly': calculate power spectrum as a polynomial of 1/z
80 %% N.B. this argument is ignored if the "freq" argument is a
81 %% vector. The default is 'poly' unless the "freq"
82 %% argument is an integer power of 2.
84 plot_type%% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
85 %% specifies the type of plot. The default is 'plot', which
86 %% means linear-linear axes. 'squared' is the same as 'plot'.
87 %% 'dB' plots "10*log10(psd)". This argument is ignored and a
88 %% spectrum is not plotted if the caller requires a returned
92 If returned values are not required by the caller, the spectrum
93 is plotted and nothing is returned.
94 psd %% [real vector] estimate of power-spectral density
95 f_out %% [real vector] frequency values
97 N.B. arburg runs in octave and matlab, and does not depend on octave-forge
98 or signal-processing-toolbox functions.
101 [1] Equation 2.28 from Steven M. Kay and Stanley Lawrence Marple Jr.:
102 "Spectrum analysis -- a modern perspective",
103 Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
108 # name: <cell-element>
113 [psd,f_out] = ar_psd(a,v,freq,Fs,range,method,plot_type)
118 # name: <cell-element>
125 # name: <cell-element>
129 [a,v,k] = arburg(x,poles,criterion)
131 Calculate coefficients of an autoregressive (AR) model of complex data
132 "x" using the whitening lattice-filter method of Burg (1968). The inverse
133 of the model is a moving-average filter which reduces "x" to white noise.
134 The power spectrum of the AR model is an estimate of the maximum
135 entropy power spectrum of the data. The function "ar_psd" calculates the
136 power spectrum of the AR model.
139 x %% [vector] sampled data
141 poles %% [integer scalar] number of poles in the AR model or
142 %% limit to the number of poles if a
143 %% valid "stop_crit" is provided.
145 criterion %% [optional string arg] model-selection criterion. Limits
146 %% the number of poles so that spurious poles are not
147 %% added when the whitened data has no more information
148 %% in it (see Kay & Marple, 1981). Recognised values are
149 %% 'AKICc' -- approximate corrected Kullback information
150 %% criterion (recommended),
151 %% 'KIC' -- Kullback information criterion
152 %% 'AICc' -- corrected Akaike information criterion
153 %% 'AIC' -- Akaike information criterion
154 %% 'FPE' -- final prediction error" criterion
155 %% The default is to NOT use a model-selection criterion
158 a %% [polynomial/vector] list of (P+1) autoregression coeffic-
159 %% ients; for data input x(n) and white noise e(n),
162 %% x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
165 v %% [real scalar] mean square of residual noise from the
166 %% whitening operation of the Burg lattice filter.
168 k %% [column vector] reflection coefficients defining the
169 %% lattice-filter embodiment of the model
172 (1) arburg does not remove the mean from the data. You should remove
173 the mean from the data if you want a power spectrum. A non-zero mean
174 can produce large errors in a power-spectrum estimate. See
176 (2) If you don't know what the value of "poles" should be, choose the
177 largest (reasonable) value you could want and use the recommended
178 value, criterion='AKICc', so that arburg can find it.
179 E.g. arburg(x,64,'AKICc')
180 The AKICc has the least bias and best resolution of the available
181 model-selection criteria.
182 (3) arburg runs in octave and matlab, does not depend on octave forge
183 or signal-processing-toolbox functions.
184 (4) Autoregressive and moving-average filters are stored as polynomials
185 which, in matlab, are row vectors.
187 NOTE ON SELECTION CRITERION
188 AIC, AICc, KIC and AKICc are based on information theory. They attempt
189 to balance the complexity (or length) of the model against how well the
190 model fits the data. AIC and KIC are biassed estimates of the asymmetric
191 and the symmetric Kullback-Leibler divergence respectively. AICc and
192 AKICc attempt to correct the bias. See reference [4].
196 [1] John Parker Burg (1968)
197 "A new analysis technique for time series data",
198 NATO advanced study Institute on Signal Processing with Emphasis on
199 Underwater Acoustics, Enschede, Netherlands, Aug. 12-23, 1968.
201 [2] Steven M. Kay and Stanley Lawrence Marple Jr.:
202 "Spectrum analysis -- a modern perspective",
203 Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
205 [3] William H. Press and Saul A. Teukolsky and William T. Vetterling and
207 "Numerical recipes in C, The art of scientific computing", 2nd edition,
208 Cambridge University Press, 2002 --- Section 13.7.
210 [4] Abd-Krim Seghouane and Maiza Bekara
211 "A small sample model selection criterion based on Kullback's symmetric
212 divergence", IEEE Transactions on Signal Processing,
213 Vol. 52(12), pp 3314-3323, Dec. 2004
217 # name: <cell-element>
221 [a,v,k] = arburg(x,poles,criterion)
226 # name: <cell-element>
233 # name: <cell-element>
237 usage: [a, v, k] = aryule (x, p)
239 fits an AR (p)-model with Yule-Walker estimates.
240 x = data vector to estimate
242 v: variance of white noise
243 k: reflection coeffients for use in lattice filter
245 The power spectrum of the resulting filter can be plotted with
246 pyulear(x, p), or you can plot it directly with ar_psd(a,v,...).
249 pyulear, power, freqz, impz -- for observing characteristics of the model
250 arburg -- for alternative spectral estimators
252 Example: Use example from arburg, but substitute aryule for arburg.
254 Note: Orphanidis '85 claims lattice filters are more tolerant of
255 truncation errors, which is why you might want to use them. However,
256 lacking a lattice filter processor, I haven't tested that the lattice
257 filter coefficients are reasonable.
261 # name: <cell-element>
265 usage: [a, v, k] = aryule (x, p)
267 fits an AR (p)-model with Yule-Walker esti
271 # name: <cell-element>
278 # name: <cell-element>
282 -- Function File: [W] = barthannwin(L)
283 Compute the modified Bartlett-Hann window of lenght L.
285 See also: rectwin, bartlett
291 # name: <cell-element>
295 Compute the modified Bartlett-Hann window of lenght L.
299 # name: <cell-element>
306 # name: <cell-element>
310 Return bessel analog filter prototype.
314 http://en.wikipedia.org/wiki/Bessel_polynomials
318 # name: <cell-element>
322 Return bessel analog filter prototype.
326 # name: <cell-element>
333 # name: <cell-element>
337 Generate a bessel filter.
338 Default is a Laplace space (s) filter.
340 [b,a] = besself(n, Wc)
341 low pass filter with cutoff pi*Wc radians
343 [b,a] = besself(n, Wc, 'high')
344 high pass filter with cutoff pi*Wc radians
346 [b,a] = besself(n, [Wl, Wh])
347 band pass filter with edges pi*Wl and pi*Wh radians
349 [b,a] = besself(n, [Wl, Wh], 'stop')
350 band reject filter with edges pi*Wl and pi*Wh radians
352 [z,p,g] = besself(...)
353 return filter as zero-pole-gain rather than coefficients of the
354 numerator and denominator polynomials.
356 [...] = besself(...,'z')
357 return a discrete space (Z) filter, W must be less than 1.
359 [a,b,c,d] = besself(...)
360 return state-space matrices
364 Proakis & Manolakis (1992). Digital Signal Processing. New York:
365 Macmillan Publishing Company.
369 # name: <cell-element>
373 Generate a bessel filter.
377 # name: <cell-element>
384 # name: <cell-element>
388 usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
389 [Zb, Za] = bilinear(Sb, Sa, T)
391 Transform a s-plane filter specification into a z-plane
392 specification. Filters can be specified in either zero-pole-gain or
393 transfer function form. The input form does not have to match the
394 output form. 1/T is the sampling frequency represented in the z plane.
396 Note: this differs from the bilinear function in the signal processing
397 toolbox, which uses 1/T rather than T.
399 Theory: Given a piecewise flat filter design, you can transform it
400 from the s-plane to the z-plane while maintaining the band edges by
401 means of the bilinear transform. This maps the left hand side of the
402 s-plane into the interior of the unit circle. The mapping is highly
403 non-linear, so you must design your filter with band edges in the
404 s-plane positioned at 2/T tan(w*T/2) so that they will be positioned
405 at w after the bilinear transform is complete.
407 The following table summarizes the transformation:
409 +---------------+-----------------------+----------------------+
410 | Transform | Zero at x | Pole at x |
411 | H(S) | H(S) = S-x | H(S)=1/(S-x) |
412 +---------------+-----------------------+----------------------+
413 | 2 z-1 | zero: (2+xT)/(2-xT) | zero: -1 |
414 | S -> - --- | pole: -1 | pole: (2+xT)/(2-xT) |
415 | T z+1 | gain: (2-xT)/T | gain: (2-xT)/T |
416 +---------------+-----------------------+----------------------+
418 With tedious algebra, you can derive the above formulae yourself by
419 substituting the transform for S into H(S)=S-x for a zero at x or
420 H(S)=1/(S-x) for a pole at x, and converting the result into the
423 H(Z)=g prod(Z-Xi)/prod(Z-Xj)
425 Please note that a pole and a zero at the same place exactly cancel.
426 This is significant since the bilinear transform creates numerous
427 extra poles and zeros, most of which cancel. Those which do not
428 cancel have a "fill-in" effect, extending the shorter of the sets to
429 have the same number of as the longer of the sets of poles and zeros
430 (or at least split the difference in the case of the band pass
431 filter). There may be other opportunistic cancellations but I will
434 Also note that any pole on the unit circle or beyond will result in
435 an unstable filter. Because of cancellation, this will only happen
436 if the number of poles is smaller than the number of zeros. The
437 analytic design methods all yield more poles than zeros, so this will
442 Proakis & Manolakis (1992). Digital Signal Processing. New York:
443 Macmillan Publishing Company.
447 # name: <cell-element>
451 usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
452 [Zb, Za] = bilinear(Sb, S
456 # name: <cell-element>
463 # name: <cell-element>
467 -- Function File: [Y I] = bitrevorder(X)
468 Reorder x in the bit reversed order
476 # name: <cell-element>
480 Reorder x in the bit reversed order
485 # name: <cell-element>
492 # name: <cell-element>
496 -- Function File: [W] = blackmanharris(L)
497 Compute the Blackman-Harris window.
499 See also: rectwin, bartlett
505 # name: <cell-element>
509 Compute the Blackman-Harris window.
513 # name: <cell-element>
520 # name: <cell-element>
524 -- Function File: [W] = blackmannuttall(L)
525 Compute the Blackman-Nuttall window.
527 See also: nuttallwin, kaiser
533 # name: <cell-element>
537 Compute the Blackman-Nuttall window.
541 # name: <cell-element>
548 # name: <cell-element>
552 -- Function File: [W] = bohmanwin(L)
553 Compute the Bohman window of lenght L.
555 See also: rectwin, bartlett
561 # name: <cell-element>
565 Compute the Bohman window of lenght L.
569 # name: <cell-element>
576 # name: <cell-element>
580 usage: w = boxcar (n)
582 Returns the filter coefficients of a rectangular window of length n.
586 # name: <cell-element>
590 usage: w = boxcar (n)
595 # name: <cell-element>
602 # name: <cell-element>
606 -- Function File: Y = buffer (X, N, P, OPT)
607 -- Function File: [Y, Z, OPT] = buffer (...)
608 Buffer a signal into a data frame. The arguments to `buffer' are
611 The data to be buffered.
614 The number of rows in the produced data buffer. This is an
615 positive integer value and must be supplied.
618 An integer less than N that specifies the under- or overlap
619 between column in the data frame. The default value of P is 0.
622 In the case of an overlap, OPT can be either a vector of
623 length P or the string 'nodelay'. If OPT is a vector, then the
624 first P entries in Y will be filled with these values. If OPT
625 is the string 'nodelay', then the first value of Y
626 corresponds to the first value of X.
628 In the can of an underlap, OPT must be an integer between 0
629 and `-P'. The represents the initial underlap of the first
632 The default value for OPT the vector `zeros (1, P)' in the
633 case of an overlap, or 0 otherwise.
635 In the case of a single output argument, Y will be padded with
636 zeros to fill the missing values in the data frame. With two output
637 arguments Z is the remaining data that has not been used in the
640 Likewise, the output OPT is the overlap, or underlap that might be
641 used for a future call to `code' to allow continuous buffering.
646 # name: <cell-element>
650 Buffer a signal into a data frame.
654 # name: <cell-element>
661 # name: <cell-element>
665 Generate a butterworth filter.
666 Default is a discrete space (Z) filter.
668 [b,a] = butter(n, Wc)
669 low pass filter with cutoff pi*Wc radians
671 [b,a] = butter(n, Wc, 'high')
672 high pass filter with cutoff pi*Wc radians
674 [b,a] = butter(n, [Wl, Wh])
675 band pass filter with edges pi*Wl and pi*Wh radians
677 [b,a] = butter(n, [Wl, Wh], 'stop')
678 band reject filter with edges pi*Wl and pi*Wh radians
680 [z,p,g] = butter(...)
681 return filter as zero-pole-gain rather than coefficients of the
682 numerator and denominator polynomials.
684 [...] = butter(...,'s')
685 return a Laplace space filter, W can be larger than 1.
687 [a,b,c,d] = butter(...)
688 return state-space matrices
692 Proakis & Manolakis (1992). Digital Signal Processing. New York:
693 Macmillan Publishing Company.
697 # name: <cell-element>
701 Generate a butterworth filter.
705 # name: <cell-element>
712 # name: <cell-element>
716 Compute butterworth filter order and cutoff for the desired response
717 characteristics. Rp is the allowable decibels of ripple in the pass
718 band. Rs is the minimum attenuation in the stop band.
720 [n, Wc] = buttord(Wp, Ws, Rp, Rs)
721 Low pass (Wp<Ws) or high pass (Wp>Ws) filter design. Wp is the
722 pass band edge and Ws is the stop band edge. Frequencies are
723 normalized to [0,1], corresponding to the range [0,Fs/2].
725 [n, Wc] = buttord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
726 Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
727 filter design. Wp gives the edges of the pass band, and Ws gives
728 the edges of the stop band.
730 Theory: |H(W)|^2 = 1/[1+(W/Wc)^(2N)] = 10^(-R/10)
731 With some algebra, you can solve simultaneously for Wc and N given
732 Ws,Rs and Wp,Rp. For high pass filters, subtracting the band edges
733 from Fs/2, performing the test, and swapping the resulting Wc back
734 works beautifully. For bandpass and bandstop filters this process
735 significantly overdesigns. Artificially dividing N by 2 in this case
736 helps a lot, but it still overdesigns.
742 # name: <cell-element>
746 Compute butterworth filter order and cutoff for the desired response
751 # name: <cell-element>
758 # name: <cell-element>
762 usage: cceps (x [, correct])
764 Returns the complex cepstrum of the vector x.
765 If the optional argument correct has the value 1, a correction
766 method is applied. The default is not to do this.
770 # name: <cell-element>
774 usage: cceps (x [, correct])
779 # name: <cell-element>
786 # name: <cell-element>
792 Returns the value of the nth-order Chebyshev polynomial calculated at
793 the point x. The Chebyshev polynomials are defined by the equations:
795 / cos(n acos(x), |x| <= 1
797 \ cosh(n acosh(x), |x| > 1
799 If x is a vector, the output is a vector of the same size, where each
800 element is calculated as y(i) = Tn(x(i)).
804 # name: <cell-element>
813 # name: <cell-element>
820 # name: <cell-element>
824 Compute chebyshev type I filter order and cutoff for the desired response
825 characteristics. Rp is the allowable decibels of ripple in the pass
826 band. Rs is the minimum attenuation in the stop band.
828 [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)
829 Low pass (Wp<Ws) or high pass (Wp>Ws) filter design. Wp is the
830 pass band edge and Ws is the stop band edge. Frequencies are
831 normalized to [0,1], corresponding to the range [0,Fs/2].
833 [n, Wc] = cheb1ord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
834 Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
835 filter design. Wp gives the edges of the pass band, and Ws gives
836 the edges of the stop band.
842 # name: <cell-element>
846 Compute chebyshev type I filter order and cutoff for the desired response
851 # name: <cell-element>
858 # name: <cell-element>
862 Compute chebyshev type II filter order and cutoff for the desired response
863 characteristics. Rp is the allowable decibels of ripple in the pass
864 band. Rs is the minimum attenuation in the stop band.
866 [n, Wc] = cheb2ord(Wp, Ws, Rp, Rs)
867 Low pass (Wp<Ws) or high pass (Wp>Ws) filter design. Wp is the
868 pass band edge and Ws is the stop band edge. Frequencies are
869 normalized to [0,1], corresponding to the range [0,Fs/2].
871 [n, Wc] = cheb2ord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
872 Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
873 filter design. Wp gives the edges of the pass band, and Ws gives
874 the edges of the stop band.
882 # name: <cell-element>
886 Compute chebyshev type II filter order and cutoff for the desired response
891 # name: <cell-element>
898 # name: <cell-element>
902 Usage: chebwin (L, at)
904 Returns the filter coefficients of the L-point Dolph-Chebyshev window
905 with at dB of attenuation in the stop-band of the corresponding
908 For the definition of the Chebyshev window, see
910 * Peter Lynch, "The Dolph-Chebyshev Window: A Simple Optimal Filter",
911 Monthly Weather Review, Vol. 125, pp. 655-660, April 1997.
912 (http://www.maths.tcd.ie/~plynch/Publications/Dolph.pdf)
914 * C. Dolph, "A current distribution for broadside arrays which
915 optimizes the relationship between beam width and side-lobe level",
916 Proc. IEEE, 34, pp. 335-348.
918 The window is described in frequency domain by the expression:
920 Cheb(L-1, beta * cos(pi * k/L))
921 W(k) = -------------------------------
926 beta = cosh(1/(L-1) * acosh(10^(at/20))
928 and Cheb(m,x) denoting the m-th order Chebyshev polynomial calculated
931 Note that the denominator in W(k) above is not computed, and after
932 the inverse Fourier transform the window is scaled by making its
933 maximum value unitary.
939 # name: <cell-element>
943 Usage: chebwin (L, at)
948 # name: <cell-element>
955 # name: <cell-element>
959 Generate an Chebyshev type I filter with Rp dB of pass band ripple.
961 [b, a] = cheby1(n, Rp, Wc)
962 low pass filter with cutoff pi*Wc radians
964 [b, a] = cheby1(n, Rp, Wc, 'high')
965 high pass filter with cutoff pi*Wc radians
967 [b, a] = cheby1(n, Rp, [Wl, Wh])
968 band pass filter with edges pi*Wl and pi*Wh radians
970 [b, a] = cheby1(n, Rp, [Wl, Wh], 'stop')
971 band reject filter with edges pi*Wl and pi*Wh radians
973 [z, p, g] = cheby1(...)
974 return filter as zero-pole-gain rather than coefficients of the
975 numerator and denominator polynomials.
977 [...] = cheby1(...,'s')
978 return a Laplace space filter, W can be larger than 1.
980 [a,b,c,d] = cheby1(...)
981 return state-space matrices
985 Parks & Burrus (1987). Digital Filter Design. New York:
986 John Wiley & Sons, Inc.
990 # name: <cell-element>
994 Generate an Chebyshev type I filter with Rp dB of pass band ripple.
998 # name: <cell-element>
1005 # name: <cell-element>
1009 Generate an Chebyshev type II filter with Rs dB of stop band attenuation.
1011 [b, a] = cheby2(n, Rs, Wc)
1012 low pass filter with cutoff pi*Wc radians
1014 [b, a] = cheby2(n, Rs, Wc, 'high')
1015 high pass filter with cutoff pi*Wc radians
1017 [b, a] = cheby2(n, Rs, [Wl, Wh])
1018 band pass filter with edges pi*Wl and pi*Wh radians
1020 [b, a] = cheby2(n, Rs, [Wl, Wh], 'stop')
1021 band reject filter with edges pi*Wl and pi*Wh radians
1023 [z, p, g] = cheby2(...)
1024 return filter as zero-pole-gain rather than coefficients of the
1025 numerator and denominator polynomials.
1027 [...] = cheby2(...,'s')
1028 return a Laplace space filter, W can be larger than 1.
1030 [a,b,c,d] = cheby2(...)
1031 return state-space matrices
1035 Parks & Burrus (1987). Digital Filter Design. New York:
1036 John Wiley & Sons, Inc.
1040 # name: <cell-element>
1044 Generate an Chebyshev type II filter with Rs dB of stop band attenuation.
1048 # name: <cell-element>
1055 # name: <cell-element>
1059 usage: y = chirp(t [, f0 [, t1 [, f1 [, form [, phase]]]]])
1061 Evaluate a chirp signal at time t. A chirp signal is a frequency
1064 t: vector of times to evaluate the chirp signal
1065 f0: frequency at time t=0 [ 0 Hz ]
1066 t1: time t1 [ 1 sec ]
1067 f1: frequency at time t=t1 [ 100 Hz ]
1068 form: shape of frequency sweep
1069 'linear' f(t) = (f1-f0)*(t/t1) + f0
1070 'quadratic' f(t) = (f1-f0)*(t/t1)^2 + f0
1071 'logarithmic' f(t) = (f1-f0)^(t/t1) + f0
1072 phase: phase shift at t=0
1075 specgram(chirp([0:0.001:5])); # linear, 0-100Hz in 1 sec
1076 specgram(chirp([-2:0.001:15], 400, 10, 100, 'quadratic'));
1077 soundsc(chirp([0:1/8000:5], 200, 2, 500, "logarithmic"),8000);
1079 If you want a different sweep shape f(t), use the following:
1080 y = cos(2*pi*integral(f(t)) + 2*pi*f0*t + phase);
1084 # name: <cell-element>
1088 usage: y = chirp(t [, f0 [, t1 [, f1 [, form [, phase]]]]])
1093 # name: <cell-element>
1100 # name: <cell-element>
1104 -- Function File: [PSI,X] = cmorwavf (LB,UB,N,FB,FC)
1105 Compute the Complex Morlet wavelet.
1110 # name: <cell-element>
1114 Compute the Complex Morlet wavelet.
1118 # name: <cell-element>
1125 # name: <cell-element>
1130 [Pxx,freq] = cohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
1132 Estimate (mean square) coherence of signals "x" and "y".
1133 Use the Welch (1967) periodogram/FFT method.
1134 Compatible with Matlab R11 cohere and earlier.
1135 See "help pwelch" for description of arguments, hints and references
1136 --- especially hint (7) for Matlab R11 defaults.
1141 # name: <cell-element>
1146 [Pxx,freq] = cohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detren
1150 # name: <cell-element>
1157 # name: <cell-element>
1161 -- Function File: convmtx (A, N)
1162 If A is a column vector and X is a column vector of length N, then
1166 gives the convolution of of A and X and is the same as `conv(A,
1167 X)'. The difference is if many vectors are to be convolved with
1168 the same vector, then this technique is possibly faster.
1170 Similarly, if A is a row vector and X is a row vector of length N,
1175 is the same as `conv(X, A)'.
1182 # name: <cell-element>
1186 If A is a column vector and X is a column vector of length N, then
1191 # name: <cell-element>
1198 # name: <cell-element>
1202 -- Function File: [ZC, ZR] = cplxreal (Z, THRESH)
1203 Split the vector z into its complex (ZC) and real (ZR) elements,
1204 eliminating one of each complex-conjugate pair.
1207 * Z = row- or column-vector of complex numbers
1208 * THRESH = tolerance threshold for numerical comparisons
1212 * ZC = elements of Z having positive imaginary parts
1213 * ZR = elements of Z having zero imaginary part
1215 Each complex element of Z is assumed to have a complex-conjugate
1216 counterpart elsewhere in Z as well. Elements are declared real if
1217 their imaginary parts have magnitude less than THRESH.
1225 # name: <cell-element>
1229 Split the vector z into its complex (ZC) and real (ZR) elements,
1234 # name: <cell-element>
1241 # name: <cell-element>
1246 [Pxx,freq] = cpsd(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
1248 Estimate cross power spectrum of data "x" and "y" by the Welch (1967)
1249 periodogram/FFT method.
1250 See "help pwelch" for description of arguments, hints and references
1254 # name: <cell-element>
1259 [Pxx,freq] = cpsd(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
1263 # name: <cell-element>
1270 # name: <cell-element>
1275 [Pxx,freq] = csd(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
1277 Estimate cross power spectrum of data "x" and "y" by the Welch (1967)
1278 periodogram/FFT method. Compatible with Matlab R11 csd and earlier.
1279 See "help pwelch" for description of arguments, hints and references
1280 --- especially hint (7) for Matlab R11 defaults.
1284 # name: <cell-element>
1289 [Pxx,freq] = csd(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
1294 # name: <cell-element>
1301 # name: <cell-element>
1305 usage y=czt(x, m, w, a)
1307 Chirp z-transform. Compute the frequency response starting at a and
1308 stepping by w for m steps. a is a point in the complex plane, and
1309 w is the ratio between points in each step (i.e., radius increases
1310 exponentially, and angle increases linearly).
1312 To evaluate the frequency response for the range f1 to f2 in a signal
1313 with sampling frequency Fs, use the following:
1314 m = 32; ## number of points desired
1315 w = exp(-j*2*pi*(f2-f1)/((m-1)*Fs)); ## freq. step of f2-f1/m
1316 a = exp(j*2*pi*f1/Fs); ## starting at frequency f1
1317 y = czt(x, m, w, a);
1319 If you don't specify them, then the parameters default to a fourier
1321 m=length(x), w=exp(-j*2*pi/m), a=1
1323 If x is a matrix, the transform will be performed column-by-column.
1327 # name: <cell-element>
1331 usage y=czt(x, m, w, a)
1336 # name: <cell-element>
1343 # name: <cell-element>
1347 -- Function File: [FHANDLE, FULLNAME] = data2fun (TI, YI)
1348 -- Function File: [ ... ] = data2fun (TI, YI,PROPERTY,VALUE)
1349 Creates a vectorized function based on data samples using
1352 The values given in YI (N-by-k matrix) correspond to evaluations
1353 of the function y(t) at the points TI (N-by-1 matrix). The data
1354 is interpolated and the function handle to the generated
1355 interpolant is returned.
1357 The function accepts property-value pairs described below.
1360 Code is generated and .m file is created. The VALUE contains
1361 the name of the function. The returned function handle is a
1362 handle to that file. If VALUE is empty, then a name is
1363 automatically generated using `tmpnam' and the file is
1364 created in the current directory. VALUE must not have an
1365 extension, since .m will be appended. Numerical value used
1366 in the function are stored in a .mat file with the same name
1370 Type of interpolation. See `interp1'.
1379 # name: <cell-element>
1383 Creates a vectorized function based on data samples using interpolation.
1387 # name: <cell-element>
1394 # name: <cell-element>
1399 Computes the discrete cosine transform of x. If n is given, then
1400 x is padded or trimmed to length n before computing the transform.
1401 If x is a matrix, compute the transform along the columns of the
1402 the matrix. The transform is faster if x is real-valued and even
1405 The discrete cosine transform X of x can be defined as follows:
1408 X[k] = w(k) sum x[n] cos (pi (2n+1) k / 2N ), k = 0, ..., N-1
1411 with w(0) = sqrt(1/N) and w(k) = sqrt(2/N), k = 1, ..., N-1. There
1412 are other definitions with different scaling of X[k], but this form
1413 is common in image processing.
1415 See also: idct, dct2, idct2, dctmtx
1419 # name: <cell-element>
1424 Computes the discrete cosine transform of x.
1428 # name: <cell-element>
1435 # name: <cell-element>
1440 Computes the 2-D discrete cosine transform of matrix x
1442 y = dct2 (x, m, n) or y = dct2 (x, [m n])
1443 Computes the 2-D DCT of x after padding or trimming rows to m and
1448 # name: <cell-element>
1453 Computes the 2-D discrete cosine transform of matrix x
1458 # name: <cell-element>
1465 # name: <cell-element>
1470 Return the DCT transformation matrix of size n x n.
1472 If A is an n x n matrix, then the following are true:
1473 T*A == dct(A), T'*A == idct(A)
1474 T*A*T' == dct2(A), T'*A*T == idct2(A)
1476 A dct transformation matrix is useful for doing things like jpeg
1477 image compression, in which an 8x8 dct matrix is applied to
1478 non-overlapping blocks throughout an image and only a subblock on the
1479 top left of each block is kept. During restoration, the remainder of
1480 the block is filled with zeros and the inverse transform is applied
1483 See also: dct, idct, dct2, idct2
1487 # name: <cell-element>
1492 Return the DCT transformation matrix of size n x n.
1496 # name: <cell-element>
1503 # name: <cell-element>
1507 usage: y = decimate(x, q [, n] [, ftype])
1509 Downsample the signal x by a factor of q, using an order n filter
1510 of ftype 'fir' or 'iir'. By default, an order 8 Chebyshev type I
1511 filter is used or a 30 point FIR filter if ftype is 'fir'. Note
1512 that q must be an integer for this rate change method.
1515 ## Generate a signal that starts away from zero, is slowly varying
1516 ## at the start and quickly varying at the end, decimate and plot.
1517 ## Since it starts away from zero, you will see the boundary
1518 ## effects of the antialiasing filter clearly. Next you will see
1519 ## how it follows the curve nicely in the slowly varying early
1520 ## part of the signal, but averages the curve in the quickly
1521 ## varying late part of the signal.
1522 t=0:0.01:2; x=chirp(t,2,.5,10,'quadratic')+sin(2*pi*t*0.4);
1523 y = decimate(x,4); # factor of 4 decimation
1524 stem(t(1:121)*1000,x(1:121),"-g;Original;"); hold on; # plot original
1525 stem(t(1:4:121)*1000,y(1:31),"-r;Decimated;"); hold off; # decimated
1529 # name: <cell-element>
1533 usage: y = decimate(x, q [, n] [, ftype])
1538 # name: <cell-element>
1545 # name: <cell-element>
1549 -- Function File: D = dftmtx (N)
1550 If N is a scalar, produces a N-by-N matrix D such that the Fourier
1551 transform of a column vector of length N is given by `dftmtx(N) *
1552 x' and the inverse Fourier transform is given by `inv(dftmtx(N)) *
1553 x'. In general this is less efficient than calling the "fft" and
1559 # name: <cell-element>
1563 If N is a scalar, produces a N-by-N matrix D such that the Fourier
1568 # name: <cell-element>
1575 # name: <cell-element>
1579 -- Function File: [Y] = diric(X,N)
1580 Compute the dirichlet function.
1582 See also: sinc, gauspuls, sawtooth
1588 # name: <cell-element>
1592 Compute the dirichlet function.
1596 # name: <cell-element>
1603 # name: <cell-element>
1607 -- Function File: Y = downsample (X, N)
1608 -- Function File: Y = downsample (X, N, OFFSET)
1609 Downsample the signal, selecting every nth element. If X is a
1610 matrix, downsample every column.
1612 For most signals you will want to use `decimate' instead since it
1613 prefilters the high frequency components of the signal and avoids
1616 If OFFSET is defined, select every nth element starting at sample
1619 See also: decimate, interp, resample, upfirdn, upsample
1625 # name: <cell-element>
1629 Downsample the signal, selecting every nth element.
1633 # name: <cell-element>
1640 # name: <cell-element>
1644 -- Function File: Y = dst (X)
1645 -- Function File: Y = dst (X, N)
1646 Computes the type I discrete sine transform of X. If N is given,
1647 then X is padded or trimmed to length N before computing the
1648 transform. If X is a matrix, compute the transform along the
1649 columns of the the matrix.
1651 The discrete sine transform X of x can be defined as follows:
1654 X[k] = sum x[n] sin (pi n k / (N+1) ), k = 1, ..., N
1663 # name: <cell-element>
1667 Computes the type I discrete sine transform of X.
1671 # name: <cell-element>
1678 # name: <cell-element>
1682 -- Function File: [CA CD] = dwt(X,LO_D,HI_D)
1683 Comupte de discrete wavelet transform of x with one level.
1688 # name: <cell-element>
1692 Comupte de discrete wavelet transform of x with one level.
1696 # name: <cell-element>
1703 # name: <cell-element>
1708 usage: [Zz, Zp, Zg] = ellip(n, Rp, Rs, Wp, stype,'s')
1710 Generate an Elliptic or Cauer filter (discrete and contnuious).
1712 [b,a] = ellip(n, Rp, Rs, Wp)
1713 low pass filter with order n, cutoff pi*Wp radians, Rp decibels
1714 of ripple in the passband and a stopband Rs decibels down.
1716 [b,a] = ellip(n, Rp, Rs, Wp, 'high')
1717 high pass filter with cutoff pi*Wp...
1719 [b,a] = ellip(n, Rp, Rs, [Wl, Wh])
1720 band pass filter with band pass edges pi*Wl and pi*Wh ...
1722 [b,a] = ellip(n, Rp, Rs, [Wl, Wh], 'stop')
1723 band reject filter with edges pi*Wl and pi*Wh, ...
1725 [z,p,g] = ellip(...)
1726 return filter as zero-pole-gain.
1728 [...] = ellip(...,'s')
1729 return a Laplace space filter, W can be larger than 1.
1731 [a,b,c,d] = ellip(...)
1732 return state-space matrices
1736 - Oppenheim, Alan V., Discrete Time Signal Processing, Hardcover, 1999.
1737 - Parente Ribeiro, E., Notas de aula da disciplina TE498 - Processamento
1738 Digital de Sinais, UFPR, 2001/2002.
1739 - Kienzle, Paul, functions from Octave-Forge, 1999 (http://octave.sf.net).
1743 # name: <cell-element>
1751 # name: <cell-element>
1758 # name: <cell-element>
1762 usage: [n,wp] = ellipord(wp,ws, rp,rs)
1764 Calculate the order for the elliptic filter (discrete)
1765 wp: Cutoff frequency
1767 rp: decibels of ripple in the passband.
1768 rs: decibels of ripple in the stopband.
1772 - Lamar, Marcus Vinicius, Notas de aula da disciplina TE 456 - Circuitos
1773 Analogicos II, UFPR, 2001/2002.
1777 # name: <cell-element>
1781 usage: [n,wp] = ellipord(wp,ws, rp,rs)
1786 # name: <cell-element>
1793 # name: <cell-element>
1797 -- Function File: m = fht ( d, n, dim )
1798 The function fht calculates Fast Hartley Transform where D is
1799 the real input vector (matrix), and M is the real-transform
1800 vector. For matrices the hartley transform is calculated along the
1801 columns by default. The options N,and DIM are similar to the
1802 options of FFT function.
1804 The forward and inverse hartley transforms are the same (except
1805 for a scale factor of 1/N for the inverse hartley transform), but
1806 implemented using different functions .
1808 The definition of the forward hartley transform for vector d,
1810 m[K] = \sum_i=0^N-1 d[i]*(cos[K*2*pi*i/N] + sin[K*2*pi*i/N]), for
1811 0 <= K < N. m[K] = \sum_i=0^N-1 d[i]*CAS[K*i], for 0 <= K < N.
1821 # name: <cell-element>
1825 The function fht calculates Fast Hartley Transform where D is the
1830 # name: <cell-element>
1837 # name: <cell-element>
1841 usage: y = filtfilt(b, a, x)
1843 Forward and reverse filter the signal. This corrects for phase
1844 distortion introduced by a one-pass filter, though it does square the
1845 magnitude response in the process. That's the theory at least. In
1846 practice the phase correction is not perfect, and magnitude response
1847 is distorted, particularly in the stop band.
1850 [b, a]=butter(3, 0.1); % 10 Hz low-pass filter
1851 t = 0:0.01:1.0; % 1 second sample
1852 x=sin(2*pi*t*2.3)+0.25*randn(size(t)); % 2.3 Hz sinusoid+noise
1853 y = filtfilt(b,a,x); z = filter(b,a,x); % apply filter
1854 plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')
1858 # name: <cell-element>
1862 usage: y = filtfilt(b, a, x)
1867 # name: <cell-element>
1874 # name: <cell-element>
1878 Set initial condition vector for filter function
1879 The vector zf has the same values that would be obtained
1880 from function filter given past inputs x and outputs y
1882 The vectors x and y contain the most recent inputs and outputs
1883 respectively, with the newest values first:
1885 x = [x(-1) x(-2) ... x(-nb)], nb = length(b)-1
1886 y = [y(-1) y(-2) ... y(-na)], na = length(a)-a
1888 If length(x)<nb then it is zero padded
1889 If length(y)<na then it is zero padded
1891 zf = filtic(b, a, y)
1892 Initial conditions for filter with coefficients a and b
1893 and output vector y, assuming input vector x is zero
1895 zf = filtic(b, a, y, x)
1896 Initial conditions for filter with coefficients a and b
1897 input vector x and output vector y
1901 # name: <cell-element>
1905 Set initial condition vector for filter function
1906 The vector zf has the same va
1910 # name: <cell-element>
1917 # name: <cell-element>
1921 usage: b = fir1(n, w [, type] [, window] [, noscale])
1923 Produce an order n FIR filter with the given frequency cutoff,
1924 returning the n+1 filter coefficients in b.
1926 n: order of the filter (1 less than the length of the filter)
1928 strictly increasing vector in range [0, 1]
1929 singleton for highpass or lowpass, vector pair for bandpass or
1930 bandstop, or vector for alternating pass/stop filter.
1931 type: choose between pass and stop bands
1932 'high' for highpass filter, cutoff at w
1933 'stop' for bandstop filter, edges at w = [lo, hi]
1934 'DC-0' for bandstop as first band of multiband filter
1935 'DC-1' for bandpass as first band of multiband filter
1936 window: smoothing window
1937 defaults to hamming(n+1) row vector
1938 returned filter is the same shape as the smoothing window
1939 noscale: choose whether to normalize or not
1940 'scale': set the magnitude of the center of the first passband to 1
1941 'noscale': don't normalize
1943 To apply the filter, use the return vector b:
1947 freqz(fir1(40,0.3));
1948 freqz(fir1(15,[0.2, 0.5], 'stop')); # note the zero-crossing at 0.1
1949 freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale'));
1953 # name: <cell-element>
1957 usage: b = fir1(n, w [, type] [, window] [, noscale])
1962 # name: <cell-element>
1969 # name: <cell-element>
1973 usage: b = fir2(n, f, m [, grid_n [, ramp_n]] [, window])
1975 Produce an FIR filter of order n with arbitrary frequency response,
1976 returning the n+1 filter coefficients in b.
1978 n: order of the filter (1 less than the length of the filter)
1979 f: frequency at band edges
1980 f is a vector of nondecreasing elements in [0,1]
1981 the first element must be 0 and the last element must be 1
1982 if elements are identical, it indicates a jump in freq. response
1983 m: magnitude at band edges
1984 m is a vector of length(f)
1985 grid_n: length of ideal frequency response function
1986 defaults to 512, should be a power of 2 bigger than n
1987 ramp_n: transition width for jumps in filter response
1988 defaults to grid_n/20; a wider ramp gives wider transitions
1989 but has better stopband characteristics.
1990 window: smoothing window
1991 defaults to hamming(n+1) row vector
1992 returned filter is the same shape as the smoothing window
1994 To apply the filter, use the return vector b:
1996 Note that plot(f,m) shows target response.
1999 f=[0, 0.3, 0.3, 0.6, 0.6, 1]; m=[0, 0, 1, 1/2, 0, 0];
2000 [h, w] = freqz(fir2(100,f,m));
2001 plot(f,m,';target response;',w/pi,abs(h),';filter response;');
2005 # name: <cell-element>
2009 usage: b = fir2(n, f, m [, grid_n [, ramp_n]] [, window])
2014 # name: <cell-element>
2021 # name: <cell-element>
2026 b = firls(N, F, A, W);
2028 FIR filter design using least squares method. Returns a length N+1
2029 linear phase filter such that the integral of the weighted mean
2030 squared error in the specified bands is minimized.
2032 F specifies the frequencies of the band edges, normalized so that
2033 half the sample frequency is equal to 1. Each band is specified by
2034 two frequencies, to the vector must have an even length.
2036 A specifies the amplitude of the desired response at each band edge.
2038 W is an optional weighting function that contains one value for each
2039 band that weights the mean squared error in that band. A must be the
2040 same length as F, and W must be half the length of F.
2042 The least squares optimization algorithm for computing FIR filter
2043 coefficients is derived in detail in:
2045 I. Selesnick, "Linear-Phase FIR Filter Design by Least Squares,"
2046 http://cnx.org/content/m10577
2050 # name: <cell-element>
2055 b = firls(N, F, A, W);
2060 # name: <cell-element>
2067 # name: <cell-element>
2071 flattopwin(L, [periodic|symmetric])
2073 Return the window f(w):
2075 f(w) = 1 - 1.93 cos(2 pi w) + 1.29 cos(4 pi w)
2076 - 0.388 cos(6 pi w) + 0.0322cos(8 pi w)
2078 where w = i/(L-1) for i=0:L-1 for a symmetric window, or
2079 w = i/L for i=0:L-1 for a periodic window. The default
2080 is symmetric. The returned window is normalized to a peak
2083 This window has low pass-band ripple, but high bandwidth.
2087 The main use for the Flat Top window is for calibration, due
2088 to its negligible amplitude errors.
2090 [1] Gade, S; Herlufsen, H; (1987) "Use of weighting functions in DFT/FFT
2091 analysis (Part I)", Bruel & Kjaer Technical Review No.3.
2095 # name: <cell-element>
2099 flattopwin(L, [periodic|symmetric])
2104 # name: <cell-element>
2111 # name: <cell-element>
2115 -- Function File: [Y H]= fracshift(X,D)
2116 -- Function File: Y = fracshift(X,D,H)
2117 Shift the series X by a (possibly fractional) number of samples D.
2118 The interpolator H is either specified or either designed with a
2119 Kaiser-windowed sinecard.
2126 # name: <cell-element>
2130 Shift the series X by a (possibly fractional) number of samples D.
2134 # name: <cell-element>
2141 # name: <cell-element>
2145 Usage: H = freqs(B,A,W);
2147 Compute the s-plane frequency response of the IIR filter B(s)/A(s) as
2148 H = polyval(B,j*W)./polyval(A,j*W). If called with no output
2149 argument, a plot of magnitude and phase are displayed.
2152 B = [1 2]; A = [1 1];
2153 w = linspace(0,4,128);
2158 # name: <cell-element>
2162 Usage: H = freqs(B,A,W);
2167 # name: <cell-element>
2174 # name: <cell-element>
2178 -- Function File: freqs_plot (W, H)
2179 Plot the amplitude and phase of the vector H.
2185 # name: <cell-element>
2189 Plot the amplitude and phase of the vector H.
2193 # name: <cell-element>
2200 # name: <cell-element>
2204 Compute peak full-width at half maximum (FWHM) or at another level of peak
2205 maximum for vector or matrix data y, optionally sampled as y(x). If y is
2206 a matrix, return FWHM for each column as a row vector.
2208 f = fwhm({x, } y {, 'zero'|'min' {, 'rlevel', rlevel}})
2209 f = fwhm({x, } y {, 'alevel', alevel})
2213 f = fwhm(x, y, 'zero')
2214 f = fwhm(x, y, 'min')
2215 f = fwhm(x, y, 'alevel', 15.3)
2216 f = fwhm(x, y, 'zero', 'rlevel', 0.5)
2217 f = fwhm(x, y, 'min', 'rlevel', 0.1)
2219 The default option 'zero' computes fwhm at half maximum, i.e. 0.5*max(y).
2220 The option 'min' computes fwhm at the middle curve, i.e. 0.5*(min(y)+max(y)).
2222 The option 'rlevel' computes full-width at the given relative level of peak
2223 profile, i.e. at rlevel*max(y) or rlevel*(min(y)+max(y)), respectively.
2224 For example, fwhm(..., 'rlevel', 0.1) computes full width at 10 % of peak
2225 maximum with respect to zero or minimum; FWHM is equivalent to
2226 fwhm(..., 'rlevel', 0.5).
2228 The option 'alevel' computes full-width at the given absolute level of y.
2230 Return 0 if FWHM does not exist (e.g. monotonous function or the function
2231 does not cut horizontal line at rlevel*max(y) or rlevel*(max(y)+min(y)) or
2232 alevel, respectively).
2234 Compatibility: Octave 3.x, Matlab
2238 # name: <cell-element>
2242 Compute peak full-width at half maximum (FWHM) or at another level of peak
2247 # name: <cell-element>
2254 # name: <cell-element>
2258 -- Function File: [Y] = gauspuls(T,FC,BW)
2259 Return the Gaussian modulated sinusoidal pulse.
2264 # name: <cell-element>
2268 Return the Gaussian modulated sinusoidal pulse.
2272 # name: <cell-element>
2279 # name: <cell-element>
2283 usage: w = gaussian(n, a)
2285 Generate an n-point gaussian convolution window of the given
2286 width. Use larger a for a narrower window. Use larger n for
2289 w = exp ( -(a*x)^2/2 )
2291 for x = linspace ( -(n-1)/2, (n-1)/2, n ).
2293 Width a is measured in frequency units (sample rate/num samples).
2294 It should be f when multiplying in the time domain, but 1/f when
2295 multiplying in the frequency domain (for use in convolutions).
2299 # name: <cell-element>
2303 usage: w = gaussian(n, a)
2308 # name: <cell-element>
2315 # name: <cell-element>
2319 usage: w = gausswin(L, a)
2321 Generate an L-point gaussian window of the given width. Use larger a
2322 for a narrow window. Use larger L for a smoother curve.
2324 w = exp ( -(a*x)^2/2 )
2326 for x = linspace(-(L-1)/L, (L-1)/L, L)
2330 # name: <cell-element>
2334 usage: w = gausswin(L, a)
2339 # name: <cell-element>
2346 # name: <cell-element>
2350 -- Function File: [Y] = gmonopuls(T,FC)
2351 Return the gaussian monopulse.
2356 # name: <cell-element>
2360 Return the gaussian monopulse.
2364 # name: <cell-element>
2371 # name: <cell-element>
2375 Compute the group delay of a filter.
2377 [g, w] = grpdelay(b)
2378 returns the group delay g of the FIR filter with coefficients b.
2379 The response is evaluated at 512 angular frequencies between 0 and
2380 pi. w is a vector containing the 512 frequencies.
2381 The group delay is in units of samples. It can be converted
2382 to seconds by multiplying by the sampling period (or dividing by
2383 the sampling rate fs).
2385 [g, w] = grpdelay(b,a)
2386 returns the group delay of the rational IIR filter whose numerator
2387 has coefficients b and denominator coefficients a.
2389 [g, w] = grpdelay(b,a,n)
2390 returns the group delay evaluated at n angular frequencies. For fastest
2391 computation n should factor into a small number of small primes.
2393 [g, w] = grpdelay(b,a,n,'whole')
2394 evaluates the group delay at n frequencies between 0 and 2*pi.
2396 [g, f] = grpdelay(b,a,n,Fs)
2397 evaluates the group delay at n frequencies between 0 and Fs/2.
2399 [g, f] = grpdelay(b,a,n,'whole',Fs)
2400 evaluates the group delay at n frequencies between 0 and Fs.
2402 [g, w] = grpdelay(b,a,w)
2403 evaluates the group delay at frequencies w (radians per sample).
2405 [g, f] = grpdelay(b,a,f,Fs)
2406 evaluates the group delay at frequencies f (in Hz).
2409 plots the group delay vs. frequency.
2411 If the denominator of the computation becomes too small, the group delay
2412 is set to zero. (The group delay approaches infinity when
2413 there are poles or zeros very close to the unit circle in the z plane.)
2415 Theory: group delay, g(w) = -d/dw [arg{H(e^jw)}], is the rate of change of
2416 phase with respect to frequency. It can be computed as:
2419 g(w) = -------------
2423 H(z) = B(z)/A(z) = sum(b_k z^k)/sum(a_k z^k).
2425 By the quotient rule,
2426 A(z) d/dw B(z) - B(z) d/dw A(z)
2427 d/dw H(z) = -------------------------------
2429 Substituting into the expression above yields:
2431 g(w) = ----------- = dB/B - dA/A
2435 d/dw B(e^-jw) = sum(k b_k e^-jwk)
2436 d/dw A(e^-jw) = sum(k a_k e^-jwk)
2437 which is just the FFT of the coefficients multiplied by a ramp.
2439 As a further optimization when nfft>>length(a), the IIR filter (b,a)
2440 is converted to the FIR filter conv(b,fliplr(conj(a))).
2441 For further details, see
2442 http://ccrma.stanford.edu/~jos/filters/Numerical_Computation_Group_Delay.html
2446 # name: <cell-element>
2450 Compute the group delay of a filter.
2454 # name: <cell-element>
2461 # name: <cell-element>
2470 # name: <cell-element>
2480 # name: <cell-element>
2487 # name: <cell-element>
2491 -- Function File: H = hilbert (F,N,DIM)
2492 Analytic extension of real valued signal
2494 `H=hilbert(F)' computes the extension of the real valued signal F
2495 to an analytic signal. If F is a matrix, the transformation is
2496 applied to each column. For N-D arrays, the transformation is
2497 applied to the first non-singleton dimension.
2499 `real(H)' contains the original signal F. `imag(H)' contains the
2500 Hilbert transform of F.
2502 `hilbert(F,N)' does the same using a length N Hilbert transform.
2503 The result will also have length N.
2505 `hilbert(F,[],DIM)' or `hilbert(F,N,DIM)' does the same along
2511 # name: <cell-element>
2515 Analytic extension of real valued signal
2520 # name: <cell-element>
2527 # name: <cell-element>
2532 Computes the inverse discrete cosine transform of x. If n is
2533 given, then x is padded or trimmed to length n before computing
2534 the transform. If x is a matrix, compute the transform along the
2535 columns of the the matrix. The transform is faster if x is
2536 real-valued and even length.
2538 The inverse discrete cosine transform x of X can be defined as follows:
2541 x[n] = sum w(k) X[k] cos (pi (2n+1) k / 2N ), n = 0, ..., N-1
2544 with w(0) = sqrt(1/N) and w(k) = sqrt(2/N), k = 1, ..., N-1
2546 See also: idct, dct2, idct2, dctmtx
2550 # name: <cell-element>
2555 Computes the inverse discrete cosine transform of x.
2559 # name: <cell-element>
2566 # name: <cell-element>
2571 Computes the inverse 2-D discrete cosine transform of matrix x
2573 y = idct2 (x, m, n) or y = idct2 (x, [m n])
2574 Computes the 2-D inverse DCT of x after padding or trimming rows to m and
2579 # name: <cell-element>
2584 Computes the inverse 2-D discrete cosine transform of matrix x
2588 # name: <cell-element>
2595 # name: <cell-element>
2599 -- Function File: Y = idst (X)
2600 -- Function File: Y = idst (X, N)
2601 Computes the inverse type I discrete sine transform of Y. If N is
2602 given, then Y is padded or trimmed to length N before computing
2603 the transform. If Y is a matrix, compute the transform along the
2604 columns of the the matrix.
2612 # name: <cell-element>
2616 Computes the inverse type I discrete sine transform of Y.
2620 # name: <cell-element>
2627 # name: <cell-element>
2631 -- Function File: m = ifht ( d, n, dim )
2632 The function ifht calculates Fast Hartley Transform where D is
2633 the real input vector (matrix), and M is the real-transform
2634 vector. For matrices the hartley transform is calculated along the
2635 columns by default. The options N, and DIM are similar to the
2636 options of FFT function.
2638 The forward and inverse hartley transforms are the same (except
2639 for a scale factor of 1/N for the inverse hartley transform), but
2640 implemented using different functions .
2642 The definition of the forward hartley transform for vector d,
2644 m[K] = 1/N \sum_i=0^N-1 d[i]*(cos[K*2*pi*i/N] + sin[K*2*pi*i/N]),
2645 for 0 <= K < N. m[K] = 1/N \sum_i=0^N-1 d[i]*CAS[K*i], for 0 <=
2656 # name: <cell-element>
2660 The function ifht calculates Fast Hartley Transform where D is the
2665 # name: <cell-element>
2672 # name: <cell-element>
2676 IIR Low Pass Filter to Multiband Filter Transformation
2678 [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt)
2679 [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt,Pass)
2681 Num,Den: numerator,denominator of the transformed filter
2682 AllpassNum,AllpassDen: numerator,denominator of allpass transform,
2683 B,A: numerator,denominator of prototype low pass filter
2684 Wo: normalized_angular_frequency/pi to be transformed
2685 Wt: [phi=normalized_angular_frequencies]/pi target vector
2686 Pass: This parameter may have values 'pass' or 'stop'. If
2687 not given, it defaults to the value of 'pass'.
2689 With normalized ang. freq. targets 0 < phi(1) < ... < phi(n) < pi radians
2691 for Pass == 'pass', the target multiband magnitude will be:
2692 -------- ---------- -----------...
2694 0 phi(1) phi(2) phi(3) phi(4) phi(5) (phi(6)) pi
2696 for Pass == 'stop', the target multiband magnitude will be:
2697 ------- --------- ----------...
2699 0 phi(1) phi(2) phi(3) phi(4) (phi(5)) pi
2702 [B, A] = butter(6, 0.5);
2703 [Num, Den] = iirlp2mb(B, A, 0.5, [.2 .4 .6 .8]);
2707 # name: <cell-element>
2711 IIR Low Pass Filter to Multiband Filter Transformation
2716 # name: <cell-element>
2723 # name: <cell-element>
2727 -- Function File: [B_OUT, A_OUT] = impinvar (B, A, FS, TOL)
2728 -- Function File: [B_OUT, A_OUT] = impinvar (B, A, FS)
2729 -- Function File: [B_OUT, A_OUT] = impinvar (B, A)
2730 Converts analog filter with coefficients B and A to digital,
2731 conserving impulse response.
2733 If FS is not specificied, or is an empty vector, it defaults to
2736 If TOL is not specified, it defaults to 0.0001 (0.1%) This
2737 function does the inverse of impinvar so that the following
2738 example should restore the original values of A and B.
2740 `invimpinvar' implements the reverse of this function.
2741 [b, a] = impinvar (b, a);
2742 [b, a] = invimpinvar (b, a);
2744 Reference: Thomas J. Cavicchi (1996) "Impulse invariance and
2745 multiple-order poles". IEEE transactions on signal processing, Vol
2748 See also: bilinear, invimpinvar
2754 # name: <cell-element>
2758 Converts analog filter with coefficients B and A to digital, conserving
2763 # name: <cell-element>
2770 # name: <cell-element>
2774 usage: [x, t] = impz(b [, a, n, fs])
2776 Generate impulse-response characteristics of the filter. The filter
2777 coefficients correspond to the the z-plane rational function with
2778 numerator b and denominator a. If a is not specified, it defaults to
2779 1. If n is not specified, or specified as [], it will be chosen such
2780 that the signal has a chance to die down to -120dB, or to not explode
2781 beyond 120dB, or to show five periods if there is no significant
2782 damping. If no return arguments are requested, plot the results.
2784 See also: freqz, zplane
2788 # name: <cell-element>
2792 usage: [x, t] = impz(b [, a, n, fs])
2797 # name: <cell-element>
2804 # name: <cell-element>
2808 usage: y = interp(x, q [, n [, Wc]])
2810 Upsample the signal x by a factor of q, using an order 2*q*n+1 FIR
2811 filter. Note that q must be an integer for this rate change method.
2812 n defaults to 4 and Wc defaults to 0.5.
2815 # Generate a signal.
2816 t=0:0.01:2; x=chirp(t,2,.5,10,'quadratic')+sin(2*pi*t*0.4);
2817 y = interp(x(1:4:length(x)),4,4,1); # interpolate a sub-sample
2818 stem(t(1:121)*1000,x(1:121),"-g;Original;"); hold on;
2819 stem(t(1:121)*1000,y(1:121),"-r;Interpolated;");
2820 stem(t(1:4:121)*1000,x(1:4:121),"-b;Subsampled;"); hold off;
2822 See also: decimate, resample
2826 # name: <cell-element>
2830 usage: y = interp(x, q [, n [, Wc]])
2835 # name: <cell-element>
2842 # name: <cell-element>
2846 usage: [B,A] = invfreq(H,F,nB,nA)
2847 [B,A] = invfreq(H,F,nB,nA,W)
2848 [B,A] = invfreq(H,F,nB,nA,W,[],[],plane)
2849 [B,A] = invfreq(H,F,nB,nA,W,iter,tol,plane)
2851 Fit filter B(z)/A(z) or B(s)/A(s) to complex frequency response at
2852 frequency points F. A and B are real polynomial coefficients of order
2853 nA and nB respectively. Optionally, the fit-errors can be weighted vs
2854 frequency according to the weights W. Also, the transform plane can be
2855 specified as either 's' for continuous time or 'z' for discrete time. 'z'
2856 is chosen by default. Eventually, Steiglitz-McBride iterations will be
2857 specified by iter and tol.
2859 H: desired complex frequency response
2860 It is assumed that A and B are real polynomials, hence H is one-sided.
2861 F: vector of frequency samples in radians
2862 nA: order of denominator polynomial A
2863 nB: order of numerator polynomial B
2864 plane='z': F on unit circle (discrete-time spectra, z-plane design)
2865 plane='s': F on jw axis (continuous-time spectra, s-plane design)
2866 H(k) = spectral samples of filter frequency response at points zk,
2867 where zk=exp(sqrt(-1)*F(k)) when plane='z' (F(k) in [0,.5])
2868 and zk=(sqrt(-1)*F(k)) when plane='s' (F(k) nonnegative)
2870 [B,A] = butter(12,1/4);
2871 [H,w] = freqz(B,A,128);
2872 [Bh,Ah] = invfreq(H,F,4,4);
2874 disp(sprintf('||frequency response error|| = %f',norm(H-Hh)));
2876 References: J. O. Smith, "Techniques for Digital Filter Design and System
2877 Identification with Application to the Violin, Ph.D. Dissertation,
2878 Elec. Eng. Dept., Stanford University, June 1983, page 50; or,
2880 http://ccrma.stanford.edu/~jos/filters/FFT_Based_Equation_Error_Method.html
2884 # name: <cell-element>
2888 usage: [B,A] = invfreq(H,F,nB,nA)
2889 [B,A] = invfreq(H,F,nB,nA,W)
2894 # name: <cell-element>
2901 # name: <cell-element>
2905 Usage: [B,A] = invfreqs(H,F,nB,nA)
2906 [B,A] = invfreqs(H,F,nB,nA,W)
2907 [B,A] = invfreqs(H,F,nB,nA,W,iter,tol,'trace')
2909 Fit filter B(s)/A(s)to the complex frequency response H at frequency
2910 points F. A and B are real polynomial coefficients of order nA and nB.
2911 Optionally, the fit-errors can be weighted vs frequency according to
2913 Note: all the guts are in invfreq.m
2915 H: desired complex frequency response
2916 F: frequency (must be same length as H)
2917 nA: order of the denominator polynomial A
2918 nB: order of the numerator polynomial B
2919 W: vector of weights (must be same length as F)
2924 w = linspace(0,4,128);
2926 [Bh,Ah] = invfreqs(H,w,1,1);
2927 Hh = freqs(Bh,Ah,w);
2928 plot(w,[abs(H);abs(Hh)])
2929 legend('Original','Measured');
2931 disp(sprintf('L2 norm of frequency response error = %f',err));
2935 # name: <cell-element>
2939 Usage: [B,A] = invfreqs(H,F,nB,nA)
2940 [B,A] = invfreqs(H,F,nB,nA,W)
2945 # name: <cell-element>
2952 # name: <cell-element>
2956 usage: [B,A] = invfreqz(H,F,nB,nA)
2957 [B,A] = invfreqz(H,F,nB,nA,W)
2958 [B,A] = invfreqz(H,F,nB,nA,W,iter,tol,'trace')
2960 Fit filter B(z)/A(z)to the complex frequency response H at frequency
2961 points F. A and B are real polynomial coefficients of order nA and nB.
2962 Optionally, the fit-errors can be weighted vs frequency according to
2964 Note: all the guts are in invfreq.m
2966 H: desired complex frequency response
2967 F: normalized frequncy (0 to pi) (must be same length as H)
2968 nA: order of the denominator polynomial A
2969 nB: order of the numerator polynomial B
2970 W: vector of weights (must be same length as F)
2973 [B,A] = butter(4,1/4);
2975 [Bh,Ah] = invfreq(H,F,4,4);
2977 disp(sprintf('||frequency response error|| = %f',norm(H-Hh)));
2981 # name: <cell-element>
2985 usage: [B,A] = invfreqz(H,F,nB,nA)
2986 [B,A] = invfreqz(H,F,nB,nA,W)
2991 # name: <cell-element>
2998 # name: <cell-element>
3002 -- Function File: [B_OUT, A_OUT] = invimpinvar (B, A, FS, TOL)
3003 -- Function File: [B_OUT, A_OUT] = invimpinvar (B, A, FS)
3004 -- Function File: [B_OUT, A_OUT] = invimpinvar (B, A)
3005 Converts digital filter with coefficients B and A to analog,
3006 conserving impulse response.
3008 This function does the inverse of impinvar so that the following
3009 example should restore the original values of A and B.
3010 [b, a] = impinvar (b, a);
3011 [b, a] = invimpinvar (b, a);
3013 If FS is not specificied, or is an empty vector, it defaults to
3016 If TOL is not specified, it defaults to 0.0001 (0.1%)
3018 Reference: Thomas J. Cavicchi (1996) "Impulse invariance and
3019 multiple-order poles". IEEE transactions on signal processing, Vol
3022 See also: bilinear, impinvar
3028 # name: <cell-element>
3032 Converts digital filter with coefficients B and A to analog, conserving
3037 # name: <cell-element>
3044 # name: <cell-element>
3048 usage: kaiser (L, beta)
3050 Returns the filter coefficients of the L-point Kaiser window with
3053 For the definition of the Kaiser window, see A. V. Oppenheim &
3054 R. W. Schafer, "Discrete-Time Signal Processing".
3056 The continuous version of width L centered about x=0 is:
3058 besseli(0, beta * sqrt(1-(2*x/L).^2))
3059 k(x) = -------------------------------------, L/2 <= x <= L/2
3066 # name: <cell-element>
3070 usage: kaiser (L, beta)
3075 # name: <cell-element>
3082 # name: <cell-element>
3086 usage: [n, Wn, beta, ftype] = kaiserord(f, m, dev [, fs])
3088 Returns the parameters needed for fir1 to produce a filter of the
3089 desired specification from a kaiser window:
3090 n: order of the filter (length of filter minus 1)
3091 Wn: band edges for use in fir1
3092 beta: parameter for kaiser window of length n+1
3093 ftype: choose between pass and stop bands
3094 b = fir1(n,Wn,kaiser(n+1,beta),ftype,'noscale');
3096 f: frequency bands, given as pairs, with the first half of the
3097 first pair assumed to start at 0 and the last half of the last
3098 pair assumed to end at 1. It is important to separate the
3099 band edges, since narrow transition regions require large order
3101 m: magnitude within each band. Should be non-zero for pass band
3102 and zero for stop band. All passbands must have the same
3103 magnitude, or you will get the error that pass and stop bands
3104 must be strictly alternating.
3105 dev: deviation within each band. Since all bands in the resulting
3106 filter have the same deviation, only the minimum deviation is
3107 used. In this version, a single scalar will work just as well.
3108 fs: sampling rate. Used to convert the frequency specification into
3109 the [0, 1], where 1 corresponds to the Nyquist frequency, fs/2.
3111 The Kaiser window parameters n and beta are computed from the
3112 relation between ripple (A=-20*log10(dev)) and transition width
3113 (dw in radians) discovered empirically by Kaiser:
3115 / 0.1102(A-8.7) A > 50
3116 beta = | 0.5842(A-21)^0.4 + 0.07886(A-21) 21 <= A <= 50
3119 n = (A-8)/(2.285 dw)
3122 [n, w, beta, ftype] = kaiserord([1000,1200], [1,0], [0.05,0.05], 11025);
3123 freqz(fir1(n,w,kaiser(n+1,beta),ftype,'noscale'),1,[],11025);
3127 # name: <cell-element>
3131 usage: [n, Wn, beta, ftype] = kaiserord(f, m, dev [, fs])
3136 # name: <cell-element>
3143 # name: <cell-element>
3147 usage: [a, v, ref] = levinson (acf [, p])
3149 Use the Durbin-Levinson algorithm to solve:
3150 toeplitz(acf(1:p)) * x = -acf(2:p+1).
3151 The solution [1, x'] is the denominator of an all pole filter
3152 approximation to the signal x which generated the autocorrelation
3155 acf is the autocorrelation function for lags 0 to p.
3156 p defaults to length(acf)-1.
3158 a=[1, x'] the denominator filter coefficients.
3159 v= variance of the white noise = square of the numerator constant
3160 ref = reflection coefficients = coefficients of the lattice
3161 implementation of the filter
3162 Use freqz(sqrt(v),a) to plot the power spectrum.
3165 [1] Steven M. Kay and Stanley Lawrence Marple Jr.:
3166 "Spectrum analysis -- a modern perspective",
3167 Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
3171 # name: <cell-element>
3175 usage: [a, v, ref] = levinson (acf [, p])
3180 # name: <cell-element>
3187 # name: <cell-element>
3191 -- Function File: Q = marcumq (A, B)
3192 -- Function File: Q = marcumq (A, B, M)
3193 -- Function File: Q = marcumq (A, B, M, TOL)
3194 Compute the generalized Marcum Q function of order M with
3195 noncentrality parameter A and argument B. If the order M is
3196 omitted it defaults to 1. An optional relative tolerance TOL may
3197 be included, the default is `eps'.
3199 If the input arguments are commensurate vectors, this function
3200 will produce a table of values.
3202 This function computes Marcum's Q function using the infinite
3203 Bessel series, truncated when the relative error is less than the
3204 specified tolerance. The accuracy is limited by that of the
3205 Bessel functions, so reducing the tolerance is probably not useful.
3207 Reference: Marcum, "Tables of Q Functions", Rand Corporation.
3209 Reference: R.T. Short, "Computation of Noncentral Chi-squared and
3210 Rice Random Variables", www.phaselockedsystems.com/publications
3216 # name: <cell-element>
3220 Compute the generalized Marcum Q function of order M with noncentrality
3225 # name: <cell-element>
3232 # name: <cell-element>
3236 -- Function File: [PSI,X] = mexihat(LB,UB,N)
3237 Compute the Mexican hat wavelet.
3242 # name: <cell-element>
3246 Compute the Mexican hat wavelet.
3250 # name: <cell-element>
3257 # name: <cell-element>
3261 -- Function File: [Y] = meyeraux(X)
3262 Compute the Meyer wavelet auxiliary function.
3267 # name: <cell-element>
3271 Compute the Meyer wavelet auxiliary function.
3275 # name: <cell-element>
3282 # name: <cell-element>
3286 -- Function File: [PSI,X] = morlet(LB,UB,N)
3287 Compute the Morlet wavelet.
3292 # name: <cell-element>
3296 Compute the Morlet wavelet.
3300 # name: <cell-element>
3307 # name: <cell-element>
3312 [Pxx,freq]=mscohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
3314 Estimate (mean square) coherence of signals "x" and "y".
3315 Use the Welch (1967) periodogram/FFT method.
3316 See "help pwelch" for description of arguments, hints and references
3320 # name: <cell-element>
3325 [Pxx,freq]=mscohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detren
3329 # name: <cell-element>
3336 # name: <cell-element>
3340 usage: [Zz, Zp, Zg] = ncauer(Rp, Rs, n)
3342 Analog prototype for Cauer filter.
3343 [z, p, g]=ncauer(Rp, Rs, ws)
3344 Rp = Passband ripple
3345 Rs = Stopband ripple
3350 - Serra, Celso Penteado, Teoria e Projeto de Filtros, Campinas: CARTGRAF,
3352 - Lamar, Marcus Vinicius, Notas de aula da disciplina TE 456 - Circuitos
3353 Analogicos II, UFPR, 2001/2002.
3357 # name: <cell-element>
3361 usage: [Zz, Zp, Zg] = ncauer(Rp, Rs, n)
3366 # name: <cell-element>
3373 # name: <cell-element>
3377 -- Function File: [W] = nuttallwin(L)
3378 Compute the Blackman-Harris window defined by Nuttall of length L.
3380 See also: blackman, blackmanharris
3386 # name: <cell-element>
3390 Compute the Blackman-Harris window defined by Nuttall of length L.
3394 # name: <cell-element>
3401 # name: <cell-element>
3405 -- Function File: [W] = parzenwin(L)
3406 Compute the Parzen window of lenght L.
3408 See also: rectwin, bartlett
3414 # name: <cell-element>
3418 Compute the Parzen window of lenght L.
3422 # name: <cell-element>
3429 # name: <cell-element>
3434 [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion)
3436 Calculate Burg maximum-entropy power spectral density.
3437 The functions "arburg" and "ar_psd" do all the work.
3438 See "help arburg" and "help ar_psd" for further details.
3441 All but the first two arguments are optional and may be empty.
3442 x %% [vector] sampled data
3444 poles %% [integer scalar] required number of poles of the AR model
3446 freq %% [real vector] frequencies at which power spectral density
3448 %% [integer scalar] number of uniformly distributed frequency
3449 %% values at which spectral density is calculated.
3452 Fs %% [real scalar] sampling frequency (Hertz) [default=1]
3455 CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
3456 Control-string arguments can be in any order after the other arguments.
3459 range %% 'half', 'onesided' : frequency range of the spectrum is
3460 %% from zero up to but not including sample_f/2. Power
3461 %% from negative frequencies is added to the positive
3462 %% side of the spectrum.
3463 %% 'whole', 'twosided' : frequency range of the spectrum is
3464 %% -sample_f/2 to sample_f/2, with negative frequencies
3465 %% stored in "wrap around" order after the positive
3466 %% frequencies; e.g. frequencies for a 10-point 'twosided'
3467 %% spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
3468 %% 'shift', 'centerdc' : same as 'whole' but with the first half
3469 %% of the spectrum swapped with second half to put the
3470 %% zero-frequency value in the middle. (See "help
3471 %% fftshift". If "freq" is vector, 'shift' is ignored.
3472 %% If model coefficients "ar_coeffs" are real, the default
3473 %% range is 'half', otherwise default range is 'whole'.
3475 method %% 'fft': use FFT to calculate power spectral density.
3476 %% 'poly': calculate spectral density as a polynomial of 1/z
3477 %% N.B. this argument is ignored if the "freq" argument is a
3478 %% vector. The default is 'poly' unless the "freq"
3479 %% argument is an integer power of 2.
3481 plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
3482 %% specifies the type of plot. The default is 'plot', which
3483 %% means linear-linear axes. 'squared' is the same as 'plot'.
3484 %% 'dB' plots "10*log10(psd)". This argument is ignored and a
3485 %% spectrum is not plotted if the caller requires a returned
3488 criterion %% [optional string arg] model-selection criterion. Limits
3489 %% the number of poles so that spurious poles are not
3490 %% added when the whitened data has no more information
3491 %% in it (see Kay & Marple, 1981). Recognised values are
3492 %% 'AKICc' -- approximate corrected Kullback information
3493 %% criterion (recommended),
3494 %% 'KIC' -- Kullback information criterion
3495 %% 'AICc' -- corrected Akaike information criterion
3496 %% 'AIC' -- Akaike information criterion
3497 %% 'FPE' -- final prediction error" criterion
3498 %% The default is to NOT use a model-selection criterion
3501 If return values are not required by the caller, the spectrum
3502 is plotted and nothing is returned.
3503 psd %% [real vector] power-spectral density estimate
3504 f_out %% [real vector] frequency values
3507 This function is a wrapper for arburg and ar_psd.
3508 See "help arburg", "help ar_psd".
3512 # name: <cell-element>
3517 [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion
3521 # name: <cell-element>
3528 # name: <cell-element>
3532 -- Function File: [ B, A ] = pei_tseng_notch ( FREQUENCIES,
3534 Return coefficients for an IIR notch-filter with one or more
3535 filter frequencies and according (very narrow) bandwidths to be
3536 used with `filter' or `filtfilt'. The filter construction is
3537 based on an allpass which performs a reversal of phase at the
3538 filter frequencies. Thus, the mean of the phase-distorted and the
3539 original signal has the respective frequencies removed. See the
3540 demo for an illustration.
3542 Original source: Pei, Soo-Chang, and Chien-Cheng Tseng "IIR
3543 Multiple Notch Filter Design Based on Allpass Filter" 1996 IEEE
3544 Tencon doi: 10.1109/TENCON.1996.608814)
3549 # name: <cell-element>
3553 Return coefficients for an IIR notch-filter with one or more filter
3558 # name: <cell-element>
3565 # name: <cell-element>
3571 Stabalize the polynomial transfer function by replacing all roots
3572 outside the unit circle with their reflection inside the unit circle.
3576 # name: <cell-element>
3585 # name: <cell-element>
3592 # name: <cell-element>
3596 usage: y=pulstran(t,d,'func',...)
3597 y=pulstran(t,d,p,Fs,'interp')
3599 Generate the signal y=sum(func(t+d,...)) for each d. If d is a
3600 matrix of two columns, the first column is the delay d and the second
3601 column is the amplitude a, and y=sum(a*func(t+d)) for each d,a.
3602 Clearly, func must be a function which accepts a vector of times.
3603 Any extra arguments needed for the function must be tagged on the end.
3606 fs = 11025; # arbitrary sample rate
3607 f0 = 100; # pulse train sample rate
3608 w = 0.001; # pulse width of 1 millisecond
3609 auplot(pulstran(0:1/fs:0.1, 0:1/f0:0.1, 'rectpuls', w), fs);
3611 If instead of a function name you supply a pulse shape sampled at
3612 frequency Fs (default 1 Hz), an interpolated version of the pulse
3613 is added at each delay d. The interpolation stays within the the
3614 time range of the delayed pulse. The interpolation method defaults
3615 to linear, but it can be any interpolation method accepted by the
3619 fs = 11025; # arbitrary sample rate
3620 f0 = 100; # pulse train sample rate
3621 w = boxcar(10); # pulse width of 1 millisecond at 10 kHz
3622 auplot(pulstran(0:1/fs:0.1, 0:1/f0:0.1, w, 10000), fs);
3626 # name: <cell-element>
3630 usage: y=pulstran(t,d,'func',.
3634 # name: <cell-element>
3641 # name: <cell-element>
3646 [spectra,freq] = pwelch(x,window,overlap,Nfft,Fs,
3647 range,plot_type,detrend,sloppy)
3648 Estimate power spectral density of data "x" by the Welch (1967)
3649 periodogram/FFT method. All arguments except "x" are optional.
3650 The data is divided into segments. If "window" is a vector, each
3651 segment has the same length as "window" and is multiplied by "window"
3652 before (optional) zero-padding and calculation of its periodogram. If
3653 "window" is a scalar, each segment has a length of "window" and a
3654 Hamming window is used.
3655 The spectral density is the mean of the periodograms, scaled so that
3656 area under the spectrum is the same as the mean square of the
3657 data. This equivalence is supposed to be exact, but in practice there
3658 is a mismatch of up to 0.5% when comparing area under a periodogram
3659 with the mean square of the data.
3661 [spectra,freq] = pwelch(x,y,window,overlap,Nfft,Fs,
3662 range,plot_type,detrend,sloppy,results)
3663 Two-channel spectrum analyser. Estimate power spectral density, cross-
3664 spectral density, transfer function and/or coherence functions of time-
3665 series input data "x" and output data "y" by the Welch (1967)
3666 periodogram/FFT method.
3667 pwelch treats the second argument as "y" if there is a control-string
3668 argument "cross", "trans", "coher" or "ypower"; "power" does not force
3669 the 2nd argument to be treated as "y". All other arguments are
3670 optional. All spectra are returned in matrix "spectra".
3672 [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf,
3673 range,plot_type,detrend,sloppy)
3674 [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf,
3675 range,plot_type,detrend,sloppy,results)
3676 Estimates confidence intervals for the spectral density.
3677 See Hint (7) below for compatibility options. Confidence level "conf"
3678 is the 6th or 7th numeric argument. If "results" control-string
3679 arguments are used, one of them must be "power" when the "conf"
3680 argument is present; pwelch can estimate confidence intervals only for
3681 the power spectrum of the "x" data. It does not know how to estimate
3682 confidence intervals of the cross-power spectrum, transfer function or
3683 coherence; if you can suggest a good method, please send a bug report.
3686 All but the first argument are optional and may be empty, except that
3687 the "results" argument may require the second argument to be "y".
3689 x %% [non-empty vector] system-input time-series data
3690 y %% [non-empty vector] system-output time-series data
3692 window %% [real vector] of window-function values between 0 and 1; the
3693 %% data segment has the same length as the window.
3694 %% Default window shape is Hamming.
3695 %% [integer scalar] length of each data segment. The default
3696 %% value is window=sqrt(length(x)) rounded up to the
3697 %% nearest integer power of 2; see 'sloppy' argument.
3699 overlap %% [real scalar] segment overlap expressed as a multiple of
3700 %% window or segment length. 0 <= overlap < 1,
3701 %% The default is overlap=0.5 .
3703 Nfft %% [integer scalar] Length of FFT. The default is the length
3704 %% of the "window" vector or has the same value as the
3705 %% scalar "window" argument. If Nfft is larger than the
3706 %% segment length, "seg_len", the data segment is padded
3707 %% with "Nfft-seg_len" zeros. The default is no padding.
3708 %% Nfft values smaller than the length of the data
3709 %% segment (or window) are ignored silently.
3711 Fs %% [real scalar] sampling frequency (Hertz); default=1.0
3713 conf %% [real scalar] confidence level between 0 and 1. Confidence
3714 %% intervals of the spectral density are estimated from
3715 %% scatter in the periodograms and are returned as Pxx_ci.
3716 %% Pxx_ci(:,1) is the lower bound of the confidence
3717 %% interval and Pxx_ci(:,2) is the upper bound. If there
3718 %% are three return values, or conf is an empty matrix,
3719 %% confidence intervals are calculated for conf=0.95 .
3720 %% If conf is zero or is not given, confidence intervals
3721 %% are not calculated. Confidence intervals can be
3722 %% obtained only for the power spectral density of x;
3725 CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
3726 Control-string arguments must be after the other arguments but can be in
3729 range %% 'half', 'onesided' : frequency range of the spectrum is
3730 %% zero up to but not including Fs/2. Power from
3731 %% negative frequencies is added to the positive side of
3732 %% the spectrum, but not at zero or Nyquist (Fs/2)
3733 %% frequencies. This keeps power equal in time and
3734 %% spectral domains. See reference [2].
3735 %% 'whole', 'twosided' : frequency range of the spectrum is
3736 %% -Fs/2 to Fs/2, with negative frequencies
3737 %% stored in "wrap around" order after the positive
3738 %% frequencies; e.g. frequencies for a 10-point 'twosided'
3739 %% spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
3740 %% 'shift', 'centerdc' : same as 'whole' but with the first half
3741 %% of the spectrum swapped with second half to put the
3742 %% zero-frequency value in the middle. (See "help
3744 %% If data (x and y) are real, the default range is 'half',
3745 %% otherwise default range is 'whole'.
3747 plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
3748 %% specifies the type of plot. The default is 'plot', which
3749 %% means linear-linear axes. 'squared' is the same as 'plot'.
3750 %% 'dB' plots "10*log10(psd)". This argument is ignored and a
3751 %% spectrum is not plotted if the caller requires a returned
3754 detrend %% 'no-strip', 'none' -- do NOT remove mean value from the data
3755 %% 'short', 'mean' -- remove the mean value of each segment from
3756 %% each segment of the data.
3757 %% 'linear', -- remove linear trend from each segment of
3759 %% 'long-mean' -- remove the mean value from the data before
3760 %% splitting it into segments. This is the default.
3762 sloppy %% 'sloppy': FFT length is rounded up to the nearest integer
3763 %% power of 2 by zero padding. FFT length is adjusted
3764 %% after addition of padding by explicit Nfft argument.
3765 %% The default is to use exactly the FFT and window/
3769 # name: <cell-element>
3774 [spectra,freq] = pwelch(x,window,overlap,Nfft,Fs,
3779 # name: <cell-element>
3786 # name: <cell-element>
3791 [psd,f_out] = pyulear(x,poles,freq,Fs,range,method,plot_type)
3793 Calculates a Yule-Walker autoregressive (all-pole) model of the data "x"
3794 and computes the power spectrum of the model. This is a wrapper for
3795 functions "aryule" and "ar_psd" which perform the argument checking.
3796 See "help aryule" and "help ar_psd" for further details.
3799 All but the first two arguments are optional and may be empty.
3800 x %% [vector] sampled data
3802 poles %% [integer scalar] required number of poles of the AR model
3804 freq %% [real vector] frequencies at which power spectral density
3806 %% [integer scalar] number of uniformly distributed frequency
3807 %% values at which spectral density is calculated.
3810 Fs %% [real scalar] sampling frequency (Hertz) [default=1]
3813 CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
3814 Control-string arguments can be in any order after the other arguments.
3817 range %% 'half', 'onesided' : frequency range of the spectrum is
3818 %% from zero up to but not including sample_f/2. Power
3819 %% from negative frequencies is added to the positive
3820 %% side of the spectrum.
3821 %% 'whole', 'twosided' : frequency range of the spectrum is
3822 %% -sample_f/2 to sample_f/2, with negative frequencies
3823 %% stored in "wrap around" order after the positive
3824 %% frequencies; e.g. frequencies for a 10-point 'twosided'
3825 %% spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
3826 %% 'shift', 'centerdc' : same as 'whole' but with the first half
3827 %% of the spectrum swapped with second half to put the
3828 %% zero-frequency value in the middle. (See "help
3829 %% fftshift". If "freq" is vector, 'shift' is ignored.
3830 %% If model coefficients "ar_coeffs" are real, the default
3831 %% range is 'half', otherwise default range is 'whole'.
3833 method %% 'fft': use FFT to calculate power spectrum.
3834 %% 'poly': calculate power spectrum as a polynomial of 1/z
3835 %% N.B. this argument is ignored if the "freq" argument is a
3836 %% vector. The default is 'poly' unless the "freq"
3837 %% argument is an integer power of 2.
3839 plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
3840 %% specifies the type of plot. The default is 'plot', which
3841 %% means linear-linear axes. 'squared' is the same as 'plot'.
3842 %% 'dB' plots "10*log10(psd)". This argument is ignored and a
3843 %% spectrum is not plotted if the caller requires a returned
3847 If return values are not required by the caller, the spectrum
3848 is plotted and nothing is returned.
3849 psd %% [real vector] power-spectrum estimate
3850 f_out %% [real vector] frequency values
3853 This function is a wrapper for aryule and ar_psd.
3854 See "help aryule", "help ar_psd".
3858 # name: <cell-element>
3863 [psd,f_out] = pyulear(x,poles,freq,Fs,range,method,plot_type)
3868 # name: <cell-element>
3875 # name: <cell-element>
3879 Usage: qp_kaiser (nb, at, linear)
3881 Computes a finite impulse response (FIR) filter for use with a
3882 quasi-perfect reconstruction polyphase-network filter bank. This
3883 version utilizes a Kaiser window to shape the frequency response of
3884 the designed filter. Tha number nb of bands and the desired
3885 attenuation at in the stop-band are given as parameters.
3887 The Kaiser window is multiplied by the ideal impulse response
3888 h(n)=a.sinc(a.n) and converted to its minimum-phase version by means
3889 of a Hilbert transform.
3891 By using a third non-null argument, the minimum-phase calculation is
3896 # name: <cell-element>
3900 Usage: qp_kaiser (nb, at, linear)
3905 # name: <cell-element>
3912 # name: <cell-element>
3916 usage: [y, xm] = rceps(x)
3917 Produce the cepstrum of the signal x, and if desired, the minimum
3918 phase reconstruction of the signal x. If x is a matrix, do so
3919 for each column of the matrix.
3922 f0=70; Fs=10000; # 100 Hz fundamental, 10kHz sampling rate
3923 a=poly(0.985*exp(1i*pi*[0.1, -0.1, 0.3, -0.3])); # two formants
3924 s=0.005*randn(1024,1); # Noise excitation signal
3925 s(1:Fs/f0:length(s)) = 1; # Impulse glottal wave
3926 x=filter(1,a,s); # Speech signal in x
3927 [y, xm] = rceps(x.*hanning(1024)); # cepstrum and min phase reconstruction
3930 Programs for digital signal processing. IEEE Press.
3931 New York: John Wiley & Sons. 1979.
3935 # name: <cell-element>
3939 usage: [y, xm] = rceps(x)
3940 Produce the cepstrum of the signal x, and if desir
3944 # name: <cell-element>
3951 # name: <cell-element>
3955 usage: y = rectpuls(t, w)
3957 Generate a rectangular pulse over the interval [-w/2,w/2), sampled at
3958 times t. This is useful with the function pulstran for generating a
3962 fs = 11025; # arbitrary sample rate
3963 f0 = 100; # pulse train sample rate
3964 w = 0.3/f0; # pulse width 3/10th the distance between pulses
3965 auplot(pulstran(0:1/fs:4/f0, 0:1/f0:4/f0, 'rectpuls', w), fs);
3971 # name: <cell-element>
3975 usage: y = rectpuls(t, w)
3980 # name: <cell-element>
3987 # name: <cell-element>
3991 -- Function File: [W] = rectwin(L)
3992 Return the filter coefficients of a rectangle window of length L.
3994 See also: hamming, hanning
4000 # name: <cell-element>
4004 Return the filter coefficients of a rectangle window of length L.
4008 # name: <cell-element>
4015 # name: <cell-element>
4019 -- Function File: [Y H]= resample(X,P,Q)
4020 -- Function File: Y = resample(X,P,Q,H)
4021 Change the sample rate of X by a factor of P/Q. This is performed
4022 using a polyphase algorithm. The impulse response H of the
4023 antialiasing filter is either specified or either designed with a
4024 Kaiser-windowed sinecard.
4026 Ref [1] J. G. Proakis and D. G. Manolakis, Digital Signal
4027 Processing: Principles, Algorithms, and Applications, 4th ed.,
4028 Prentice Hall, 2007. Chap. 6
4030 Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
4031 Discrete-time signal processing, Signal processing series,
4037 # name: <cell-element>
4041 Change the sample rate of X by a factor of P/Q.
4045 # name: <cell-element>
4052 # name: <cell-element>
4056 -- Function File: [R, P, F, M] = residued (B, A)
4057 Compute the partial fraction expansion (PFE) of filter H(z) =
4058 B(z)/A(z). In the usual PFE function `residuez', the IIR part
4059 (poles P and residues R) is driven _in parallel_ with the FIR part
4060 (F). In this variant (`residued') the IIR part is driven by the
4061 _output_ of the FIR part. This structure can be more accurate in
4062 signal modeling applications.
4064 INPUTS: B and A are vectors specifying the digital filter H(z) =
4065 B(z)/A(z). Say `help filter' for documentation of the B and A
4066 filter coefficients.
4069 * R = column vector containing the filter-pole residues
4070 * P = column vector containing the filter poles
4071 * F = row vector containing the FIR part, if any
4072 * M = column vector of pole multiplicities
4075 Say `test residued verbose' to see a number of examples.
4077 For the theory of operation, see
4078 <http://ccrma.stanford.edu/~jos/filters/residued.html>
4080 See also: residue residued
4086 # name: <cell-element>
4090 Compute the partial fraction expansion (PFE) of filter H(z) = B(z)/A(z).
4094 # name: <cell-element>
4101 # name: <cell-element>
4105 -- Function File: [R, P, F, M] = residuez (B, A)
4106 Compute the partial fraction expansion of filter H(z) = B(z)/A(z).
4108 INPUTS: B and A are vectors specifying the digital filter H(z) =
4109 B(z)/A(z). Say `help filter' for documentation of the B and A
4110 filter coefficients.
4113 * R = column vector containing the filter-pole residues
4114 * P = column vector containing the filter poles
4115 * F = row vector containing the FIR part, if any
4116 * M = column vector of pole multiplicities
4119 Say `test residuez verbose' to see a number of examples.
4121 For the theory of operation, see
4122 <http://ccrma.stanford.edu/~jos/filters/residuez.html>
4124 See also: residue residued
4130 # name: <cell-element>
4134 Compute the partial fraction expansion of filter H(z) = B(z)/A(z).
4138 # name: <cell-element>
4145 # name: <cell-element>
4151 xt = sampled2continuous( xn , T, t )
4153 Calculate the x(t) reconstructed
4154 from samples x[n] sampled at a rate 1/T samples
4157 t is all the instants of time when you need x(t)
4158 from x[n]; this time is relative to x[0] and not
4161 This function can be used to calculate sampling rate
4162 effects on aliasing, actual signal reconstruction
4163 from discrete samples.
4167 # name: <cell-element>
4173 xt = sampled2continuous( xn , T, t )
4175 Calculate the x(t) reconstruc
4179 # name: <cell-element>
4186 # name: <cell-element>
4190 -- Function File: [Y] = sawtooth(T)
4191 -- Function File: [Y] = sawtooth(T,WIDTH)
4192 Generates a sawtooth wave of period `2 * pi' with limits `+1/-1'
4193 for the elements of T.
4195 WIDTH is a real number between `0' and `1' which specifies the
4196 point between `0' and `2 * pi' where the maximum is. The function
4197 increases linearly from `-1' to `1' in `[0, 2 * pi * WIDTH]'
4198 interval, and decreases linearly from `1' to `-1' in the interval
4199 `[2 * pi * WIDTH, 2 * pi]'.
4201 If WIDTH is 0.5, the function generates a standard triangular wave.
4203 If WIDTH is not specified, it takes a value of 1, which is a
4204 standard sawtooth function.
4209 # name: <cell-element>
4213 Generates a sawtooth wave of period `2 * pi' with limits `+1/-1' for
4218 # name: <cell-element>
4225 # name: <cell-element>
4229 usage: [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop)
4231 Transform band edges of a generic lowpass filter (cutoff at W=1)
4232 represented in splane zero-pole-gain form. W is the edge of the
4233 target filter (or edges if band pass or band stop). Stop is true for
4234 high pass and band stop filters or false for low pass and band pass
4235 filters. Filter edges are specified in radians, from 0 to pi (the
4238 Theory: Given a low pass filter represented by poles and zeros in the
4239 splane, you can convert it to a low pass, high pass, band pass or
4240 band stop by transforming each of the poles and zeros individually.
4241 The following table summarizes the transformation:
4243 Transform Zero at x Pole at x
4244 ---------------- ------------------------- ------------------------
4245 Low Pass zero: Fc x/C pole: Fc x/C
4246 S -> C S/Fc gain: C/Fc gain: Fc/C
4247 ---------------- ------------------------- ------------------------
4248 High Pass zero: Fc C/x pole: Fc C/x
4249 S -> C Fc/S pole: 0 zero: 0
4251 ---------------- ------------------------- ------------------------
4252 Band Pass zero: b ± sqrt(b^2-FhFl) pole: b ± sqrt(b^2-FhFl)
4253 S^2+FhFl pole: 0 zero: 0
4254 S -> C -------- gain: C/(Fh-Fl) gain: (Fh-Fl)/C
4255 S(Fh-Fl) b=x/C (Fh-Fl)/2 b=x/C (Fh-Fl)/2
4256 ---------------- ------------------------- ------------------------
4257 Band Stop zero: b ± sqrt(b^2-FhFl) pole: b ± sqrt(b^2-FhFl)
4258 S(Fh-Fl) pole: ±sqrt(-FhFl) zero: ±sqrt(-FhFl)
4259 S -> C -------- gain: -x gain: -1/x
4260 S^2+FhFl b=C/x (Fh-Fl)/2 b=C/x (Fh-Fl)/2
4261 ---------------- ------------------------- ------------------------
4262 Bilinear zero: (2+xT)/(2-xT) pole: (2+xT)/(2-xT)
4263 2 z-1 pole: -1 zero: -1
4264 S -> - --- gain: (2-xT)/T gain: (2-xT)/T
4266 ---------------- ------------------------- ------------------------
4268 where C is the cutoff frequency of the initial lowpass filter, Fc is
4269 the edge of the target low/high pass filter and [Fl,Fh] are the edges
4270 of the target band pass/stop filter. With abundant tedious algebra,
4271 you can derive the above formulae yourself by substituting the
4272 transform for S into H(S)=S-x for a zero at x or H(S)=1/(S-x) for a
4273 pole at x, and converting the result into the form:
4275 H(S)=g prod(S-Xi)/prod(S-Xj)
4277 The transforms are from the references. The actual pole-zero-gain
4278 changes I derived myself.
4280 Please note that a pole and a zero at the same place exactly cancel.
4281 This is significant for High Pass, Band Pass and Band Stop filters
4282 which create numerous extra poles and zeros, most of which cancel.
4283 Those which do not cancel have a "fill-in" effect, extending the
4284 shorter of the sets to have the same number of as the longer of the
4285 sets of poles and zeros (or at least split the difference in the case
4286 of the band pass filter). There may be other opportunistic
4287 cancellations but I will not check for them.
4289 Also note that any pole on the unit circle or beyond will result in
4290 an unstable filter. Because of cancellation, this will only happen
4291 if the number of poles is smaller than the number of zeros and the
4292 filter is high pass or band pass. The analytic design methods all
4293 yield more poles than zeros, so this will not be a problem.
4297 Proakis & Manolakis (1992). Digital Signal Processing. New York:
4298 Macmillan Publishing Company.
4302 # name: <cell-element>
4306 usage: [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop)
4311 # name: <cell-element>
4318 # name: <cell-element>
4322 F = sgolay (p, n [, m [, ts]])
4323 Computes the filter coefficients for all Savitzsky-Golay smoothing
4324 filters of order p for length n (odd). m can be used in order to
4325 get directly the mth derivative. In this case, ts is a scaling factor.
4327 The early rows of F smooth based on future values and later rows
4328 smooth based on past values, with the middle row using half future
4329 and half past. In particular, you can use row i to estimate x(k)
4330 based on the i-1 preceding values and the n-i following values of x
4331 values as y(k) = F(i,:) * x(k-i+1:k+n-i).
4333 Normally, you would apply the first (n-1)/2 rows to the first k
4334 points of the vector, the last k rows to the last k points of the
4335 vector and middle row to the remainder, but for example if you were
4336 running on a realtime system where you wanted to smooth based on the
4337 all the data collected up to the current time, with a lag of five
4338 samples, you could apply just the filter on row n-5 to your window
4339 of length n each time you added a new sample.
4341 Reference: Numerical recipes in C. p 650
4343 See also: sgolayfilt
4347 # name: <cell-element>
4351 F = sgolay (p, n [, m [, ts]])
4352 Computes the filter coefficients for all Savi
4356 # name: <cell-element>
4363 # name: <cell-element>
4367 y = sgolayfilt (x, p, n [, m [, ts]])
4368 Smooth the data in x with a Savitsky-Golay smoothing filter of
4369 polynomial order p and length n, n odd, n > p. By default, p=3
4370 and n=p+2 or n=p+3 if p is even.
4372 y = sgolayfilt (x, F)
4373 Smooth the data in x with smoothing filter F computed by sgolay.
4375 These filters are particularly good at preserving lineshape while
4376 removing high frequency squiggles. Particularly, compare a 5 sample
4377 averager, an order 5 butterworth lowpass filter (cutoff 1/3) and
4378 sgolayfilt(x, 3, 5), the best cubic estimated from 5 points:
4380 [b, a] = butter(5,1/3);
4381 x=[zeros(1,15), 10*ones(1,10), zeros(1,15)];
4382 plot(sgolayfilt(x),"r;sgolayfilt;",...
4383 filtfilt(ones(1,5)/5,1,x),"g;5 sample average;",...
4384 filtfilt(b,a,x),"c;order 5 butterworth;",...
4385 x,"+b;original data;");
4391 # name: <cell-element>
4395 y = sgolayfilt (x, p, n [, m [, ts]])
4396 Smooth the data in x with a Savitsky-
4400 # name: <cell-element>
4407 # name: <cell-element>
4411 -- Function File: [PSI,X] = shanwavf (LB,UB,N,FB,FC)
4412 Compute the Complex Shannon wavelet.
4417 # name: <cell-element>
4421 Compute the Complex Shannon wavelet.
4425 # name: <cell-element>
4432 # name: <cell-element>
4436 -- Function File: Y = sigmoid_train(T, RANGES, RC)
4437 Evaluates a train of sigmoid functions at T.
4439 The number and duration of each sigmoid is determined from RANGES.
4440 Each row of RANGES represents a real interval, e.g. if sigmod `i'
4441 starts at `t=0.1' and ends at `t=0.5', then `RANGES(i,:) = [0.1
4442 0.5]'. The input RC is a array that defines the rising and
4443 falling time constants of each sigmoids. Its size must equal the
4446 Run `demo sigmoid_train' to some examples of the use of this
4453 # name: <cell-element>
4457 Evaluates a train of sigmoid functions at T.
4461 # name: <cell-element>
4468 # name: <cell-element>
4472 -- Function File: [B, A] = sos2tf (SOS, BSCALE)
4473 Convert series second-order sections to direct form H(z) =
4477 * SOS = matrix of series second-order sections, one per row:
4478 SOS = [B1.' A1.'; ...; BN.' AN.'], where
4479 `B1.'==[b0 b1 b2] and A1.'==[1 a1 a2]' for section 1, etc.
4480 b0 must be nonzero for each section.
4481 See `filter()' for documentation of the second-order
4482 direct-form filter coefficients Bi and Ai.
4484 * BSCALE is an overall gain factor that effectively scales the
4485 output B vector (or any one of the input Bi vectors).
4487 RETURNED: B and A are vectors specifying the digital filter H(z) =
4488 B(z)/A(z). See `filter()' for further details.
4490 See also: tf2sos zp2sos sos2pz zp2tf tf2zp
4496 # name: <cell-element>
4500 Convert series second-order sections to direct form H(z) = B(z)/A(z).
4504 # name: <cell-element>
4511 # name: <cell-element>
4515 -- Function File: [Z, P, G] = sos2zp (SOS, BSCALE)
4516 Convert series second-order sections to zeros, poles, and gains
4520 * SOS = matrix of series second-order sections, one per row:
4521 SOS = [B1.' A1.'; ...; BN.' AN.'], where
4522 `B1.'==[b0 b1 b2] and A1.'==[1 a1 a2]' for section 1, etc.
4523 b0 must be nonzero for each section. See `filter()' for
4524 documentation of the second-order direct-form filter
4525 coefficients Bi and Ai.
4527 * BSCALE is an overall gain factor that effectively scales any
4528 one of the input Bi vectors.
4531 * Z = column-vector containing all zeros (roots of B(z))
4532 * P = column-vector containing all poles (roots of A(z))
4533 * G = overall gain = B(Inf)
4536 [z,p,g] = sos2zp([1 0 1, 1 0 -0.81; 1 0 0, 1 0 0.49])
4537 => z = [i; -i; 0; 0], p = [0.9, -0.9, 0.7i, -0.7i], g=1
4539 See also: zp2sos sos2tf tf2sos zp2tf tf2zp
4545 # name: <cell-element>
4549 Convert series second-order sections to zeros, poles, and gains (pole
4554 # name: <cell-element>
4561 # name: <cell-element>
4565 usage: [S [, f [, t]]] = specgram(x [, n [, Fs [, window [, overlap]]]])
4567 Generate a spectrogram for the signal. This chops the signal into
4568 overlapping slices, windows each slice and applies a Fourier
4569 transform to determine the frequency components at that slice.
4571 x: vector of samples
4572 n: size of fourier transform window, or [] for default=256
4573 Fs: sample rate, or [] for default=2 Hz
4574 window: shape of the fourier transform window, or [] for default=hanning(n)
4575 Note: window length can be specified instead, in which case
4576 window=hanning(length)
4577 overlap: overlap with previous window, or [] for default=length(window)/2
4580 S is complex output of the FFT, one row per slice
4581 f is the frequency indices corresponding to the rows of S.
4582 t is the time indices corresponding to the columns of S.
4583 If no return value is requested, the spectrogram is displayed instead.
4586 x = chirp([0:0.001:2],0,2,500); # freq. sweep from 0-500 over 2 sec.
4587 Fs=1000; # sampled every 0.001 sec so rate is 1 kHz
4588 step=ceil(20*Fs/1000); # one spectral slice every 20 ms
4589 window=ceil(100*Fs/1000); # 100 ms data window
4590 specgram(x, 2^nextpow2(window), Fs, window, window-step);
4592 ## Speech spectrogram
4593 [x, Fs] = auload(file_in_loadpath("sample.wav")); # audio file
4594 step = fix(5*Fs/1000); # one spectral slice every 5 ms
4595 window = fix(40*Fs/1000); # 40 ms data window
4596 fftn = 2^nextpow2(window); # next highest power of 2
4597 [S, f, t] = specgram(x, fftn, Fs, window, window-step);
4598 S = abs(S(2:fftn*4000/Fs,:)); # magnitude in range 0<f<=4000 Hz.
4599 S = S/max(S(:)); # normalize magnitude so that max is 0 dB.
4600 S = max(S, 10^(-40/10)); # clip below -40 dB.
4601 S = min(S, 10^(-3/10)); # clip above -3 dB.
4602 imagesc(t, f, flipud(log(S))); # display in log scale
4604 The choice of window defines the time-frequency resolution. In
4605 speech for example, a wide window shows more harmonic detail while a
4606 narrow window averages over the harmonic detail and shows more
4607 formant structure. The shape of the window is not so critical so long
4608 as it goes gradually to zero on the ends.
4610 Step size (which is window length minus overlap) controls the
4611 horizontal scale of the spectrogram. Decrease it to stretch, or
4612 increase it to compress. Increasing step size will reduce time
4613 resolution, but decreasing it will not improve it much beyond the
4614 limits imposed by the window size (you do gain a little bit,
4615 depending on the shape of your window, as the peak of the window
4616 slides over peaks in the signal energy). The range 1-5 msec is good
4619 FFT length controls the vertical scale. Selecting an FFT length
4620 greater than the window length does not add any information to the
4621 spectrum, but it is a good way to interpolate between frequency
4622 points which can make for prettier spectrograms.
4624 After you have generated the spectral slices, there are a number of
4625 decisions for displaying them. First the phase information is
4626 discarded and the energy normalized:
4628 S = abs(S); S = S/max(S(:));
4630 Then the dynamic range of the signal is chosen. Since information in
4631 speech is well above the noise floor, it makes sense to eliminate any
4632 dynamic range at the bottom end. This is done by taking the max of
4633 the magnitude and some minimum energy such as minE=-40dB. Similarly,
4634 there is not much information in the very top of the range, so
4635 clipping to a maximum energy such as maxE=-3dB makes sense:
4637 S = max(S, 10^(minE/10)); S = min(S, 10^(maxE/10));
4639 The frequency range of the FFT is from 0 to the Nyquist frequency of
4640 one half the sampling rate. If the signal of interest is band
4641 limited, you do not need to display the entire frequency range. In
4642 speech for example, most of the signal is below 4 kHz, so there is no
4643 reason to display up to the Nyquist frequency of 10 kHz for a 20 kHz
4644 sampling rate. In this case you will want to keep only the first 40%
4645 of the rows of the returned S and f. More generally, to display the
4646 frequency range [minF, maxF], you could use the following row index:
4648 idx = (f >= minF & f <= maxF);
4650 Then there is the choice of colormap. A brightness varying colormap
4651 such as copper or bone gives good shape to the ridges and valleys. A
4652 hue varying colormap such as jet or hsv gives an indication of the
4653 steepness of the slopes. The final spectrogram is displayed in log
4654 energy scale and by convention has low frequencies on the bottom of
4657 imagesc(t, f, flipud(log(S(idx,:))));
4661 # name: <cell-element>
4665 usage: [S [, f [, t]]] = specgram(x [, n [, Fs [, window [, overlap]]]])
4670 # name: <cell-element>
4677 # name: <cell-element>
4681 -- Function File: S = square(T, DUTY)
4682 -- Function File: S = square(T)
4683 Generate a square wave of period 2 pi with limits +1/-1.
4685 If DUTY is specified, the square wave is +1 for that portion of
4689 duty cycle = ------------------
4692 See also: cos, sawtooth, sin, tripuls
4698 # name: <cell-element>
4702 Generate a square wave of period 2 pi with limits +1/-1.
4706 # name: <cell-element>
4713 # name: <cell-element>
4717 -- Function File: [NUM, DEN] = ss2tf (A, B, C, D)
4718 Conversion from transfer function to state-space. The state space
4724 is converted to a transfer function:
4734 # name: <cell-element>
4738 Conversion from transfer function to state-space.
4742 # name: <cell-element>
4749 # name: <cell-element>
4753 -- Function File: [POL, ZER, K] = ss2zp (A, B, C, D)
4754 Converts a state space representation to a set of poles and zeros;
4755 K is a gain associated with the zeros.
4761 # name: <cell-element>
4765 Converts a state space representation to a set of poles and zeros; K is
4770 # name: <cell-element>
4777 # name: <cell-element>
4781 -- Function File: [SOS, G] = tf2sos (B, A)
4782 Convert direct-form filter coefficients to series second-order
4785 INPUTS: B and A are vectors specifying the digital filter H(z) =
4786 B(z)/A(z). See `filter()' for documentation of the B and A filter
4789 RETURNED: SOS = matrix of series second-order sections, one per
4791 SOS = [B1.' A1.'; ...; BN.' AN.'], where
4792 `B1.'==[b0 b1 b2] and A1.'==[1 a1 a2]' for section 1, etc.
4793 b0 must be nonzero for each section (zeros at infinity not
4794 supported). BSCALE is an overall gain factor that effectively
4795 scales any one of the Bi vectors.
4800 [sos,g] = tf2sos(B,A)
4804 1.00000 0.61803 1.00000 1.00000 0.60515 0.95873
4805 1.00000 -1.61803 1.00000 1.00000 -1.58430 0.95873
4806 1.00000 1.00000 -0.00000 1.00000 0.97915 -0.00000
4810 See also: sos2tf zp2sos sos2pz zp2tf tf2zp
4816 # name: <cell-element>
4820 Convert direct-form filter coefficients to series second-order sections.
4824 # name: <cell-element>
4831 # name: <cell-element>
4835 -- Function File: [A, B, C, D] = tf2ss (NUM, DEN)
4836 Conversion from transfer function to state-space. The state space
4841 is obtained from a transfer function:
4846 The state space system matrices obtained from this function will
4847 be in observable companion form as Wolovich's Observable Structure
4853 # name: <cell-element>
4857 Conversion from transfer function to state-space.
4861 # name: <cell-element>
4868 # name: <cell-element>
4872 -- Function File: [ZER, POL, K] = tf2zp (NUM, DEN)
4873 Converts transfer functions to poles-and-zero representations.
4875 Returns the zeros and poles of the system defined by NUM/DEN. K
4876 is a gain associated with the system zeros.
4881 # name: <cell-element>
4885 Converts transfer functions to poles-and-zero representations.
4889 # name: <cell-element>
4896 # name: <cell-element>
4901 [Pxx,freq] = tfe(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
4903 Estimate transfer function of system with input "x" and output "y".
4904 Use the Welch (1967) periodogram/FFT method.
4905 Compatible with Matlab R11 tfe and earlier.
4906 See "help pwelch" for description of arguments, hints and references
4907 --- especially hint (7) for Matlab R11 defaults.
4911 # name: <cell-element>
4916 [Pxx,freq] = tfe(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
4921 # name: <cell-element>
4928 # name: <cell-element>
4933 [Pxx,freq]=tfestimate(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)
4935 Estimate transfer function of system with input "x" and output "y".
4936 Use the Welch (1967) periodogram/FFT method.
4937 See "help pwelch" for description of arguments, hints and references.
4941 # name: <cell-element>
4946 [Pxx,freq]=tfestimate(x,y,Nfft,Fs,window,overlap,range,plot_type,detr
4950 # name: <cell-element>
4957 # name: <cell-element>
4961 usage: w = triang (L)
4963 Returns the filter coefficients of a triangular window of length L.
4964 Unlike the bartlett window, triang does not go to zero at the edges
4965 of the window. For odd L, triang(L) is equal to bartlett(L+2) except
4966 for the zeros at the edges of the window.
4970 # name: <cell-element>
4974 usage: w = triang (L)
4979 # name: <cell-element>
4986 # name: <cell-element>
4990 usage: y = tripuls(t, w, skew)
4992 Generate a triangular pulse over the interval [-w/2,w/2), sampled at
4993 times t. This is useful with the function pulstran for generating a
4996 skew is a value between -1 and 1, indicating the relative placement
4997 of the peak within the width. -1 indicates that the peak should be
4998 at -w/2, and 1 indicates that the peak should be at w/2.
5001 fs = 11025; # arbitrary sample rate
5002 f0 = 100; # pulse train sample rate
5003 w = 0.3/f0; # pulse width 3/10th the distance between pulses
5004 auplot(pulstran(0:1/fs:4/f0, 0:1/f0:4/f0, 'tripuls', w), fs);
5010 # name: <cell-element>
5014 usage: y = tripuls(t, w, skew)
5019 # name: <cell-element>
5026 # name: <cell-element>
5030 -- Function File: W = tukeywin (L, R)
5031 Return the filter coefficients of a Tukey window (also known as the
5032 cosine-tapered window) of length L. R defines the ratio between
5033 the constant section and and the cosine section. It has to be
5034 between 0 and 1. The function returns a Hanning window for R egals
5035 0 and a full box for R egals 1. By default R is set to 1/2.
5037 For a definition of the Tukey window, see e.g. Fredric J. Harris,
5038 "On the Use of Windows for Harmonic Analysis with the Discrete
5039 Fourier Transform, Proceedings of the IEEE", Vol. 66, No. 1,
5040 January 1978, Page 67, Equation 38.
5045 # name: <cell-element>
5049 Return the filter coefficients of a Tukey window (also known as the
5054 # name: <cell-element>
5061 # name: <cell-element>
5065 -- Function File: Y = upsample (X, N)
5066 -- Function File: Y = upsample (X, N, OFFSET)
5067 Upsample the signal, inserting n-1 zeros between every element.
5069 If X is a matrix, upsample every column.
5071 If OFFSET is specified, control the position of the inserted
5072 sample in the block of n zeros.
5074 See also: decimate, downsample, interp, resample, upfirdn
5080 # name: <cell-element>
5084 Upsample the signal, inserting n-1 zeros between every element.
5088 # name: <cell-element>
5095 # name: <cell-element>
5099 -- Function File: [W] = welchwin(L,C)
5100 Returns a row vector containing a Welch window, given by
5101 W(n)=1-(n/N-1)^2, n=[0,1, ... L-1]. Argument L is the length of
5102 the window. Optional argument C specifies a "symmetric" window
5103 (the default), or a "periodic" window.
5105 A symmetric window has zero at each end and maximum in the middle;
5106 L must be an integer larger than 2. `if c=="symmetric", N=(L-1)/2'
5108 A periodic window wraps around the cyclic interval [0,1, ... L-1],
5109 and is intended for use with the DFT (functions fft(),
5110 periodogram() etc). L must be an integer larger than 1. `if
5111 c=="periodic", N=L/2'.
5113 See also: blackman, kaiser
5119 # name: <cell-element>
5123 Returns a row vector containing a Welch window, given by
5124 W(n)=1-(n/N-1)^2, n=[
5128 # name: <cell-element>
5135 # name: <cell-element>
5139 -- Function File: W = window (F, N, OPTS)
5140 Create a N-point windowing from the function F. The function F can
5141 be for example `@blackman'. Any additional arguments OPT are
5142 passed to the windowing function.
5147 # name: <cell-element>
5151 Create a N-point windowing from the function F.
5155 # name: <cell-element>
5162 # name: <cell-element>
5166 -- Function File: [Y] = wkeep(X,L,OPT)
5167 Extract the elements of x of size l from the center, the right or
5173 # name: <cell-element>
5177 Extract the elements of x of size l from the center, the right or the
5182 # name: <cell-element>
5189 # name: <cell-element>
5193 -- Function File: [Y] = wrev(X)
5194 Reverse the order of the element of the vector x.
5196 See also: flipud, fliplr
5202 # name: <cell-element>
5206 Reverse the order of the element of the vector x.
5210 # name: <cell-element>
5217 # name: <cell-element>
5221 usage: [R, lag] = xcorr (X [, Y] [, maxlag] [, scale])
5223 Estimate the cross correlation R_xy(k) of vector arguments X and Y
5224 or, if Y is omitted, estimate autocorrelation R_xx(k) of vector X,
5225 for a range of lags k specified by argument "maxlag". If X is a
5226 matrix, each column of X is correlated with itself and every other
5229 The cross-correlation estimate between vectors "x" and "y" (of
5230 length N) for lag "k" is given by
5231 R_xy(k) = sum_{i=1}^{N}{x_{i+k} conj(y_i),
5232 where data not provided (for example x(-1), y(N+1)) is zero.
5235 X [non-empty; real or complex; vector or matrix] data
5237 Y [real or complex vector] data
5238 If X is a matrix (not a vector), Y must be omitted.
5239 Y may be omitted if X is a vector; in this case xcorr
5240 estimates the autocorrelation of X.
5242 maxlag [integer scalar] maximum correlation lag
5243 If omitted, the default value is N-1, where N is the
5244 greater of the lengths of X and Y or, if X is a matrix,
5245 the number of rows in X.
5247 scale [character string] specifies the type of scaling applied
5248 to the correlation vector (or matrix). is one of:
5249 'none' return the unscaled correlation, R,
5250 'biased' return the biased average, R/N,
5251 'unbiased' return the unbiassed average, R(k)/(N-|k|),
5252 'coeff' return the correlation coefficient, R/(rms(x).rms(y)),
5253 where "k" is the lag, and "N" is the length of X.
5254 If omitted, the default value is "none".
5255 If Y is supplied but does not have the ame length as X,
5256 scale must be "none".
5259 R array of correlation estimates
5260 lag row vector of correlation lags [-maxlag:maxlag]
5262 The array of correlation estimates has one of the following forms.
5263 (1) Cross-correlation estimate if X and Y are vectors.
5264 (2) Autocorrelation estimate if is a vector and Y is omitted,
5265 (3) If X is a matrix, R is an matrix containing the cross-
5266 correlation estimate of each column with every other column.
5267 Lag varies with the first index so that R has 2*maxlag+1
5268 rows and P^2 columns where P is the number of columns in X.
5269 If Rij(k) is the correlation between columns i and j of X
5270 R(k+maxlag+1,P*(i-1)+j) == Rij(k)
5271 for lag k in [-maxlag:maxlag], or
5272 R(:,P*(i-1)+j) == xcorr(X(:,i),X(:,j)).
5273 "reshape(R(k,:),P,P)" is the cross-correlation matrix for X(k,:).
5278 # name: <cell-element>
5282 usage: [R, lag] = xcorr (X [, Y] [, maxlag] [, scale])
5287 # name: <cell-element>
5294 # name: <cell-element>
5298 -- Function File: C = xcorr2 (A)
5299 -- Function File: C = xcorr2 (A, B)
5300 -- Function File: C = xcorr2 (..., SCALE)
5301 Compute the 2D cross-correlation of matrices A and B.
5303 If B is not specified, it defaults to the same matrix as A, i.e.,
5304 it's the same as `xcorr(A, A)'.
5306 The optional argument SCALE, defines the type of scaling applied
5307 to the cross-correlation matrix (defaults to "none"). Possible
5311 Scales the raw cross-correlation by the maximum number of
5312 elements of A and B involved in the generation of any element
5317 Scales the raw correlation by dividing each element in the
5318 cross-correlation matrix by the number of products A and B
5319 used to generate that element
5323 Normalizes the sequence so that the largest cross-correlation
5324 element is identically 1.0.
5328 No scaling (this is the default).
5330 See also: conv2, corr2, xcorr
5336 # name: <cell-element>
5340 Compute the 2D cross-correlation of matrices A and B.
5344 # name: <cell-element>
5351 # name: <cell-element>
5355 usage: [c, lag] = xcov (X [, Y] [, maxlag] [, scale])
5357 Compute covariance at various lags [=correlation(x-mean(x),y-mean(y))].
5360 Y: if specified, compute cross-covariance between X and Y,
5361 otherwise compute autocovariance of X.
5362 maxlag: is specified, use lag range [-maxlag:maxlag],
5363 otherwise use range [-n+1:n-1].
5365 'biased' for covariance=raw/N,
5366 'unbiased' for covariance=raw/(N-|lag|),
5367 'coeff' for covariance=raw/(covariance at lag 0),
5368 'none' for covariance=raw
5369 'none' is the default.
5371 Returns the covariance for each lag in the range, plus an
5372 optional vector of lags.
5376 # name: <cell-element>
5380 usage: [c, lag] = xcov (X [, Y] [, maxlag] [, scale])
5385 # name: <cell-element>
5392 # name: <cell-element>
5396 -- Function File: X0 = zerocrossing (X, Y)
5397 Estimates the points at which a given waveform y=y(x) crosses the
5398 x-axis using linear interpolation.
5400 See also: fzero, roots
5406 # name: <cell-element>
5410 Estimates the points at which a given waveform y=y(x) crosses the
5415 # name: <cell-element>
5422 # name: <cell-element>
5426 -- Function File: [SOS, G] = zp2sos (Z, P)
5427 Convert filter poles and zeros to second-order sections.
5430 * Z = column-vector containing the filter zeros
5431 * P = column-vector containing the filter poles
5432 * G = overall filter gain factor
5435 * SOS = matrix of series second-order sections, one per row:
5436 SOS = [B1.' A1.'; ...; BN.' AN.'], where
5437 `B1.'==[b0 b1 b2] and A1.'==[1 a1 a2]' for section 1, etc.
5438 b0 must be nonzero for each section.
5439 See `filter()' for documentation of the second-order
5440 direct-form filter coefficients Bi and %Ai, i=1:N.
5442 * BSCALE is an overall gain factor that effectively scales any
5443 one of the Bi vectors.
5446 [z,p,g] = tf2zp([1 0 0 0 0 1],[1 0 0 0 0 .9]);
5447 [sos,g] = zp2sos(z,p,g)
5450 1.0000 0.6180 1.0000 1.0000 0.6051 0.9587
5451 1.0000 -1.6180 1.0000 1.0000 -1.5843 0.9587
5452 1.0000 1.0000 0 1.0000 0.9791 0
5457 See also: sos2pz sos2tf tf2sos zp2tf tf2zp
5463 # name: <cell-element>
5467 Convert filter poles and zeros to second-order sections.
5471 # name: <cell-element>
5478 # name: <cell-element>
5482 -- Function File: [A, B, C, D] = zp2ss (ZER, POL, K)
5483 Conversion from zero / pole to state space.
5488 Vectors of (possibly) complex poles and zeros of a transfer
5489 function. Complex values must come in conjugate pairs (i.e.,
5490 x+jy in ZER means that x-jy is also in ZER).
5493 Real scalar (leading coefficient).
5500 The state space system, in the form:
5508 # name: <cell-element>
5512 Conversion from zero / pole to state space.
5516 # name: <cell-element>
5523 # name: <cell-element>
5527 -- Function File: [NUM, DEN] = zp2tf (ZER, POL, K)
5528 Converts zeros / poles to a transfer function.
5533 Vectors of (possibly complex) poles and zeros of a transfer
5534 function. Complex values must appear in conjugate pairs.
5537 Real scalar (leading coefficient).
5542 # name: <cell-element>
5546 Converts zeros / poles to a transfer function.
5550 # name: <cell-element>
5557 # name: <cell-element>
5561 usage: zplane(b [, a]) or zplane(z [, p])
5563 Plot the poles and zeros. If the arguments are row vectors then they
5564 represent filter coefficients (numerator polynomial b and denominator
5565 polynomial a), but if they are column vectors or matrices then they
5566 represent poles and zeros.
5568 This is a horrid interface, but I didn't choose it; better would be
5569 to accept b,a or z,p,g like other functions. The saving grace is
5570 that poly(x) always returns a row vector and roots(x) always returns
5571 a column vector, so it is usually right. You must only be careful
5572 when you are creating filters by hand.
5574 Note that due to the nature of the roots() function, poles and zeros
5575 may be displayed as occurring around a circle rather than at a single
5578 The transfer function is
5580 B(z) b0 + b1 z^(-1) + b2 z^(-2) + ... + bM z^(-M)
5581 H(z) = ---- = --------------------------------------------
5582 A(z) a0 + a1 z^(-1) + a2 z^(-2) + ... + aN z^(-N)
5584 b0 (z - z1) (z - z2) ... (z - zM)
5585 = -- z^(-M+N) ------------------------------
5586 a0 (z - p1) (z - p2) ... (z - pN)
5588 The denominator a defaults to 1, and the poles p defaults to [].
5592 # name: <cell-element>
5596 usage: zplane(b [, a]) or zplane(z [, p])