Home

Higher-Order Spectral Analysis Toolbox User's Guide

1. svals of c2 40 0 Poo OR Qa 0 ji oo 100 0 10 100 50 0 50 100 0 0 2 2 20 UK 3 20 40 40 100 50 0 50 100 100 50 0 50 100 0 E 20 INS cz 5 20 CN 40 40 100 50 0 50 100 100 50 0 50 100 0 0 E zo eS T er 40 20 100 50 0 50 100 100 50 0 50 100 Figure 1 14 Singular Values of Spatial Correlation Matrix and Estimated Angular Spectra doa 1 77 1 Tutorial load doal spec theta doa ymat 0 5 1 2 The singular value plot indicates the possible presence of two sources the estimated angular spectra show peaks at around the true bearings of 15 and 25 degrees notice that the spatially correlated Gaussian noise source has been virtually suppressed by using fourth order cumulants The vertical lines and circles denote the true bearings the cross denotes the virtual bearing of the spatially correlated noise The bearings estimated by ESPRIT should be 15 0066 and 25 1361 svals of c4 400 200 0 1 0 5 10 0 E A o b wee 50 100 50 0 50 100 0 us Med 40 100 50 0 50 100 0 gt E 20 EIN 40 100 50 0 50 100 ar i M 50 n 100 100 50 0 50 100 0 20 C 40 100 50 0 50 100 0 x 50 100 50 0 50 100 BUM 20 e vs 40 100 50 0 50 100 Figure 1 15 Singular Values of Spatial Cumulant Matrix and
2. N I E e L 5 2000 ee j gt SSS o semen D Lo 1500 4 z 500F z 4 0 i 1 i i 1 i 0 20 40 60 80 100 120 time in ms Figure 1 44 Laughter Data Spectrogram The spectrogram indicates that the speech signal is essentially harmonic in nature it also shows that the frequencies appear to be approximately harmonically related Since the STFT is resolution limited we can use parametric methods such as har mest to obtain accurate estimates we should of course keep in mind the nonstationary nature of the data 1 123 1 Tutorial 1 124 The 1 1 panel of Figure 1 45 shows the singular values of the covariance matrix as estimated by har mest we set nl ag 25 and nf f t 512 The singular values are essentially zero after p 12 hence we selected an order of 12 and obtained the various estimates shown in the figure Wecan useroots tocomputethe roots of the AR polynomials estimated by the AR and Min Norm methods angle roots ar 2 pi yields the normalized frequenciesandabs roots ar yieldstheradius thetop two panels of Figure 1 47 show the locations of the roots of thetwo AR polynomials The roots closest to the unit circle were found to be approximately at 0 0702 0 1329 0 1962 Hz for the AR method less dominant roots at 0 2480 0 3706 and 0 2012 for the Min Norm method the dominant roots are at frequencies 0 0710 0 1945 0 1272 and 0 137
3. where S is the power spectrum of process x n and Syy is the cross spectrum of the processes x n and y n This relationship assumes that x n is symmetrically distributed The quadratic part q k l is estimated in the frequency domain from Syxx 1 02 2 2Q 01 05 55 091 55 2 Sy 2 5 1 w2 E y n which is obtained under the assumption that q k I q I k and that x n is Gaussian The cross bispectrum Sy 1 2 is estimated via the direct FFT method using function bi specdx nl pow bi specdx 1 Tick L J The estimation of transfer functions of quadraticsystems Technometrics 3 563 67 Nov 1961 pickpeak Purpose Syntax Description Picks peaks subject to a separation criterion loc val pickpeak spec npicks rdi ff pi ckpeak picksthelargest npi cks peaks in the data vector spec such that the peaks are separated by at least r di f f samples The default values arenpi cks 2andrdiff 5 Ifspec isa matrix each column is treated independently loc is the matrix of locations indices of the picked peaks the kth column corresponds to the kth column of s pec val isthe matrix of the amplitudes of the picked peaks the kth column corresponds to the kth column of s pec A Oin location i j of array oc or a NaN in array val indicates that the jth data vector has less than i peaks with a separation of rdi f f or more 2 69 qpcgen Purpose Syntax Description
4. 1 128 1 63 1 Tutorial 1 64 where 6 is the bearing of the kth source and d is the location of the nth source with reference to some arbitrary origin For a uniformly spaced linear array with spacing d we can set dp nd The observed signals are thus given by P dsin Q Ym N Y calm exp j2nn 5 8 Wm N m 1 M a 1 129 Comparing 1 129 with 1 125 we see that the term dsin 6 A plays the role of frequency and oy m plays the role of the complex amplitude aexp jO_ m Timesamples indexed by n in 1 125 correspond to spatial samples or sensor numbers in 1 129 The realization number m in 1 125 is interpreted as the snapshot number in 1 129 In 1 124 we assume 0 lt f lt 1in order to ensure that thereis noaliasing Similarly in 1 129 we assume that A gt 2d to ensure that there is no aliasing a common assumption is A 2d half wavelength spacing Thus a solution to the harmonic retrieval problem in 1 125 also applies tothe DOA problem in 1 129 In the former we estimate power spectra as functions of frequency in the latter we estimate angular spectra as functions of angle 0 Thefollowing table states the duality between the harmonic retrieval and the DOA problems Harmonic Retrieval DOA Realization Snapshot Time sample Sensor sample Frequency f Angular frequency d sin 6 A Random phases Random complex amplitudes 0O f 1 d lt A 2 Temporal correlation matrix Spat
5. 1 86 The 2 D FT of 1 164 yields the bispectrum S 304 A gt Na 5 gt o4 05055 i l 84 Fed fib 804 fio fi 304 FDA fi 804 f19 80 5 fi 3504 FDSA f 84 fe fr 4804 580 f 804 75 5002 fi 3904 FDSA fio 804 f19 8 3 fio 3904 f5 80 f 804 fiA f nues It is evident that the nonredundant region of the bispectrum shows peaks only at the phase and frequency coupled bifrequencies fis Ti TheFT of the diagonal slice in 1 165 on the other hand displays peaks at each of thethree frequencies involved in the phase coupling Consider with Ng 1 We have omitting unnecessary superscripts 3 1 C3x T 1 5019203 Y cos 2n fx er 1 167 where f3 f2 f Note that the diagonal slice of the third order cumulant is expressed as the sum of three harmonics From our discussions on AR models for harmonics we note that C3 t t satisfies a self driven AR 6 model whose roots are at exp j2zf k 1 2 3 If we estimate the AR polynomial we can compute the parametric bispectrum S3 Np 12 A AQP A n4 n2 1 168 The parametric bispectrum estimator is implemented in routineqpctor N onlinear Processes Examples load qpc arvec bSpec qpctor zmat 18 12 Maximum of bispectrum B 0 1484 0 1016 4239 You should see the display on Figure 1 19 Note that f4 f 2 0 25 which corresponds to the third peak in the amplitude spectrum
6. Assume that we have good perfect estimates of the power spectrum and consider the sample estimate of the squared bicoherence A 2 Sxxx f 4 f2 I 2 bi i ee oe a eee i ee pe icxxx f 1 fj Sox f4 f2 S2x f1 S2x f2 1 30 Polyspectra and Linear Processes t has been established in 6 that sample estimates of the bispectrum using conventional methods are asymptotically complex Gaussian additionally the estimates at different frequencies are uncorrelated provided the frequency separation is greater than thereciprocal of the effective window length If Sxxx is Gaussian distributed we know that IS is chi squared distributed with two degrees of freedom If Soo fq f2 20 then the statistic in 1 30 is a central chi squared r v with two degrees of freedom The squared bicoherence is summed over the P points in the nonredundant region details are given in 24 The resulting statistic S is x distributed with 2P degrees of freedom Hence it is easy to devise a statistical test to determine whether the observed S is consistent with a central chi squared distribution this consistency is reported as probability of false alarm value that is the probability that we will be wrongin assuming that the data have a nonzero bispectrum If this probability is small say 0 95 we accept the assumption of zero bispectrum that is we cannot reject the Gaussianity assumption This test is implemented in routine gl stat Let
7. If y isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset to zero andsamp seg is set tothe row dimension Cumulants are estimated from each realization and then averaged The estimated AR and MA parameters are returned in the vectors avec and bvec respectively 2 11 armarts Algorithm See Also 2 12 The signal is assumed to be described by x n p q Y a k x n k b k u n k ke k 0 x n g n y n where u n is i i d non Gaussian and g n is Gaussian noise The AR parameters are obtained as the least squares solution to the normal equations given by a k R m k 0 m gt q a k C3 m k p 0 m gt q 0 k k a k C m k p 0 0 m gt q k 0 0 p where a 0 21 m q 1 maxlag andp q p q Since one must have atleast m q 1 q p 1 2 it follows that the cumulant lags involved must include q 1 p q p this leads to the default value of max ag If the estimated AR parameters are exact then the residual time series is an MA q process p q p yin ackyy n k b kyu n k 5 a k g n k 0 k 0 k 0 hence the MA parameters can now be determined via the algorithms described in maest Note that the even if g n is white the additive noiseiny n is no longer white We assume that b 0 h 0 1 arrcest maest cumest armaqs armarts References 1 Giannakis G B and J M Mendel
8. X n p 1 90 Matrix is Hermitian but not Toeplitz and the relationship in 1 69 is not satisfied for finite N 1 53 1 Tutorial Inthe forward backward least squares F BLS problem we minimizeP y Pf P assuming cy k ap P k This leads tothe deterministic normal equations 1 91 where is given by N H 2 H L x n x n X n x n Tz 1 92 where m H X n x n p x n 1 93 If is nonsingular we can solve 1 91 directly if it is singular then the Min Norm solution is usually used The solution a may be interpreted as providing an AR p model for the data hence we can compute the spectrum of the process x n in terms of the estimated AR parameters Additional details may be found in 34 35 Adaptive Linear Prediction The Burg algorithm provides a recursive solution to minimizing the mean square error over a timeinterval of length n and is appropriate when the data may be considered to be stationary over that timeinterval Often the data are nonstationary in this case we can estimate power spectra over a sliding window of fixed length however as we have seen before the resulting estimate will have poor resolution and high variance An alternativeis to fit a parametric model and to allow the model parameters to vary with time F or example we may fit a time varying AR model that is we perform time varying linear prediction The deterministic normal equatio
9. morder z n q 0 y n q 1 sign y n 2 y n 3 y n 4 y n 3s n y n where for 0 lt A lt 1 s n s n 1 ay2 n s 1 y 1 and for 1 n s n Y y k k 1 1 Swami A and M Mendel Adaptive Cumulant Based Estimation of ARMA Parameters Proc Amer Control Conf ACC 88 Atlanta GA 2114 19 J une 1988 2 57 maest Purpose Syntax Description Algorithm 2 58 Estimates MA parameters using the modified GM method bvec maest y q bvec maest y q norder samp seg overlap flag maest estimates the parameters of the MA q process y via the modified GM method based on cumulants of order nor der and the autocorrelation y is the data vector or matrix q istheMA order norder should be 3 or 4 the default value is 3 Variablessamp seg overlap andf lag control the manner in which sample cumulants are esti mated samp seg specifies the number of samples per segment the default valueis the length of the time series overlap specifies the percentage overlap between segments maximum allowed value is 99 default valueis O Iff I ag is biased then biased sample esti mates are computed def aul t ifthe first letter is not b unbiased estimates are computed Ify isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset tozero andsamp seg is set totherow dimension Cumulants are estimated from each realiza
10. 40 amp 40 5 F c30 PN lt 30 co co EE 20 e 20 10 10 0 2 0 1 0 0 1 0 2 0 2 0 1 0 0 1 0 2 frequency frequency ws WS 60 60 u 50 u 50 E o amp 40 40 E F c30 lt 30 20 20 10 10 0 2 0 1 0 0 1 0 2 frequency 0 2 0 1 0 0 1 0 2 frequency Figure 1 20 Wigner Spectra wig2 Signal s2 is a harmonic with a nominal frequency of 0 05 Hz which has been multiplied by a Gaussian window and is centered at n 220 Note that the WS of the signal is concentrated around its nominal center frequency note also the interaction between the energies at the positive and negative frequencies giving rise to the interference terms around D C By using the analytic version of thesignal the negative frequency terms are suppressed this also suppresses the distortion terms around D C 1 93 1 Tutorial 1 94 Signal 3 is a harmonic with a nominal frequency of 0 15 Hz which has been multiplied by a Gaussian window and is centered at n 50 Note that the WS of the signal is concentrated around its nominal center frequency In the last panel we compute the WS of the sum signal s4 2 453 note the cross terms in the WS some of these can be eliminated by using routinewi g2c Examples load wi gdat wi g2c 54 Compare the display in Figure 1 21 with the 2 2 panel of Figure 1 20 note that the cross terms have been suppressed CWS s 0 05 60r B o T o T time in samples
11. Function Tables Higher Order Spectrum Estimation Parametric Methods Function Purpose ar maqs ARMA parameter estimation q slice algorithm armarts ARMA parameter estimation residual time series algorithm armasyn Generates ARMA synthetics arorder AR order estimation AR and ARMA processes arrcest AR parameter estimation bispect Computes theoretical bispectrum of an ARMA process cumt rue Computes true cumulants of ARMA processes maest MA parameter estimation maorder MA order estimation rpiid Generates sequence of i i d random variables trispect Computes theoretical trispectrum of an ARMA process Quadratic Phase Coupling QPC Function Purpose qpcgen Generates synthetics for the QPC problem qpctor Parametric QPC detection 2 Reference Second Order Volterra Systems Function Purpose nl gen Computes output of a second order Volterra system nl pow Parameter estimation arbitrary inputs Powers method nltick Parameter estimation Gaussian inputs Tick s method Harmonic Retrieval Function Purpose har mest Estimates frequencies of sinusoidal signals harmgen Generates synthetics for harmonic retrieval probl em Time Delay Estimation TDE Function Purpose t de TDE based on cross cumulants t deb TDE based on cross bispectrum t degen Generates synthetics for the TDE problem tder TDE based on cross correlation Array Processing Direction of Arrival DOA
12. X N wg Pen 1 125 where m denotes the realization number If only a single realization is available then the phases are nonrandom constants Since 2cos 0 e e we note that a signal consisting of p real harmonics can be handled by the model in 1 125 with 2p frequencies f K 1 p The problem is to determine p and then to estimate the frequencies amplitudes and phases of each harmonic Intuitively we expect the power spectrum of x n to consist of peaks at f f perhaps masked by the spectrum of the noise The classical power spectrum estimators Blackman T ukey and Welch can be used only if the frequencies are well separated that is min fm fa gt 1 L where L is the length of the estimated autocorrelation sequence in the Blackman Tukey method or the block size in the Welch method hence parametric models are useful for short data lengths The model in 1 125 is similar to the model for the observed data in the array processing direction of arrival DOA problem This problem arises in various applications such as radar sonar geophysics and analysis of EEG and ECG signals An array of N independent sensors provides samples of the data at different points in space Data are collected simultaneously from the set of sensors the set collected at a single instant of time is called a snapshot and is a vector of N observations Thesignal recorded at each sensor consists of sensor Harmonic Processes and DO A
13. bicepsf 1 Pan R and C L Nikias The complex cepstrum of higher order cumulants and non minimum phase system identification IEEE Trans ASSP Vol 36 pp 186 205 F eb 1988 bicepsf Purpose Syntax Description Estimates impulse response via frequency domain bicepstral method hest ceps bicepsf y nlag samp seg overlap flag nfft wind bi cepsf estimates the impulse response of the linear process y using the bicepstrum F ast Fourier Transform FFT method y is the data vector or matrix nl ag specifies the number of cumulant lags to be computed the third order cumulants of y C3 m n will be estimated for nlag lt m n lt nlag This parameter must be specified A useful rule of thumb istosetnl ag tonsamp 10 wherensamp isthelength of thetimeseries y Variablessamp seg overlap andflag control the manner in which sample cumulants are esti mated samp seg Specifies the number of samples per segment or record Its default value is the row dimension of y if y is a row vector the column dimension is used as the default value overlap specifies the percentage overlap between segments The default value is 0 The allowed range is 0 99 Ify isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset to zero andsamp seg is set tothe row dimension Ifflag is b then biased sample estimates of cumulants are estimated default ifthe first letter is not b
14. denotes the ensemble expectation operator The power spectrum is formally defined as the Fourier Transform FT of the autocorrelation sequence the Wiener K hintchine theorem P Y Ryanexpoj2nfm PET 1 2 where f denotes the frequency An equivalent definition is given by Pyx f s ELX f X f 1 3 where X f is the Fourier Transform of x n X f Y x n exp j2nfn Diss 1 4 A sufficient but not necessary condition for the existence of the power spectrum is that the autocorrelation be absolutely summable The power spectrum is real valued and nonnegative that is P f 0 if x n is real valued then the power spectrum is also symmetric that is P y f P xxl f As we shall see next higher order moments are natural generalizations of the autocorrelation and cumulants are specific nonlinear combinations of these moments Thefirst order cumulant of a stationary process is the mean C4 E amp t The higher order cumulants areinvariant to a shift of mean hence it is convenient to define them under the assumption of zero mean if the process has nonzero mean we subtract the mean and then apply the following definitions to the resulting process The second third and fourth order cumulants of a zero mean stationary process are defined by 4 Polyspectra and Linear Processes C3 k E x n x n k 1 5 x C3 k 1 E x n x n k x n 1 1 6 Cy k l m E x n x n k x n 1 x
15. eliminates cross terms in the WT of multicomponent signals by applying the Choi Williams filter x is the signal and must be a vector nfft specifies the FFT length and hence the frequency resolution this parameter must be greater than four times the length of the signal x in order to avoid aliasing The default value is the power of 2 equal to or just greater than four times the signal length sigma isthe parameter in the Choi Williams window and controls the amount of cross term suppression very large values result in no Suppression whereas very small values result in loss of signal terms The default value is 0 05 flag ifflag is nonzero and x is real valued the analytic form of x is used instead of x the default valueis 1 If x t is a real valued signal with Hilbert transform y t then theanalyticsignal is defined as the complex valued signal z t x t jy t wx isthe computed WT The rows correspond to time samples and the columns to frequencies waxis isa vector whose entries are the frequencies associated with the columns of the WT For a complex signal the Wigner trispectrum can be defined in two different ways here we implement the symmetric Wigner trispectrum 1 Definethe instantaneous fourth order product ra t t1 t2 73 x t t x t t T1 X t T v2 x t T 73 where t 14 12 13 4 Define the smoothing kernel 0 11 15 T3 exp e r r ra signia wig4c See Also
16. noise as well as signals from p sources assumed to be uncorrelated with one another The signal recorded at the nth sensor and mth snapshot is p Ym N Y AkSy m Tk n Wm N k 1 where s m is the signal emitted by the kth source and tk n is the signal delay with respect to some arbitrary reference point at the nth sensor Typically the sensors are assumed to be isotropic and the source signals are assumed to be plane waves originating from the far field of the array The signals are assumed to be narrow band and as such they are characterized by a single frequency fo or a single wavelength X c fo where cis the speed of signal propagation The received signal can now be written as p ym n a m exp jo k n wq k 1 1 126 where o4 m is the amplitude of thekth source signal at snapshot m and 6 k n is the phase delay of the kth source signal at the nth sensor and depends upon the array geometry In a circular array the sensors are assumed to be equally distributed in angle on the circumference of a circle with radius R Let the bearings 6 be measured with respect to the radius vector passing through a sensor The phase delay is then given by 23 _ 2nR 2nn 2nn o k n AER cos 252 9 cos Er 1 127 For a linear array the bearing of the plane wave is defined as the angle between the plane wave front and the normal to the array The phase delay is given by 2nd o k n X sin
17. 001 ng 20 f 50 o 0 01 for signal s2 and T 001 ng 50 f 150 o 0 01 for signal s3 Signal s4 is the sum of s2 ands3 Note that signalss2 ands3 are harmonics at frequencies 0 05 and 0 15 Hz which have been multiplied by a Gaussian shaped window function Signals s1 s2 s3 and s4 are displayed in Figure 1 52 si s2 1 1 0 8 0 5 6 04 0 02 0 5 9 0 2 1 0 4 0 10 20 30 40 0 20 40 60 80 s3 S4 S2 s3 1 15 1 0 5 0 5 0 0 0 5 0 5 1 zl 1 5 0 20 40 60 80 0 20 40 60 80 Figure 1 52 Signals in file wigdat mat 1 138 References References Allen J B Short time spectral analysis synthesis and modification by discrete Fourier transform IEEE Trans ASSP Vol 25 pp 235 38 J une 1977 Blackman R B and J W Tukey TheM easurement of Power Spectra New York Dover 1959 Boudreaux Bartels G Time F requency Signal Processing Algorithms Analysis and Synthesis Using Wigner Distributions Ph D dissertation Rice University 1984 Brillinger D R An Introduction to Polyspectra Ann Math Stat Vol 36 pp 1351 74 1965 Brillinger D R TimeSeries Data Analysis and Theory McGraw Hill 1981 Brillinger D R and M Rosenblatt Computation and interpretation of k th order spectra in Spectral Analysis of TimeSeries B Harris ed pp 189 232 J ohn Wiley 1967 Cadzow J A Spectral Estimation An Overdetermined Rational M ode
18. C ase Studies Table 1 2 Summary Statistics for Lynx Data Differenced Mean u E x t Variance c E x t u 2 Skewness E x t Yo Kurtosis E x t u yo 3 27 673 1 409 e 06 0 282 1 365 Routinehar mest indicates an order of three see panel 1 1 of Figure 1 36 and the power spectra are displayed in Figure 1 36 We see a peak at DC and two other peaks corresponding to a cycle of 9 63 years x 10 sval cum 2 2 NENNEN NI 0 0 10 20 30 40 pisar o o c ml 50 ar o P min norm 100 0 5 0 0 5 Figure 1 36 Lynx Data Power Spectra We used gl stat totest for Gaussianity and linearity the test results are Test statistic for Gaussianity is 196 752 with df 28 Pfa 0 Linearity test R estimated 6 8468 lambda 11 299 R theory 9 2282 N 5 1 115 1 Tutorial 1 116 The test indicates that the data are non Gaussian as expected recall that the skewness is around unity and that the histogram is indicative of an exponential distribution The time series length is only 114 samples and gl stat bases its estimate of the interquartile range on only 5 samples hence the test for linearity is not conclusive The bispectrum of the data was estimated via bi speci we used 25 lags and the default window The resulting estimate is shown in Figure 1 37 The bispectrum shows sharp peaks at aro
19. DOA 2 44 examples AR order determination 1 34 AR parameter estimation 1 32 ARMA parameter estimation 1 34 bicepstrum based IR estimation 1 39 1 41 bicoherence estimation 1 20 computing true cumulants 1 43 cross bicoherence estimation 1 20 cross bispectrum 1 19 cumulant esti mation 1 14 cumulation estimation 1 14 DOA estimation 1 77 Gaussianity linearity tests 1 24 harmonic retrieval 1 75 Levinson recursion 1 50 MA order determination 1 37 MA parameter esti mation 1 30 Matsuoka U Irych algorithm 1 42 QPC detection 1 87 RIV double lattice form 1 60 RIV transversal form 1 57 speech signal 1 122 sunspot data 1 108 time delay estimation 1 103 1 105 1 107 trench recursion 1 50 Volterra system identification 1 82 1 83 Wigner bispectrum 1 96 smoothed 1 97 Wigner spectrum 1 93 smoothed 1 94 Wigner trispectrum 1 99 smoothed 1 100 Index F FBLS 1 53 deterministic formulation 1 53 forward prediction problem 1 47 forward backward least squares problem 1 54 frequency coupling 1 128 frequency esti mation 1 65 2 50 G Gaussianity test 1 22 2 47 gl dat mat 1 135 gl stat 2 47 GM equations 1 29 guided tour 2 54 H harm mat 1 135 harmest 2 50 harmgen 2 53 harmonic retrieval 1 62 1 64 AR models 1 66 ARMA models 1 66 cumulant based method 1 74 minimum norm method 1 69 MUSIC 1 68 Pisarenko s method 1 67 synthetics 2 53 help 2 55 higher order spectra 1 2 higher order statistics 1 2 motivations 1 10 ho
20. Figure 1 5 Cross Bicoherence Estimate bicoherx 1 21 1 Tutorial 1 22 Testing for Linearity and Gaussianity n subsequent sections we discuss algorithms for the estimation of parameters of non Gaussian linear processes of course these routines assume that the data are linear and non Gaussian We call a process y n linear if it can be represented by y n Y h guin amp k where u n is assumed to be i i d If u n is Gaussian non Gaussian we say that y n is linear Gaussian non Gaussian How do we know that the data are non Gaussian and that they are additionally linear Hinich 24 has developed algorithms to test for non skewness loosely called Gaussianity and linearity The basic idea is that if the third order cumulants of a process are zero then its bispectrum is zero and hence its bicoherence is also zero If the bispectrum is not zero then the process is non Gaussian if the process is linear and non Gaussian then the bicoherenceis a nonzero constant see 1 19 1 20 and 1 14 with x y z Thus we have a hypothesis testing problem for non Gaussianity non zero bispectrum H1 the bispectrum of y n is nonzero HO the bispectrum of y n is zero If hypothesis H1 holds we can test for linearity that is we have a second hypothesis testing problem H1 the bicoherence of y n is not constant HO the bicoherence of y n is a constant If hypothesis HO holds the process is linear
21. Reference Detection of quadratic phase coupling using the TOR method arvec bspec qpctor y arvec bspec qpctor y maxlag arorder nfft samp seg overlap flag qpctor detects quadratically phase coupled harmonics using the TOR method y is the data matrix each column of y is assumed to correspond to a different realization maxl ag specifies the maximum number of third order cumulant lags C3 m m to be used ar order specifies the AR order to use lf this parameter is not specified or is not positive a plot of the singular values of the cumulant matrix is displayed and you are prompted to choose an order nfft specifies theFFT length its default valueis 64 samp seg specifies thenumber of samples per segment the default valueisthe length of the time series or the row dimension if y is a matrix overlap specifies the percentage overlap between segments maximum allowed value is 99 default valueis 0 the parameter is ignored if y is a matrix Ifflag is biased then biased sample estimates of cumulants are computed default if thefirst letter is not b unbiased estimates are computed arvec is the vector of estimated AR parameters bspec is the estimated parametric bispectrum It is an nf ft 2 by nfft array whose upper left handcornercorresponds to the origin in the bispectral plane the frequenciesincrease downwards and to the right Details of the algorithm are given in the Tutorial 1 Raghuvee
22. 1 2 e f wind is a 2 D matrix it is assumed to specify the 2 D smoother directly The bispectrum estimate averaged across records is smoothed by convolving with the 2 D window function The window function should be real and nonnegative samp seg specifies thenumber of samples per segment The default valueis set such that eight possibly overlapped records are obtained 2 27 bispecd Algorithm See Also References 2 28 overlap specifies the percentage overlap between segments The default value is 50 Theallowed range is 0 99 If y isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap is set to zero and samp seg is set tothe row dimension bspec isthe estimated bispectrum it isan nf ft by nfft array with origin at the center and axes pointing down and to the right waxis isthe set of frequencies associated with the bispectrum in bs pec Thus theith row or column of bs pec corresponds tothefrequency waxi s i i 1 nf f t Frequencies are normalized that is the sampling frequeng is assumed to be unity A contour plot of the magnitude of the estimated bispectrum is displayed The data y are segmented into possibly overlapping records the mean is removed from each record and the FFT computed the bispectrum of the kth record is computed as B m n Xy m X n X m n where Xy denotes the FFT of the kth record where denotes the FFT of the kth
23. 1 k 0 1 66 1 47 1 Tutorial 1 48 This leads to the normal equations 1 67 1 68 where 3i is the p 1 x p 1 autocorrelation matrix with m n entry Ry n m It follows immediately that f b c k a p k k 20 p P P 1 69 The Levinson Durbin recursion implemented in routinetrench provides an efficient order recursive solution of the normal equations it recursively solves 1 67 for p 1 2 The recursion is given by for m 1 p 32 m 1 Am 1 2 am 1 K Rxx K m k 0 1 70 Am 1 Tim A ue 1 71 Pin Pm 1l m m 1 m 1 72 Am k am 1 K I ma m 1 m k k 0 m 1 73 with initial conditions Pg R xx 0 Ag Rx 1 and ag 0 1 In practice we replace Rxx k by its biased sample estimate The resulting solution is linear Prediction M odels guaranteed to be stable if biased estimates are used Necessary and sufficient conditions for the stability of the estimate are I 1 Vip The recursion involves computation of theT s which are also called the reflection coefficients a necessary and sufficient condition for stability is Ty 1 Trench Recursion If the matrix 3i is Toeplitz but not Hermitian symmetric then 1 69 does not hold this for example is the case when the so called higher order Y ule Walker equations p L ayK Ry m K 0 me q L q p k 0 1 74 are used or when the cumulant based normal equations are used for examp
24. Identification of non minimum phase systems using higher order statistics IEEE Trans ASSP Vol 37 pp 360 77 Mar 1989 2 Swami A and J M Mendel Identifiability of the AR parameters of an ARMA process using cumulants IEEE Trans on Automatic Control Vol 37 pp 268 73 Feb 1992 2 13 armasyn Purpose Syntax Description 2 14 Generates ARMA synthetics zmat armasyn zmat armasyn default armasyn generates the time series described by z k E u k gt where u k is thei i d input to the ARMA p q model whose AR p and MA q polynomials are given by A z and B z respectively g k is signal independent noise generated via Ba k g k A w k where w k is an i i d sequence and the polynomials A z and B z determine the spectrum of the noise Ifar masyn isinvoked without any input arguments then you are prompted for all the parameters samples per realization number of realizations ARMA parameters for the signal process pdf of the input driving noise u k noise variance ARMA parameters for the observation noise process g k pdf of the noise w n signal to noise ratio If the function is invoked as ar masyn def aul t wherethe variabledef aul t may take on any value s then the default settings are used Time series y in filear 1 mat was generated via y r masyn 1 zmat returns the generated data each column corresponds to a different realization
25. L Time Frequency Distributions A Review Proc IEEE pp 941 81 J uly 1989 3 Choi H and WJ Williams Improved Time F requency Representation of Multicomponent Signals Using Exponential Kernels IEEE Trans ASSP pp 862 71 J une 1989 wig3 Purpose Syntax Description Algorithm See Also Computes the f4 f slice of the Wigner bispectrum wx waxis wig3 x nfft flag Computes the diagonal slice of the Wigner bispectrum WB x is the signal and must be a vector nfft specifies the FFT length and hence the frequency resolution this parameter must be greater than threetimes thelength of the signal x in order to avoid aliasing The default value is the power of 2 equal to or just greater than three times the signal length flag ifflag is nonzero and x is real valued the analytic form of x is used instead of x the default valueis 1 If x t is a real valued signal with Hilbert transform y t then theanalyticsignal is defined as the complex valued signal z t x t jy t wx is the computed WB The rows correspond toti me samples and the columns to frequencies waxis isa vector whose entries are the frequencies associated with the columns of the WB Define the instantaneous triple product r3 t tz t2 X t at at2 x t Bez atz x t at t r2 where a 1 3 and p 2 3 The WB is then given by j2 15 242 W n fis f2 fdzideze nin fat rc Ti T2 In this routine we comput
26. arorder Purpose Syntax Description Algorithm Estimates the AR order of an AR or ARMA process p arorder y norder pmax qmax flag Estimates the AR order of an AR or ARMA process using second third or fourth order cumulants y is the observed ARMA process and must be a vector norder specifies the cumulant order s to be used valid values are 2 3 and 4 a value of 3 4 indicates that order determination should be based on the correlation as well as the third order fourth order cumulants The default value is 3 pmax specifies the maximum expected AR order the default value is 10 qmax specifies the maximum expected MA order the default value is 10 flag iffl ag is1 theinternally chosen AR order is returned in p otherwise the plot of the singular values of the cumulant matrix is displayed on the graphics window and you are prompted to choosethe order The default value is 1 p is the estimated AR order Let p and q denote the maximum expected values of the AR and MA orders the parameters pmax and qmax respectively For convenience let c m p cum y n y n 4 m y n p y n y n k 22 note that we are suppressing k 3 of the lags all of which are set to zero Also let pee Then the singular values of the matrix C s are computed wherek 2 3 or 4 If norder is specified as 3 4 the matrices C5 and C3 C4 are concatenated 2 15 arorder See Als
27. e A unity value for wi nd results in no windowing If wind lt 0 the default value of 5 will be used If wi nd is a vector it is assumed to specify a 1 D window from which a 2 D window is computed W m n 2 w m w n w m n 1H2 If wi nd isa 2 D matrix it is assumed to specify the 2 D smoother directly The bispectrum estimate averaged across records is smoothed by convolving with the 2 D window function The window function should be real and nonnegative 2 79 tdeb Algorithm See Also Reference samp seg specifies the number of samples per segment The default valueisthe power of 2 equal to or greater than 4 max delay 1 overlap specifies the percentage overlap between segments The default value is 50 Theallowed range is 0 99 If y isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset tozero andsamp seg is set totherow dimension del ay is the estimated delay of signal y with respect to signal x h is the estimated third order hologram defined below it is a vector of length nf f t corresponding to indices of nfft 2 tonfft 2 1 wherenfft istheFFT length Define auto and cross bispectra via B ox 01 02 E K w1 X 2 X 1 02 3 Bxyx 91 02 E X 01 Y 2 X 1 4 92 The third order hologram h t is then defined by Bux 91 05 yx r 2 h t do do epdor 5 cea x The absolute value of the hologram should display a strong peak at the loca
28. m o T L L L L L L L 0 25 0 2 0 15 0 1 0 05 0 0 05 0 1 0 15 0 2 frequency Figure 1 21 Wigner Spectrum With Choi Williams Smoothing wig2c Wigner Bispectrum The notion of Wigner spectra can be generalized to Wigner bispectra WB and Wigner trispectra WT by considering the FTs of appropriately defined instantaneous triple and quadruple products Time Frequency Distributions Definethe instantaneous triple product r2 t T4 05 x t at at5 x t Bt at5 x t at Bt 3 1 2 1 2 1 2 1 2 1 180 where a 1 3 and B 2 3 The WB is then given by Won fz f2 d 9 67 tt t T2 1 181 This algorithm is implemented in routine wi g3 where the slice f4 f gt is computed Notethat theoriginal signal should be sampled at twicethe Nyquist rate in order to avoid aliasing It is shown in that the frequency axes are scaled by the factor of 2 3 which can be easily fixed Asinthecaseof the WS the WB also suffers from the problems of cross terms one can attempt to suppress these terms by appropriate filtering Define the smoothing kernel 2 2 2 0 t t5 exp 0 0 Ty T2 p 0 1 15 0 1 182 The filtered WB is then given by Win fy f de4d due 177105 feng anth g nuo r3 t Ty 15 P G 14 T2 This algorithm is implemented in routine wi g3c where the slice f f2 is computed It should also be noted that the application of the Choi Williams filter to the WB d
29. real harmonics and the spectrum appears to peak at the correct frequencies 0 1 and 0 2 Hz as well as at 0 15 Hz which is due to the Gaussian noise sval cum 2 pxx 10 0 3 R PIN 0 40 0 5 10 15 0 5 0 0 5 0 0 R o D 2 20 5 pisar eo Oc mi S eo C zl E i 5 0 0 5 0 5 0 0 5 0 0 E S 50 T E 100 50 0 5 0 0 5 0 5 0 0 5 Figure 1 13 Singular Values of Correlation Matrix and Estimated Spectra harmest 1 76 Harmonic Processes and DO A Examples load doal spec theta doa ymat 0 5 1 3 2 1 You should see the display in Figure 1 14 note that the angular spectra are plotted in dB scale that is 10logjo spec The singular value plot indicates the possible presence of three sources The observed signal consists of signals from two sources at 15 and 25 degrees and is contaminated with spatially correlated noise that acts as a virtual source at a bearing of 30 degrees The bearings estimated by ESPRIT should be 15 4634 25 8851 and 29 9855 As expected the Eigenvector ML AR MUSIC and Min Norm methods resolvethetwo sources and the virtual source duetothe noise whereas the beamformer does not TheMUSIC and Min Norm estimates are virtually identical in this example The vertical lines and circles denote the true bearings the cross denotes the virtual bearing of the spatially correlated noise
30. standardize fi gure 1 data and histograms subplot 211 plot l length sp sp grid subplot 212 hist sp grid Summary stats cl mean sp c2 cumest sp 2 c3 cumest sp 3 c2 31 2 c4 cumest sp 4 c2 2 fprintf Mean g n cl fprintf Variance g n C2 fprintf Skewness normalized gin c3 fprintf Kurtosis normalized g n c4 figure 2 specgramsp 512 8000 hammi ng 256 240 oorr II power spectra and cum 4 spectra figure 3 px2 a21 a22 harmest sp 25 12 biased 512 2 figure 4 px4 a41 a42 harmest sp 25 8 biased 512 4 r21 roots a21 r22 roots a22 r4l roots a41 r42 roots a42 figure 5 subplot 221 polar angle r21 abs r21 x grid subplot 222 polar angle r22 abs r22 x grid subplot 223 polar angle r41 abs r41 x grid subplot 224 polar angle r42 abs r42 x grid gistat sp 0 51 256 gl test orm III bispectrum and diagonal slice figure 6 Blaf w bispecd sp 256 1 100 0 figure 7 dB abs diag Blaf plot w dB grid title diagonal slice loc val pickpeak dB 3 disp w loc figure 8 dqpctor sp 25 10 512 100 0 QPC test C ase Studies Pitfalls and Tricks of the Trade As the previous examples illustrate we should first use univariate statistics mean median skewness kurtosis histograms and second order statistics SOS Beforeusingalgorithms based on higher or
31. the true values of the cumulant C3y m 0 will beidentically zero if m gt q When thetrue cumulants are replaced by sample estimates the estimated values of c3 q 1 0 i 20 will not be identically zero a statistical test is used to determine whether the estimated values are close to zero The test is based on estimating the theoretical variance of the sample estimates of c3 m 0 Sample esti mates C3y m 0 and their variances are estimated for m ranging from qmi n toqmax For an MA q process the asymptotic variance of the sample estimate of cs q 1 0 is given by E i UN Y yt q 2 6y q 1 0 3 1 N 2q o q 1 3 N i 3j 1 q Y DyG j q 1 y q 1 0 2 61 maorder See Also Reference 2 62 where N is thelength of thetime series The sample estimates are asymptotically normal and unbiased hence the threshold t in Pr esy m 1 0 t m 1 1 pfa is given by t m 1 inverf 1 pfa 42 o m 1 whereer fi nv isthe MATLAB inverseerror function Let mg denote the largest value of m in the rangeqmi n toqmax for which csy m 1 0 gt t m 1 so that the hypothesis of M A mo model fails Then the MA order is declared to be mg 1 In the example mg is 2 and the MA order was declared to be 3 This is a statistical test and pf a specifies the fraction of the time that the test results will be wrong In other words in a Monte Carlo simulation of 1000 trials
32. theory 9 1222 N 16 The probability of false alarm Pfa in rejecting the Gaussian hypothesis is 0 that is we are virtually certain that the data have nonzero bispectrum thus the test indicates that the data are non Gaussian as expected recall that the histogram shows that the univariate pdf is non Gaussian and that the skewness is around unity The estimated and theoretical values of R the interquartile range of the estimated bicoherence values are not very different from one another hence the tests do not show any strong evidence of nonlinearity The bispectrum of the data was estimated via bi speci we used 25 lags and the default window The resulting estimate is shown in Figure 1 31 The bispectrum shows sharp peaks at around 0 0 1 and 0 1 0 1 and symmetric locations and is indicative of possible quadratic coupling The peak at 0 0 1 is most probably an artifact due to the fact that the data areall positive valued Bispectrum estimated via the indirect method T T T T l 1 l I l i l 1 0 5 0 4 0 3 0 2 0 1 o 0 1 0 2 0 3 0 4 Figure 1 31 Sunspot Data Bispecrum C ase Studies Differencing the data is a standard trick in exploratory data analysis we can use MATLAB sdi f f function to compute the first order differences The differenced data and the corresponding histogram are shown in Figure 1 32 Summary statistics are shown in thelast column of Table 1 1 150 T 0 50
33. unbiased estimates are computed only the first letter of f ag is checked Cumulants are estimated from each record and then averaged across the set of records nfft specifies the FFT sizeto be used for computing the bispectrum and the bicepstrum longer FFT lengths yield better estimates but also require more storage two arrays of sizenfft by nf tt The default valueis 128 if nf ft is smaller than 2 n ag 41 the power of 2 just larger than 2 nl ag 1 will be used wi nd specifies the lag domain smoothing window If wi nd is 0 the Parzen window will be applied Otherwise the usual unit hexagonal window will be applied The Parzen window is the default 2 21 bicepsf Algorithm See Also Reference 2 22 The 1 D Parzen window is defined by 2 3 ps I I Im L 2 L L dy m AG Im L 2x m xL a Im L whereL nl ag Theactual window applied tothe estimated cumulants is given by W m n dy m dy n dy m n The unit hexagonal window is given by W m n d m d n d m n where d m 1 m nl ag hest is the estimated impulse response h n n 2 nfft 2 nfft 2 1 Samples of h n n lt 0 will have significant values if the original linear system is not minimum phase e g the original system is an ARMA model some of whose zeros or poles lie outside the unit circle ceps is the estimated complex cepstrum h n n nfft 2 nfft 2 1 Note that the method has an inherent
34. you should expect the test results to be wrong 1000 pf a times Additionally note that the test condition is necessary but not sufficient arorder 1 Giannakis G B andJ M Mendel Cumulant based order determination of non Gaussian ARMA models IEEE Trans ASSP pp 1411 21 Aug 1990 matul Purpose Syntax Description Algorithm See Also References Nonparametric magnitude and phase retrieval using the M atsuoka U Irych algorithm hest matul bspec The phase and log magnitude of the transfer function are estimated via the Matsuoka UIrych algorithm and then converted to the time domain impulse response bspec isthe bispectrum array such as that computed by bi speci orbispecd hest is the estimated impulse response The phase unwrapping algorithm in 2 is used to resolve the phase ambiguity in the basic algorithm reported in 1 bispecd bispeci biceps 1 Matsuoka T and T J Ulrych Phase estimation using the bispectrum Proc IEEE Vol 72 pp 1403 11 1984 2 Rangoussi M and G B Giannakis FIR modeling using log bispectra Weighted least squares algorithms and performance analysis IEEE Trans Cir Sys Vol 38 pp 281 96 1991 2 63 nlgen Purpose Syntax Description See Also Computes the output of a second order Volterra system y nlgen x h q Computes the output of a second order Volterra system x is the input to the nonlinear system it may be a
35. 10 0 4 10 0 0 2 0 4 0 6 0 4 0 2 0 0 2 0 4 Figure 1 51 Transfer Functions of a Second Order Volterra Model nl2 mat This file contains data to test function nl pow This data file was generated in the same way as was nl 1 mat except that the process x was chosen to be i i d and Laplacian Data Files qpc mat A quadratically phase coupled synthetic The data consist of four harmonics corrupted by white Gaussian noise The frequencies of the harmonics are A4 0 1 Az 0 15 A3 20 25 and Aq 0 40 Hz For each realization the phases of the first second and fourth harmonics 4 62 and 4 were chosen randomly the phase of the third harmonic was set to 5 2 Note that the triplet of harmonics with frequencies 0 1 0 15 0 25 exhibits both phase and frequency coupling the triplet with frequencies 0 15 0 25 0 40 is frequency coupled but not phase coupled Sixty four independent realizations each consisting of 64 samples were generated The amplitudes of all four harmonics wereunity and white Gaussian noise with a variance of 1 5 was added tothesignal Theset of realizationsis contained in the matrixz mat the columns correspond to independent realizations riv mat This data file contains three vectors y zw andzc each consisting of 1024 samples A zero mean exponentially distributed i i d sequence was passed through the AR filter 1 1 5 0 8 to obtain y Additive white Gaussian noise was added to y to obtain
36. 100 150 200 250 300 0 L 11 L L 1 f L 60 40 20 0 20 40 60 80 100 120 Figure 1 32 Sunspot Data First Difference 1 111 1 Tutorial The singular value plot displayed by har mest see the 1 1 panel in Figure 1 33 indicates either four or six significant singular values using p 6 we obtain the power spectra shown in the remaining panels of Figure 1 33 Again using roots we estimate the periods to be 5 47 and 10 50 years AR method We used gl stat totest for Gaussianity and linearity the test results are Test statistic for Gaussianity is 250 1965 with df 70 Pfa 0 Linearity test R estimated 13 5335 lambda 6 4449 R theory 7 0433 N 16 The test indicates that the data are non Gaussian as expected recall that the skewness is around unity sval cum 2 pxx 10000 0 0 50 0 10 20 30 40 0 5 0 0 5 0 0 2 o 20 2 20 40 40 0 5 0 0 5 0 5 0 0 5 0 0 S 50 WP 20 Poul Vw o 100 40 0 0 0 5 0 5 0 0 5 0 0 E o 50 T 20 E 4 100 40 0 5 0 0 5 0 5 0 0 5 Figure 1 33 Differenced Sunspot Data Power Spectra 1 112 C ase Studies The bispectrum of the differenced data in Figure 1 34 has sharp peaks around 094 094 and indicates coupling between the periods at 5 3 years and 10 6 years Bispectrum estimated via the indirect
37. 12 The discrete time algorithm is given next 13 25 Let the instantaneous cross correlation be defined by ryy m n x n m y n m 1 176 wheren is identified with time and m with lag The WCS is then defined by 1 91 1 Tutorial 1 92 Wyy f n Yr n m expojzfm a 1 177 The WS is obtained when x n y n The original signal must be sampled at twice the Nyquist rate or faster in order to avoid aliasing In practice the frequency variable f is also discretized f k K whereK controls the frequency resolution The WCS in 1 177 can also be implemented via two FFT algorithms following the approach in 46 Both approaches demand the same computational and storage complexity This algorithm in 1 177 is implemented in routine wi g2 The ambiguity function AF can be computed as AF m 8 Y ryan n exp j2xn9 9 1 178 The AF is multiplied by the Choi Williams filter 2 w m 0 exp m0 0 1 179 A 2 D FT 8 ton and mto f yields the filtered WS Cross terms can usually be suppressed via this approach but with a concomitant loss of resolution In practice 0 is a discretized frequency grid This algorithm is implemented in routine wi g2c Time Frequency Distributions Examples cl f load wigdat Subplot 221 wig2 s2 subplot 222 wig2 s2 subplot 223 wig2 s3 Subplot 224 wig2 s4 You should see the display in Figure 1 20 ws WS 60 60 u 50 o 50 E o
38. 124 liein the signal subspace and are orthogonal to the noise subspace that is ey v 0 i p l M k 1 p hence with ange j2nf j2xM f T ef d g e TE one can search for the set of frequencies such that e f v 20 k p 1 M When sample estimates of the autocorrelation are used these relationships will hold only approximately In MUSIC Multiple Signal Classification one determines the frequencies of the harmonics by looking for peaks in the spectrum defined by M H 2 1 PC Y wwe wv MED 1 138 where w k 1 Usually the correlation matrix used in the F BLS method is used When w k 1 A we get the so called Eigenvector solution and when W k 8 k M we get the Pisarenko solution These algorithms are implemented in routines har mest and doa Harmonic Processes and DO A Minimum Norm Method The minimum norm M in N orm method is also an eigen space method Motivated by the normal equations for the harmonic signal model one seeks the AR vector a which is orthogonal to the signal subspace recall the discussion of the Pisarenko method L et us decompose the signal subspace as T Qs Vy Vpl E X 1 139 The condition Qla 0 can be rewritten as T X w x z 1 140 where w a 1 a p Equation 1 140 represents p equations in M 1 unknowns Hence if p M 1 the set of equations is underdetermined and the solution is not unique A unique solution is
39. 2 m M where M isthe FFT length We can collect 1 158 at the various frequencies into a set of equations that are linear in the unknowns namely H f and Q f4 f2 50 N onlinear Processes Let r m mod 2 and let m m r 2 m m r2 Form 20 M 2 let the unknown parameter vector be b m Hom Q m m Q m 1 m 1 qs m A and A Z SF X oN 3 gt Z lt A X m X m X M Xm DX M 2 x Then we have the system of linear equations Y m 2 A b m which is overdetermined we can obtain the least squares solution for each frequency grid point m This algorithm is implemented in routinen pow Examples load nll hl ql nl pow x y 128 You should see the display on Figure 1 17 linear TF quadratic TF 10 z 0 2 107 O 1 107 o 10 O lt 4 0 2 10 oO 0 2 0 4 0 6 0 2 O 1 o O 1 0 2 T fi linear part IR quadratic part IR 0 4 sol Fae Zi 0 3 25 CD E 0 2 20 eu A YE isl LP 7 15 o M TT 5 EM dA LO A ae 50 100 150 10 20 30 t1 t1 Figure 1 17 Transfer Functions Estimated by nlpow 1 83 1 Tutorial 1 84 load ni2 h2 q2 nl pow x y 128 You should see the display on Figure 1 18 linear TF linear part IR OW Wee o 50 100 ti 150 quadratic TF Figure 1 18 Transfer Functions Estimated by nlpow As expected
40. 4f We can test whether the bispectrum of the data are statistically nonzero by usinggl stat Thetest results usingc par m 0 51 nf f t 256 were Test statistic for Gaussianity is 71 3231 with df 48 Pfa 0 Linearity test R estimated 2 3216 lambda 2 376 R theory 4 472 N 14 Since the Pfa is close to 0 we can be virtually certain that the data have nonzero bispectrum and hence are non Gaussian The estimated and theoretical R values are not close to each other indicating that the data are not linear 1 125 1 Tutorial sval cum 4 o 50 eig eo 0 5 a pisar E 400 0 5 a 3 1 s 100 0 5 0 0 Figure 1 46 a 0 5 d 0 5 7 4 ar 0 5 J 0 5 J Figure 1 47 Laughter Data 1 126 5 pxx mi l m e 40 0 5 min norm 50 0 5 Laughter Data Cumulant Spectra eK Ky gt pe PEE x x Root Locations 0 5 C ase Studies In order to check for quadratic frequency coupling we can estimate the nonparametric bispectrum either via bi speci orbispeci or the parametric bispectrum via qpct or The bispectrum estimated via bi speci with nf f t 256 s e gs a mp 100 over ap O wi nd L is shown in Figure 1 48 Sincethe data contain several harmonic components the bispectrum displays several impulses in the case of speech signals w
41. 7200 0 0 0 0 0 0000 The backward AR coefficient matrix c 4 should be 0 0 0 0 0 0045 0 0 0 0 2612 0 2599 0 0 0 1620 0 1141 0 1137 0 0 0786 0 0454 0 0278 0 0261 0 1956 0 1776 0 1588 0 0868 0 0900 1 0000 1 0000 1 0000 1 0000 1 0000 and the prediction error variances p4 should be 0 8892 0 8810 0 8417 1 0000 1 0000 The forward reflection coefficients gf 4 are 0 2278 0 1163 0 2757 0 7200 0 0000 and the backward reflection coefficients gb4 are 0 1956 0 0786 0 1620 0 2612 0 0045 Since the Toeplitz matrix is not symmetric gf 4zgb4 andaz4fli pud c4 Note alsothat the prediction error variances may increase with the prediction order One can find the AR parameters by minimizing the variance of the AR compensated process but not by minimizing the skewness Although the correlation and cumulant sequences are derived from the same AR model the lower order AR model fits are quite different 1 52 linear Prediction M odels Deterministic Formulation of FBLS In the deterministic formulation of the linear prediction problem rather than minimizethe varianceof the prediction errors weminimizethe sum of squared errors that is N P Y jew n p 1 1 85 N P Y lem n p 1 1 86 This leads to the deterministic normal equation 1 87 c 9 p Pe 1 88 where the deterministic correlation matrix is given by N o y x n x n TEES 1 89 x n x n x n 1
42. A and J M Mendel Cumulant based approach to the harmonic retrieval and related problems IEEE Trans ASSP Vol 39 pp 1099 1109 May 1991 see also Proc CASSP 88 pp 2264 67 1988 3 Cadzow J A Spectral Estimation An Overdetermined Rational M odel Equation Approach Proc IEEE Vol 70 pp 907 38 1982 4 Haykin S AdaptiveFilter Theory New J ersey Prentice Hall pp 464 70 2nd ed 1991 5 Kumaresan R and D W Tufts Estimatingtheangles of arrival of multiple plane waves IEEE Trans AES Vol 19 pp 134 39 1983 harmgen Purpose Syntax Description Generates harmonics in colored Gaussian noise zmat harmgen z mat harmgen default har mgen generates independent realizations of the signal p y n Y ak COS 21 n 4 N k 1 Ifhar mgen is invoked without any input arguments then you are prompted for the length of the realizations the number of realizations the number of harmonics p their frequencies and amplitudes oj and the variance of the additive colored Gaussian noise g n The additive colored Gaussian noise g n is generated by passing a white Gaussian noise sequence through an user specified ARMA filter you are prompted for these ARMA parameters For each realization the phases are chosen randomly from an uniform distribution Note that noise free realizations can be obtained by specifying a valueof zerofor the noise variance when prompted I
43. Function Purpose doa Estimates DOA from a linear array of sensors based on spatial cross cumulant or covariance matrix doagen Generates synthetics for the DOA problem 24 Function Tables Adaptive Linear Prediction Function Purpose i vcal Computes instrumental variable processes rivdl Adaptive LP using double lattice filter rivtr Adaptive LP using transversal filter Impulse Response IR Magnitude and Phase Retrieval Function Purpose biceps Estimates IR complex cepstrum lag domain bicepsf Estimates IR amp complex cepstrum FFT method mat ul Estimates Fourier phase and magnitude of a signal using the Matsuoka Ulrych algorithm Time Frequency Estimates Function Purpose wi g2 Wigner spectrum wi g2c Wigner spectrum with Choi Williams filtering wi g3 Wigner bispectrum wi g3c Wigner bispectrum with Choi Williams filtering wi g4 Wigner trispectrum wi g4c Wigner trispectrum with Choi Williams filtering 2 Reference Utilities Function Purpose hosahelp hprony pickpeak tls trench One line help of Higher Order Sectral Analysis Toolbox commands Modeling of transient signals via Prony s method Picks peaks subject to a separation criterion Total Least Squares solution Trench recursion for non symmetric Toeplitz matrix Demo Function Purpose hosademo Guided tour of the Higher Order Sectral Analysis Toolbox 2 6 M i
44. Invariance Techniques ESPRIT is an eigen space method developed by Kung Rao amp al 30 and by Roy and Kailath 57 however the basic idea behind this algorithm is considerably different from that in MUSIC and Min Norm methods In ESPRIT the basic idea is to decompose the sensor array into two possibly overlapping subarrays and then to usethe cross correlation between the two subarrays to esti mate the source bearings The discussion here is in terms of the DOA problem but as we have seen earlier the treatment is applicable to the harmonic retrieval problem as well Assume that we have two arrays such that elements of one array are at a constant displacement A with respect to corresponding elements of the other array For example if we have an uniformly spaced linear array of N elements we can partition this array into two subarrays each consisting of L elements with a displacement of A N L between the two subarrays that is with sensors 1 N A in thefirst subarray and sensors A 1 N in the second array Let ym N and Z n n 21 L denote the signals at the elements of the two subarrays that is p Ym n Y a m exp jo k n wg k l 1 143 j2nAd sin o A z m Y o m exp jock npe Vm N k 1 1 144 Harmonic Processes and DO A The mth snapshot can be written in vector form as Ym Aag Wm 1 145 Z A a V i pi 1 146 where is a diagonal matrix with k k entr
45. References The filtered WT is then given by Won fy f f3 fff dodridrzdrzdue 1770 m tf feo j2nt0 j2nu0 el d T ra t t t5 T 0 tq Ta T3 In this routine we compute the slice f4 f f3 Note that the original signal should be sampled at twice the Nyquist rate in order to avoid aliasing Also note that the frequency axes are scaled by the factor of 1 2 In this routine we undo the scaling so that you can find the peaks at the expected frequencies It should be noted that the Wigner trispectrum is real valued wig2 wig2c wig3 wig3c wig4 1 Choi H and WJ Williams Improved Time F requency Representation of Multicomponent Signals Using Exponential Kernels IEEE Trans ASSP pp 862 71 J une 1989 2 F onollosa J R and C L Nikias Analysis of finite energy signals using higher order moments and spectra based time frequency distributions Signal Processing Vol 36 pp 315 28 1994 2 97 2 98 A adaptive filter double lattice 2 72 RIV 2 74 adaptive linear prediction 1 54 ambiguity function 1 90 1 92 AR method DOA 2 44 harmonic retrieval 2 50 AR models 1 31 order determination 1 34 2 15 parameter estimation 2 17 parameter identifiability 1 31 arl mat 1 134 ARMA models 1 32 AR order estimation 2 15 AR parameter estimation 2 17 residual time series 2 12 armal mat 1 135 armags 1 33 2 8 armarts 2 11 armasyn 2 14 arorder 2 15 arrcest 2 17 autocorrelation 1 6 B backward predi
46. Y C m k 1 im cum3x cum4x cumest 2 35 cum3x Purpose Syntax Description 2 36 Computes the third order cross cumulant of three signals cvec cumix x y z maxlag samp seg overlap flag kl Computes the third order cross cumulant of thethreesignals x y and z which should have identical dimensions maxlag specifies the maximum lag of the cumulant to be computed its default valueis O samp seg specifies the number of samples per segment Its default value is the row dimension of y if y is a row vector the column dimension is used as the default value overlap Spedifies the percentage overlap between segments the allowed range is 0 99 and the default value is O Iff I ag is biased then biased sample estimates are computed default if the first letter is not b unbiased estimates are computed Parameter k1 controls which 1 D slice of the cross cumulant is computed the default value is 0 By varying k1 we can obtain the entire third order cross cumulants cvec will contain the sample estimates of E x n ug Y n m py z n k pz m maxl ag 2 s maxlag Here u denotes the mean of process x and the superscript denotes complex conjugation If x y z are matrices columns are assumed to correspond to different realizations in this caseoverlap isset to zero andsamp seg is set tothe row dimension the cross cumulant is estimated from each realization and then averaged ac
47. are relatively new As such there is no guarantee that a particular Higher Order Spectral Analysis Toolbox routine will work well on your data United Signals Systems Inc will be updating and upgrading the Higher Order Spectral Analysis Toolbox from time to time to incorporate new routines and provide users with guidance on the applicability of existing routines as more experience is obtained through their use The Tutorial has numerous examples that reinforce the theory and demonstrate how to use the toolbox functions All of the data files used by these examples are included in your Higher Order Spectral Analysis Toolbox distribution diskette and are described in the section on Data Files We encourage you totry out the examples yourself Additional examples may be found in the demo which can be invoked via hosademo A later section in the Tutorial demonstrates how to deal with real data 1 Tutorial Polyspectra and Linear Processes 1 4 In this section we will define cumulants polyspectra and various other related statistics such as bicepstra and bi coherence We will discuss tests for linearity and Gaussianity and we will develop cumulant based algorithms for estimating the parameters of linear e g ARM A processes Introduction The notion of decomposing a signal intoits harmonic components dates back to the analysis of the motion of planets the music of the spheres as the Pythagorians called it phases of th
48. be used with great caution Wi g2 wi g2c wi g4 wi g4c 1 Choi H and WJ Williams Improved Time F requency Representation of Multicomponent Signals Using Exponential Kernels IEEE Trans ASSP pp 862 71 J une 1989 2 Fonollosa J R and C L Nikias Wigner Higher Order Moment Spectra Definitions Properties Computation and Applications to Transient Signal Detection IEEE Trans SP J an 1993 2 93 wig4 Purpose Syntax Description Algorithm Computes the f f f3 slice of the Wigner trispectrum wx waxis wig4 x nfft flag Estimates the f f f3 slice of the Wigner trispectrum WT of a signal using the ti me domain approach x is the signal and must be a vector nfft specifies the FFT length and hence the frequency resolution this parameter must be greater than four times the length of the signal x in order to avoid aliasing The default value is the power of 2 equal to or just greater than four times the signal length flag ifflag is nonzero and x is real valued the analytic form of x is used instead of x the default valueis 1 If x t is a real valued signal with Hilbert transform y t then theanalyticsignal is defined as the complex valued signal z t x t jy t wx isthe computed WT The rows correspond to time samples and the columns to frequencies waxis isa vector whose entries arethe frequencies associated with the columns of the WT For a complex
49. demanding that the prediction error sequence be orthogonal to an instrumental process derived from the data 59 64 We put normal in quotes to emphasize these differences Identifiability of the AR parameters is guaranteed by choosing p q p andm q 1 q p we may use more slices p or more lags m 67 In practice we use sample estimates of the cumulants This algorithm is implemented in routinear marts The pureAR case corresponds to q 0 and is implemented in routinearrcest In both routines you can use cumulant orders 2 3 or 4 It is also possible to simultaneously solve for the normal equations based on cumulant orders 2 and 3 or 2 and 4 1 31 1 Tutorial Examples Try the following load arl ar l sarrcest y 2 0 2 12 128 ar 2 arrcest y 2 0 3 12 128 ar 3 arrcest y 2 0 4 12 128 ar 4 sarrcest y 2 0 3 12 128 ar 5 arrcest y 2 0 4 12 128 disp ar 1 0000 1 0000 1 0000 1 0000 1 0000 1 4636 1 5559 1 4963 1 4912 1 4755 0 7664 0 8779 0 8917 0 7973 0 7927 The true parameters are 1 1 5 0 8 The five columns correspond to AR esti mates based on 1 second order 2 third order 3 fourth order 4 combined second order and third order and 5 combined second order and fourth order cumulants N ote that combined use of autocorrelation and cumulants may give better results when thesignal to noise ratio SNR is high In the case of low or moderate SNR the correla
50. discrete time stationary random process can be expressed in the form x n y n z n such that 1 Processes y n and z n are uncorrelated with one another 2 Process y n has a causal linear process representation y n Y h k u n k k 0 where h 0 1 hee gh Ado lt and u n is a white noise process and 3 z n issingular that is it can be predicted perfectly with zero variance from its past An example of a singular process is the harmonic process s n aexp j2zfn A process with z n 0 has a purely continuous spectrum additionally a strictly band limited process is also singular Since real world signals cannot be strictly band limited we may think of the Wold decomposition as decomposing a process into a linear process which has a continuous spectrum and a harmonic process which has a line spectrum It is also important to note that the theorem only states that u t is uncorrelated it does not state that u t is i i d higher order white For example u n might be the output of an all pass system whose input is an i i d process We need higher order statistics to determine whether or not u t is i i d or merely uncorrelated Other motivations for using higher order statistics HOS are discussed throughout this Tutorial 1 Tutorial 1 6 Definitions The autocorrd ation function or sequence of a stationary process x n is defined by Ry m E x n x n m 1 1 where E
51. impulse response is then estimated via P a k C3 q k n scheme CRIS ED inh at or via p palk Cay q k n 0 99 0 C4 q k 0 0 These equations can be readily verified by using 1 16 through 1 18 and 1 35 The MA parameters are then obtained via 1 35 which is repeated here p b n be a k h n k n 1 q k 0 1 49 1 33 1 Tutorial 1 34 An interesting point is that we can simultaneously solve for the AR and IR parameters This algorithm is implemented in routinear mags Weighted versions of 1 42 through 1 44 and nonlinear cumulant matching algorithms for the simultaneous estimation of AR and MA parameters are discussed in 16 and 71 Examples load armal avec bvec z rmaqs y 2 1 3 10 128 Here we used third order cumulants and the q slice algorithm to estimate the parameters of a non Gaussian ARMA process The estimated parameters should beavec 1 0 8057 0 6910 andbvec 1 1 9116 Thetrue AR and MA parameters were 1 0 8 0 65 and 1 2 respectively avec bvec armarts y 2 1 3 12 128 Here we used third order cumulants and the residual time series algorithm to estimate the parameters of a non Gaussian ARMA process The estimated parameters should beavec 1 0 7540 0 6465 andbvec 1 1 5665 The true AR and MA parameters were 1 0 8 0 65 and 1 2 respectively AR Order Determination Let p and q denote the maximum expected va
52. m n 1 2 e f wind isa2 D matrix it is assumed to specify the 2 D smoother directly The bispectrum estimate averaged across records is smoothed by convolving with the 2 D window function samp seg specifies thenumber of samples per segment The default valueis set such that eight possibly overlapped records are obtained overlap specifies the percentage overlap between segments The default value is 50 The allowed range is 0 99 If y is a matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset tozero andsamp seg is set to the row dimension 2 29 bispecdx Algorithm See Also References 2 30 flag acontour plot of the estimated cross bispectrum will be displayed only iffl ag is nonzero the default value is 1 bspec is the estimated cross bispectrum It is annfft by nfft array with origin at the center and axes pointing down and to the right waxis is the set of frequencies associated with the cross bispectrum in bs pec thus the ith row or column of bs pec corresponds to the frequency waxis i i l nf ft Frequencies are normalized that is the sampling frequency is assumed to be unity The cross bispectrum definition used in this routine is given by Byz 01 Q5 E X 04 Y 05 Z 04 2 and is the 2 D Fourier transform of the cross cumulant defined by Cyyz m n E x t myy t n z t j For a complex process the cross cumulant m
53. method T T T T T T O 4r s 0 3 i d j i 4 0 2 4 Figure 1 34 Differenced Sunspot Data Bispectrum In this example differencing the data helped to clarify the estimates of the power spectrum as well as the bispectrum Our tests indicatethat the data are non Gaussian and show evidence of a fundamental period of around 10 6 years as well as a harmonic at around 5 3 years this may beindicative of quadratic nonlinearities The various figures were generated via load sunspot dat sp sunspot 2 eda sp eda diff sp Function eda is listed in Table 1 3 1 113 1 Tutorial Canadian Lynx Data Examples This time series consists of the annual number of Canadian lynx trapped in the Mackenzie River district of Northwest Canada for the years 1821 1934 listings of the data may be found in 56 and 52 where this time series is discussed in detail and additional references are also given We used eda to generate Figure 1 35 through Figure 1 40 The histogram in Figure 1 35 is indicative of an exponential distribution and the data show apparent periodicity Summary statistics are shown in Table 1 2 the last column corresponds to first differences 8000 6000 4 4000 2000 0 L L 0 20 40 60 80 100 120 60 l 40 4 20r y A 0 I I Cc a 0 1000 2000 3000 4000 5000 6000 7000 Figure 1 35 Lynx Data 1 114
54. obtained by choosing the Min Norm solution that is 1 x xq X w 1 Xx X x S s s s 1 141 This solution called the minimum norm solution was developed by Kumaresan and Tufts and is implemented in routines har mest and doa Pisarenko s method uses a single noise subspace eigenvector and hence is not as robust as the general MUSIC estimator The Min Norm method has been reported to have a smaller bias than the MUSIC algorithm for estimating the frequencies An overview of eigen methods is given in 26 A modification of the FBLS algorithm has been reported by Kumaresan and Tufts 23 where the correlation matrix is replaced by throwing out the noise subspace eigenvectors that is if A and vy denote the eigenvalues and eigenvectors of 9 then p x H Y ANKE kz 1 142 1 69 1 Tutorial 1 70 This procedure of replacing by is sometimes called the low rank SVD approximation The algorithm then proceeds as in FBLS Kumaresan and Tufts 29 have established that in the noiseless case under the assumption p lt M lt N p 2 p of the roots of A z lie on the unit circle and correspond to the true frequencies the remaining M p roots are uniformly distributed in angle inside the unit circle these roots are called the extraneous zeros High resolution estimates are obtained by choosing large values for the predictor order that is M 3N 4 ESPRIT Estimation of Signal Parameters via Rotational
55. orders 1 through M the kth column corresponds to the AR k model cmat is the matrix of estimated backward AR prediction vectors for orders 1 through M the kth column corresponds to the AR Kk model pf isthe final prediction error variance for orders 1 through M gamf is the vector of forward reflection coefficients gamb is the vector of backward reflection coefficients The standard Trench algorithm 2 is implemented the Levinson Durbin algorithm 1 is obtained whenr c 1 Haykin S Adaptive Filter Theory New J ersey Prentice Hall pp 198 205 2nd ed 1991 2 Zohar S Toeplitz Matrix Inversion the Algorithm of W F Trench J Assoc Comp Mach Vol 16 pp 592 601 1968 2 85 trispect Purpose Syntax Description Algorithm See Also References 2 86 Computes a 2 D slice of the theoretical trispectrum of an ARMA process tspec waxis trispect ma ar nfft f3 A 2 D sliceof thetheoretical trispectrum corresponding to an ARMA process is computed ma is the MA parameter vector and must be specified ar isthe AR parameter vector its default valueis 1 nfft istheFFT length touse its default value is 512 f 3 specifies the fixed value of the third frequency f3 of thetrispectrum S3 f1 f2 f3 The default value is 0 the nominal range is 0 5 0 5 values outside this range are folded back into it tspec is the 2 D slice of the trispectrum corresponding to the ARMA model w
56. record The bispectral estimates are averaged across records and an optional frequency domain smoother specified by parameter wi nd is applied bispeci 1 Subba Rao T and M Gabr An Introduction to Bispectral Analysis and Bilinear TimeSeries Models pp 42 43 New York Springer Verlag 1984 2 Nikias C L and M R Raghuveer Bispectrum estimation A digital signal processing framework Proc IEEE Vol 75 pp 869 91 J uly 1987 bispecdx Purpose Syntax Description Cross bispectrum estimation using the direct FF T method bspec waxis bispecdx x y z nfft wind samp seg overlap flag The cross bispectrum of the three processes x y andz Byyz 1 2 is estimated via the direct FFT method x y andz should havethe same dimensions nfft specifies the FFT length to be used for computing the cross bispectrum the nominal default valueis 128 the actual FFT size used will be max samp seg nfft wi nd specifies the frequency domain smoothing window If wi nd is a scalar the Rao Gabr window 2 2 2 wom n B man em mme T N of length wi nd will be used hereis half the FFT length nf f t and is the set of points m n satisfying 2 wind 2 2 m n mn lt 5 nfft 2 e A unity value for wi nd results in no windowing e f wi nd lt 0 the default value of 5 will be used e f wi nd is a vector it is assumed to specify a 1 D window from which a 2 D window is computed W m n w m w n w
57. recursion 2 85 tricks 1 131 trispect 2 86 trispectrum 1 8 theoretical 2 86 Wigner 2 94 2 96 V variance 1 65 Volterra non Gaussian inputs 1 82 Volterra models arbitrary inputs 2 65 computing output 2 64 Gaussian inputs 2 67 Volterra system 1 80 W wig2 2 87 wi g2c 2 88 wig3 2 90 wig3c 2 92 wig4 2 94 1 5 Index 1 6 wig4c 2 96 wigdat mat 1 137 Wigner bispectrum 1 94 2 90 smoothed 2 92 Wigner cross spectrum 1 90 Wigner spectrum 1 90 2 87 smoothed 2 88 Wigner trispectrum 1 98 2 94 sliced 1 98 smoothed 2 96 Wigner Ville distribution 1 89 window function 1 16 1 90 Wold s decomposition 1 5
58. spectrum and Pw f PX f and P f are the averaged estimates of the power spectra of wx andy bicoher bispecdx 1 Subba Rao T and M Gabr An Introduction to Bispectral Analysis and Bilinear TimeSeries Models pp 42 43 New York Springer Verlag 1984 2 Nikias C L and A Petropulu Higher Order Spectra Analysis A Nonlinear Signal Processing F ramework N ew J ersey Prentice Hall 1993 bispecd Purpose Syntax Description Bispectrum estimation using the direct F FT based method bspec waxis bispecd y bspec waxis bispecd y nfft wind samp seg overlap The bispectrum of the process y is estimated via the direct FFT based method y is the data vector or matrix nfft specifies the FFT length to be used for computing the bispectrum the nominal default value is 128 if nf f t is smaller than samp seg the power of 2 just larger than samp seg will beused wi nd specifies the frequency domain smoothing window If wi nd is a scalar the Rao Gabr window of length wi nd will be used This window is defined by 2 2 2 wom n eee mme T N whereN is half the FFT length nf ft and G is the set of points m n satisfying wind 2 2 m n mnes 2 nfft 2 e A unity value for wi nd results in no windowing e f wi nd lt 0 the default value of 5 will be used e f wi nd is a vector it is assumed to specify a 1 D window from which a 2 D window is computed W m n 2 w m w n w m n
59. symmetric distributed their skewnesses are zero hence only the second component contributes to the bispectrum finally since u4 n is Gaussian all of its cumulants of order greater than two are zero hence it does not contribute to the fourth order cumulant Even though 3 2 0 in the above example it is easy to construct examples of H gt so that the bispectrum is zero 63 Hence a zero bispectrum is not inconsistent with nonzero skewness of the driving i i d process In the context of harmonic processes of the form p s n a cos oyn oy k 1 the phrases frequency coupling and phase coupling are often used in the same sense the former refers to relationships of the form 3 705 and the latter to 03 0 gt 04 Generally but not always the two relationships go hand in hand Further in the case of harmonic processes the power spectrum and parametric methods based on the autocorrelation are usually sufficient unless the data contain narrow band Gaussian components Higher order cumulants and spectra are useful to isolate specific types of coupling quadratic cubic etc Several of the M files in the Higher Order Spectral Analysis Toolbox allow the user to segment the data that is cumulants are estimated for each segment and are then averaged across the set of segments such segmentation usually speeds up calculations at a slight loss in statistical efficiency note that segmented estimates can never be bet
60. than or equal to 4 max delay 1 overlap Specifies the percentage overlap between segments The default value is 50 Theallowed range is 0 99 Ifx y are matrices samp_seg isset tothe row dimension of x andoverl ap isset to zero nfft specifies the FFT length to use the default value is the power of 2 just greater than samp seg delay is the estimated delay of signal y with respect to signal x rxy is the estimate of the windowed autocorrelation described below it is a vector of length nf f t corresponding to indices of nf ft 2 tonfft 2 1 Let S denote the cross spectrum between the two signals x and y and let Sy and Sy denote the auto spectra of x and y The squared coherence function is defined by 2 xy 9 ee S 0 S 0 tder See Also References The optimal ML window is then 1 Cyy W o and the windowed cross correl ation Rxy m istheinverse F ourier transform of W o S o Estimates of Sy Sy and Cy are obtained via the Signal Processing Toolbox routinespectrum R1 the data is segmented and spectrum estimates from individual segments are averaged the segment length is taken to be the smallest power of 2 larger than twicethe maximum delay Since good esti mates demand a large number of segments it is critical that the lengths of the time series x and y be much larger than max del ay The estimated delay will range from nf ft 2 tonfft 2 1 The raw delay d is estimate
61. the noisy signal z w whose SNR is 10 dB Colored Gaussian noise generated by passing a white Gaussian sequence through the AR filter 1 0 0 49 was added to y to obtain zc also at a SNR of 10 dB tdel mat A synthetic for time delay estimation The signal at the first sensor is i i d zero mean single sided exponentially distributed with unity variance and skewness of two Thesignal at the second sensor is a delayed version of the signal at the first sensor delay 16 samples The first signal was corrupted by white Gaussian noise to obtain a SNR of 0 dB The signal at the second sensor was corrupted by colored Gaussian noise obtained by passing the noise at the first sensor through the MA filter 1 2 3 4 5 6 5 4 3 2 1 Note that the two noise signals have a strong spatial correlation at a delay of five samples The signal length is samples Vectorss 1 ands2 contain thesimulated signals at the two sensors tprony mat This mat file contains the 256 sample complex transient x and is used to test the function hpr ony here x n exp jz 2 n 0 1 j0 847 2expfm 4 n 0 2 j0 647 wi gdat mat This filecontains four signals s1 s2 s3 ands4 described below which are used to test the functions wi g2 wi g2c wi g3 wi g3c wi g4 andwi g4c Signal s1 is areal harmonic given by sl n cos 2nfn f 0 1 1 137 1 Tutorial Signalss2 ands3 are generated via T N No s n exp 2n A cos 2zfnT 20 with T
62. the number of samples per segment or record The default value of samp_seg is the length of the time series overlap specifies the percentage overlap between segments The default value is 0 The allowed range is 0 99 Ify isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset tozero andsamp seg is set to the row dimension f I ag should be either biased unbiased for biased unbiased sample estimates of cumulants By default biased estimates are computed If thefirst letter is not b unbiased estimates are computed nf ft specifies the FFT length to be used for computing the bispectrum the default valueis 128 if nt ft is smaller than 2 n ag 1 the power of 2 just larger than 2 nl ag 1 will be used wi nd specifies the lag domain smoothing window If wi nd is 0 the Parzen window will be applied Otherwise the usual unit hexagonal window will be applied The Parzen window is the default 2 31 bispeci Algorithm See Also References 2 32 The 1 D Parzen window is defined by 2 3 1 6 e m L 2 dp m a 1 Im L 2 lt m lt L 0 Im gt L whereL nl ag Theactual window applied tothe estimated cumulants is given by W m n 7 dg m dp n dp m n The unit hexagonal window is given by W m n d m d n d m n where d m 1 m nl ag bspec isthe estimated bispectrum it is an nf f t by nfft array with origin at the center and axes po
63. the plots in the two figures are quite similar A description of the signals x and y in thefilesn 1 mat andnl 1 mat may be found in the section Data Files The true linear and quadratic impulse responses are shown in Figure 1 51 Quadratic Phase Coupling Phase coupling occurs due to nonlinear interactions between harmonic components Three harmonics with frequencies fy and phases K 1 2 3 are said to be quadratically phase coupled if f 3 f f2 and 63 o 09 Quadratic phase coupling coupling at sum and difference frequencies occurs when a signal is passed through a square law device for example and may be detected from the bispectrum 53 Consider the signal x n which is a mixture of harmonics with independent phases and quadratically phase coupled sinusoids that is N onlinear Processes N Na N x n Y ajcos 2nfjn 6 E E oj cos2xfin 6 i l i lk 1 1 163 where fh 3 fib 9 o Orb the f s areall distinct and 74 0 4 and i are all i i d and uniformly distributed over z 7 Then the third order cumulant of x n is given by 66 C3x 4 T1 N 1x24 0 8 22 Q 30 50 3 x i 1 cos 2rf 4 2xfit cos 2nf 47 2nf T gt cos 2nf iht 2nf T gt Cos 2x fib t4 2f T cos 2nf t4 2nfj4t gt cos 2xfist T 2nfi515 1 164 The diagonal slice C3 t t is given by N 3 a 3 1 Cou 5 ci Y cos 2nf14 i 1 k 1 1j 1 165 1 85 1 Tutorial
64. the time series x and y be much larger than the expected maximum delay An initial estimate d of the delay D is given by the location of the peak of R m A three point interpolation may be used to improve the delay estimate 39 20 1 R d R d 1 MEE R d 1 2 R d R d 1 Theoptimal maximum likelihood window estimator is implemented in routine tder Time Delay Estimation Examples load tdel delayztder s1 s2 30 64 0 64 Estimated delay 15 9987 You should see the display on Figure 1 26 and the Estimated delay Estimated 15 997 The two signals are corrupted by spatially correlated noise the noise correlation shows a strong peak at a delay of 5 samples the signal delay is 16 samples TDER Windowed Rxy delay 16 T T T 0 8r 4 0 6r 4 0 4r 0 2r 0 4 L L L L L L L 80 60 40 20 0 20 40 60 80 Figure 1 26 Time Delay Estimated Using Cross Correlation tder A Cross Cumulant Based Method The cross correlation based method assumes that the sensor noises are uncorrelated If the noise processes are correlated it may not be possible to detect the peak of Ras x since it may be masked by Rw Wy t see 1 189 If the signals are non Gaussian and the noise processes are Gaussian we can usethird order cumulants even if the noise processes are correlated 1 103 1 Tutorial 1 104 Let P be the maximum expected delay and assume that the delay D is an integer Then 3
65. tothe AR model order determination problem In the presence of additive noise process x n in 1 124 satisfies the special ARM A p p model p p Y a m y n m Y a m w n m em ea 1 132 where w n is additive noise not the driving noise of an ARMA model Now we can use any of the standard ARMA parameter estimation techniques to esti mate the a k s Multiplying 1 132 by y n k and taking expectations we obtain p 2 Y a m Ry k m oya k We 1 133 where we have assumed that the additive noise is white hence with k gt p we obtain a set of equations from which wecan estimatethe a k s these equations are sometimes called the extended normal equations The analysis extends to the real harmonic case where we obtain an AR 2p or ARMA 2p 2p model with roots at z e K 1 p Sincetheroots are all on z 21 it turns out that a k a 2p k k 20 2p Harmonic Processes and DO A Pisarenko s Method TheM xM autocorrelation matrix of the process y n in 1 124 may be written as Ry Ryxt Sal SDS owl qus where S is a Vandermonde matrix with m n element exp j2n m 1 fn m 1 p and D is the diagonal matrix of the sinusoidal powers D diag Pj Pp Since the frequencies are assumed to be distinct matrix S has full rank p M 2 p It follows immediately that the p largest eigenvalues of the matrix Ry are given by P m k 1 p the remaining M p eigenvalues a
66. vector or a matrix if itis a matrix columns are assumed to correspond to different realizations and each column is processed separately h is the impulse response of the linear part it should be a vector q istheimpulse response of the quadratic part it should be a matrix Note that functions nl ti ck and nl pow assume that q is symmetric y is the output of the nonlinear system and will have the same dimensions as x itis computed via Lt M 1N 1 y n Y h k x n k Y ak bxi xn 0 k 0 k 01 0 where L is the length of n and q is M by N The data files ni 1 mat and nl 2 mat were generated via this function nltick nlpow nlpow Purpose Syntax Description Algorithm See Also Nonparametric second order Volterra System Identification method for arbitrary input signals h q nl pow x y nfft Implements the nonparametric method of Powers amp al 1 for the identification of a second order Volterra system given both the input and output signals x istheinput tothe nonlinear system y is the output process it must have the same dimensions asx Ify isa matrix columns are assumed to correspond to independent realizations nfft isthe FFT length to be used the default value is the power of 2 greater than the length of the time series row dimension for matrices h is the estimated impulse response of the linear part q is the estimated impulse response of the quadratic part The frequency
67. via ymat doagen 1 ymat is an nsamp by msens array wherensamp is the number of samples per sensor record and ms ens isthenumber of sensors Thekth column contains the signal observed at the kth sensor gistat Purpose Syntax Description Computes decision statistics for Hinich s Gaussianity and linearity tests sg sl glstat y sg sl glstat y cparm nfft gl stat estimates the decision statistics for Hinich s Gaussianity and linearity tests The bispectrum is estimated using the direct method and a frequency domain 2 D boxcar smoother is applied The power spectrum is estimated via the direct method and a boxcar smoother is applied The bicoherence is then estimated The Gaussianity test actually zero skewness test basically involves deciding whether or not the estimated bicoherence is zero The linearity test involves deciding whether or not the estimated bicoherence is constant y is the time series should be a vector cparmistheresolution parameter it should lie between 0 5 and 1 0 if asingle record is used the default value is 0 51 Increasing c par m decreases the variance of the smoothed bispectral and spectral estimates but at the expense of poorer resolution nfft isthe FFT length to be used the default length is 128 If the length of y is greater than nf ft y is segmented into records of length nf ft The boxcar window length M isthe value obtained by rounding off nf t t
68. where fo 0 0625 Hz In the case of self coupling 73 the bispectrum can be used to verify the presence of higher order couplings as well we have fo fy 2fo and fo 2fo 3fo fo t3fo 4fo etc which lead to the indicated peaks In this example we used gl st at to confirm that the data are non Gaussian and nonlinear Weusedhar mest andqpct or to estimate the frequencies of the fundamental and the harmonics and to confirm frequency coupling We used roots todetermine the frequencies from the parameter vector estimated via the AR method we can usepi ck peak tolocate the peaks in the various spectra estimated by har mest notethat the accuracy of the peak locations is limited by the FFT length this should not be confused with resolving power Since the signal is rich in harmonics we found it more useful to look at a slice the diagonal slice of the bispectrum Our final conclusions are that the speech segment in question is non Gaussian nonlinear and contains harmonics of the form fok kfo with fy 500 Hz C ase Studies Power Spectrum Singular values o 0 2 0 4 0 6 frequency Estimated Parametric Bispectrum 9 o a T L 1 o 0 05 0 1 0 1 0 2 0 25 0 3 0 35 0 4 fi Figure 1 50 Parametric bispectrum via qpctor 1 0 45 1 129 1 Tutorial 1 130 Table 1 4 Commands Used to Analyze Laughter Data load laughter mat laughter data y y 1 1400 sp y mean y std y
69. which sample cumulants are esti mated samp seg specifies thenumber of samples per segment the default valueisthe length of thetime series overlap specifies the percentage overlap between segments the allowed range is 0 99 and the default valueis O Iff ag is biased then biased sample estimates are computed default if the first letter is not b unbiased estimates are computed If y isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset tozero andsamp seg is set tothe row dimension Cumulants are estimated from each realization and then averaged The estimated AR parameters arereturned in thep 1 element column vector avec 2 17 arrcest Algorithm In the noiseless case the AR parameters are obtained as the least squares solution to the normal equations given by p Y a kK R m k 0 m gt q k 0 p a k C3 m k p 0 m gt q k 0 p Y a k C m k p 0 0 m gt q k 0 TheAR identifiability conditions 1 2 require that p qo P do for any qo and m q i 1 q i p i 2 0 this leads tothe default value of max ag If the additive noise is white Gaussian or non Gaussian the autocorrelation based equations hold only for m gt max p q If the additive noise is non Gaussian and white the cumulant based equations hold only for m gt max p q in these cases choose mi nl ag gt max p q See Also armarts armaqs cumest References 1 Swami
70. with entries C i j 2 C4y i j 0 0 or C i j Cay i j Also let C 2 VSV denote the eigen decomposition where S is the diagonal matrix of eigenvalues A k and V isthe matrix of eigenvectors Vy K 1 maxl ag Let e o Lexp jo exp j max ag 1 o denote the FFT vector and let p denote the chosen order the parameter p order Then the power spectral estimates are obtained as follows E rmt k p 1 where 1 MUSIC w k 1 A k Eigenvector 5 kK m Pisarenko where 6 k is the Kronecker delta function The AR power spectrum is obtained as follows First a rank approximation of the matrix C is obtained as C VSV where is obtained from S by setting A k 0 k p 1 M TheAR parameter vector is then obtained as the solutionto Ca 0 the method in 5 is used and the solution is forced to have unity modulus 2 51 harmest References 2 52 The ML Capon solution is given by ET Pave KK kz1i 2 1 Let V denote the matrix of eigenvectors corresponding to the p largest eigenvalues of R Partition matrix V as v i V r The AR parameter vector for the minimum norm solution is given by 1 a INE V v1 l v v The power spectrum is given by 2 S 1 p ado exp jko k 0 1 Pan R and C L Nikias Harmonic decomposition methods incumulant domains Proc IEEE ICASSP 88 pp 2356 59 1988 2 Swami
71. 15 model based methods 1 15 non parametric methods 1 15 parametric estimators ARMA models 1 26 power spectrum estimator Burg estimator 1 72 Capon s ML estimator 1 72 MVD estimator 1 72 Prony s method 2 56 QPC 1 88 detection 2 71 synthetics 2 70 qpc mat 1 137 qpcgen 2 70 qpctor 2 71 q slice method 1 33 quadratic phase coupling 1 84 quick help 2 55 R random sequence generator 2 76 recursive instrumental variable RIV algorithm 1 56 recursive least squares RLS algorithm 1 56 reflection coefficients 1 49 1 60 residual time series 1 32 Index resolution 1 65 RIV double lattice form 1 58 transversal form 1 56 RIV algorithm 1 56 ri v mat 1 137 rivdl 2 72 rivtr 2 74 RLS algorithm 1 47 rpiid 2 76 S self driving AR model 1 66 signal subspace 1 68 Skewness 1 7 speech signals 1 122 sunspot data 1 108 synthetic generator harmonics in noise 2 53 synthetics 1 65 system identification non parametric 2 19 2 21 2 63 T TDE cross bispectral method 2 80 cross correlation method 1 101 cross cumulant method 1 101 2 77 ML window cross correlation method 2 82 synthetics 2 81 using hologram 1 105 tde 2 77 tdel mat 1 137 tdeb 2 79 tdegen 2 81 tder 2 82 TFD s Cohen dass 1 89 time delay estimation problem 1 101 time frequency distribution 1 89 tls 2 84 total least squares 2 84 tprony mat 1 137 transient signals 1 89 transients modeling 2 56 trench 2 85 Trench recursion 1 49 trench
72. 2 70 Generates synthetics for the quadratically phase coupled harmonics in colored Gaussian noise problem zmat qpcgen zmat qpcgen default qpcgen generates independent realizations of the signal 3 y n L Aki COS 2nA in OL x COS 2n n OK g n b p where A3 242 A1 and O43 702 6k kL 6 amp 2 Ok are mutually independent and uniformly distributed over 0 27 If qpcgen isinvoked without any input arguments then you are prompted for the length of the realizations the number of realizations the number of phase coupled triplets p their frequencies 4 and 4 2 and amplitudes o as well as the corresponding parameters in the uncoupled case the terms with an overbar and the variance of the colored Gaussian noise g n The colored noise is generated by passing a white Gaussian noise sequence through an user specified ARMA filter you are prompted for these ARMA parameters The phases 6 4 6 2 and i i 1 2 3 are chosen randomly for each realization Note that noise free realizations can be obtained by specifying a value of zero for the noise variance If the function is invoked as qpcgen def aul t wherethe variabledefaul t may take on any value s then the default settings are used see Examples below Each column of z mat corresponds to a different realization The matrixz mat in the file qpc mat can be regenerated via zmat qpcgen 1 qpctor Purpose Syntax Description Algorithm
73. 356 59 New York 1988 3 Roy R and T Kailath ESPRIT Estimation of signal parameters via rotational invariancetechniques IEEE Trans ASSP Vol 37 pp 984 95 J uly 1989 4 Swami A and J M Mendel Cumulant based approach to the harmonic retrieval and related problems IEEE Trans ASSP Vol 39 pp 1099 1109 May 1991 2 45 doagen Purpose Syntax Description 2 46 Generates synthetics for the direction of arrival DOA problem using a uniform linear array ymat doagen ymat doagen default doagen generates synthetics for the DOA problem The sensor array is assumed to belinear and equispaced uniform z mat returns the generated data each column corresponds to a different realization If doagen is invoked without any input arguments then you are prompted for all parameters the source bearings number of sensors ms ens sensor spacing number of samples per sensor record ns a mp pdf of source signal variance of additive noise pdf of additive noise and ARMA parameters for the noise spectrum The sensor array is assumed to be linear and uniformly spaced If theuser specified noise variance is greater than zero the sensor signals are standardized to unity variance before the additive noise is added to them If the function is invoked as doagen def aul t wherethe variabledefaul t may take on any value s then the default settings are used Matrix ymat in filedoal mat was generated
74. 41 we use the autocorrelation C5y n n 1 q hence the algorithm can handle additive white noise but not colored noise Examples load mal bvec maest y 3 3 128 The esti mated parameters should be 1 0 9714 0 3814 0 7759 Thetrue parameters are 1 0 9 0 385 0 771 and the signal is contaminated with white Gaussian noise bvec maest y 3 4 256 The estimated parameters should be 1 0 9608 0 4482 0 7343 Polyspectra and Linear Processes AR Models The cumulants of the noisy process satisfy the normal equations p Y a K Cg m k 0 m q ae 1 42 p a k C3 m k p 0 m gt q k 0 1 43 p Y a k Cyy m k p t O0 m gt q k 0 1 44 which can be verified by expressing the cumulants in terms of the IR h n using 1 16 through 1 18 and by noting that yt pack yh n k b n Therelationship between linear prediction LP and AR models is discussed in detail in the section titled Polyspectra and Linear Processes Here we will just point out that the least squares solution to the LP problem is given by 1 42 with m 1 p these equations are called the normal equations because the resulting prediction error sequence is orthogonal to the data However when q gt 0 or when additive noise is present in the data the normal equations yield inconsistent estimates of the AR parameters E quations 1 42 through 1 44 yield consistent estimates and can also be derived by
75. 7 These frequencies are normalized that is to find the correct frequencies we should multiply them by the sampling frequency which was 8000 Hz sval cum 2 pxx 20 0 0 100 0 10 20 30 0 5 0 0 5 eig a im eo music o pisar o ec ml e o o m i 1 o 0 5 0 0 5 0 5 0 0 5 0 0 E o 50 T 20 E 100 40 0 5 0 0 5 0 5 0 0 5 Figure 1 45 Laughter Data Power Spectra C ase Studies We can also use fourth order cumulants to esti mate the frequencies Notice that the singular values of the cumulant matrix shown in the 1 1 panel of Figure 1 46 look quite different from the singular values of the covariance matrix shown in the 1 1 panel of Figure 1 45 in particular the singular values of the cumulant matrix appear to be essentially zero after p 8 a possible explanation is that some of the harmonics are well modeled as narrow band Gaussian and some as narrow band non Gaussian Various spectra corresponding to order p 8 are shown in Figure 1 46 the root locations of the AR polynomials are shown in the bottom panels of Figure 1 47 The dominant roots of the AR polynomial in the AR method are located at 0 1297 and 0 1964 Hz a less dominant pole is at 0 2678 Hz In the case of the Min Norm method the dominant roots are located at 0 0616 0 1307 0 1960 and 0 2622 suggesting that the signal has harmonics at approximately 2f o 3fo and
76. 7 f4 f2 is nonzero everywhere Why Do We Need Higher Order Statistics Motivation to use cumulants and polyspectra of order k gt 2 is given by the following my 2 mi my3 e f z n x n y n and x n and y n are mutually independent processes then Ckz M Cj Imi Cyy my e If x n is Gaussian then C m 0 k 2 Hence if z n 2x n w n where w n is Gaussian and independent of x n then for k 72 Ciz my Ci my Thus we can recover the higher order cumulants of a non Gaussian signal even in the presence of colored Gaussian noise Let x n bea linear process that is x n y hou k where u n is i i d Then it follows that C k Youd h h n k n 1 16 C3 k I Tau h h n k h n l n 1 17 C k l m h n h n k h n D h n m 4x Youd 1 18 Polyspectra and Linear Processes Sj YoylH COT 1 19 53x f f2 Ya JH FH f2 H Fy f2 1 20 S Jf H H H H fi f5 ax fr fo f3 Yau FDH GEB F3 H F1 fot f3 1 21 where Yku Cku 0 Note that the power spectrum does not carry any information about the phase of H f n contrast if u n is non Gaussian this phase information can be recovered from the higher order polyspectra Thus the standard minimum phase assumption which is necessary when the process is Gaussian or only second order statistics are used may be dropped Any process can always be considered to be a linear process
77. 9 40 P yin aci x n i wn Lee 1 190 where a n 0 n z D and a D 1 Consider the third order cumulants Cy p ELy n x n 1 x n p age Cyyx t P E X N X N 7 x n p aes Substituting 1 190 into 1 191 we obtain P Cx P Y aC g t i p i aed 1 193 Using this equation for various values of p and t we get a system of linear equations in the a i s namely Coa Cy The estimated delay is the index n which maximizes a n A low rank approximation of the cumulant matrix Cxxx may be used to improve the robustness to noise This algorithm is implemented in routinet de Time Delay Estimation Examples load tdel delay avec tde sl1 s2 30 128 Estimated delay 16 You should see the display in Figure 1 27 and the Estimated delay Estimated delay 16 The two signals are corrupted by spatially correlated noise the cross correlation between the noises at the two sensors has a strong peak at a delay of 5 samples the signal delay is 16 samples TDE parameter vector delay 16 1 T T T 0 6 J 30 20 10 0 10 20 30 Delay in samples Figure 1 27 Time Delay Estimated Using Cross Cumulants tde A Hologram Based Method Thealgorithm described in the previous section can also beimplemented in the frequency domain Define auto and cross bispectra Bio f 1 f2 EUG DXGX G f2 m Sox r f2 Bxyz 3 fo EUG XQ 23 Sfi pye 1 105 1 Tuto
78. A and J M Mendel ARMA parameter estimation using only output cumulants IEEE Trans ASSP Vol 38 pp 1257 65 J uly 1990 2 Swami A and J M Mendel I dentifiability of the AR parameters of an ARMA process using cumulants IEEE Trans on Automatic Control Vol 37 pp 268 73 Feb 1992 2 18 biceps Purpose Syntax Description Estimates impulse response via lag domain bicepstral method hest ceps a b minh maxh biceps y p q hest ceps a b minh maxh biceps y p q samp seg overlap flag Ih bi ceps estimates the impulse response of the linear process y using the bicepstrum lag domain method Variables p and q denote the number of causal and anticausal cepstral parameters to be estimated The total length of the complex cepstrum is p q t1 Variablessamp seg overlap and lag control the manner in which sample cumulants are esti mated samp seg Specifies the number of samples per segment Its default valueis the row dimension of y if y is a row vector the column dimension is used as the default value overlap specifies the percentage overlap between segments the allowed range is 0 99 and the default valueis O Ify isa matrix the columns are assumed to correspond to independent realizations or records in this caseoverlap is set to zero andsamp seg is set tothe row dimension Ifflag is b then biased sample estimates of cumulants are estimated default if the first letter i
79. C3 k k Cy k m Cy I K m Cy k m 1 C44 k I k m k Hence the fundamental region of support is not the entire k D plane For example for k 2 C gt k k 0 specifies C gt k everywhere It is easily shown that the nonredundant region for C3 k l is the wedge k 1 0 lt I lt k lt and for C4x k l m it is the cone k l m O0 lt msl lt k lt oo The kth order polyspectrum is defined as the FTs of the corresponding cumulant sequence SX Y Caeci k zi 1 8 3 2xf k j2nfl Sin yy eager ee Kee Vane 1 9 Z j2 k Sax f p gt fo f3 gt Y Ca k me j2n fik fal f3m Relies 1 10 which are respectively the power spectrum the bispectrum and the trispectrum Note that the bispectrum is a function of two frequencies whereas the trispectrum is a function of three frequencies In contrast with the power spectrum which is real valued and nonnegative bispectra and trispectra are complex valued Polyspectra and Linear Processes For a real valued process symmetry properties of cumulants carry over to symmetry properties of pol yspectra The power spectrum is symmetric SX f S2 f The symmetry properties of the bispectrum are given by 56 S3x fy fo S3x fo f 534Gp 7 f3 f2 MEENE E 1 11 Hence a nonredundant region of support for the bispectrum is the triangle with vertices 0 0 1 3 1 3 and 1 2 0 recall that we have assumed a normalized sampli
80. Classification Example 000 e eee eee 1 120 Laughter Data irete ai scd ER AERIAL TIERRA 1 122 Pitfalls and Tricks of the Trade 0c cee eee eee 1 131 Data Files 0 0 0 ccc eens 1 134 References uou aene Mid ends s 1 139 Reference Function Tables esee 2 2 Higher Order Spectrum Estimation Conventional Methods 2 2 Higher Order Spectrum Estimation Parametric Methods 2 3 Quadratic Phase Coupling QPC 0c eee eee eee 2 3 Second Order Volterra Systems 0 020 cece eee 2 4 Harmonic Retrieval 0 00 eee eens 2 4 Time Delay Estimation TDE 20 0c eee eae 2 4 Array Processing Direction of Arrival DOA 2 4 Adaptive Linear Prediction llle 2 5 Impulse Response IR Magnitude and Phase Retrieval 1 0 cece cee eee eens 2 5 Time Frequency Estimates 000 c eee eee 2 5 Utilities ova Ses weed aba wend Ses Ove esa Move ora imi 2 6 Demo tad eek EEA ACCEDI ne a oe ee 2 6 Miscellaneous oil ge Ps We DE one 2 7 Promptingz sts fi Sate NA Selah A ANS RE ER ea Re SY 2 7 G ided COUP 2 2 nu mA mm meh Mera cen Am Re mm ha 2 7 Addenda ee m usw d RM RE Ud Wake ERAEN 2 7 About the Authors About the Authors About the Authors vi Ananthram Sw ami Ananthram Swami received his B Tech M S and Ph D degrees in electrical engineering from the Indian Institute of Technology at Bombay Rice Univer
81. D k 01 0 Y Y ak DE x n x n xin p k xin p 1 388 1 161 which follows since x t is symmetrically distributed Assuming Gaussianity we can rewrite the last equation in the frequency domain as Syxx fa fa 20 25 G 059 44Q02 02 8 f1 f Ety T 1 162 which is obtained under the assumption that q k I q I k The cross bispectrum Sy f1 f2 can be estimated via routinebi specdx andthe spectra can be estimated via routines pect rum we can then estimate the quadratic transfer function Q f4 f2 from 1 162 This algorithm is implemented in routinen ti ck 70 Notethat the input process is assumed to be Gaussian 1 81 1 Tutorial 1 82 Examples load nll h q nltick x y 128 1 You should see the display on Figure 1 16 A description of the signals x and y innl 1 mat is given in the section Data Files The true linear and quadratic impulse responses are shown in Figure 1 51 linear TF quadratic TF o 0 2 0 4 0 6 0 4 0 2 0 0 2 0 4 f linear part IR 0 4 0 3 0 2 0 1 50 100 150 Figure 1 16 Transfer Functions Estimated by nltick Solution Using FIs In 50 another algorithm to estimate the parameters of the model in 1 157 was developed without assuming that x n is Gaussian indeed x n may be deterministic Given the input and output sequences x n and y n we can compute their FTs X f and Y f on a suitable grid fm
82. D 1 153 We also want to minimize the contribution due to sinusoids at frequencies other than fy this can be achieved by minimizing the output variance in 1 152 subject to the constraint in 1 153 Theresulting estimator is given by Pau te coR decor 1 154 where e 1 e 727 ej2xP 1H This estimator which is variously called Capon s maxi mum likelihood estimator or the Minimum Variance Distortionless estimator 8 27 is not a true power spectral density estimator For example if x n consists of a harmonic at frequency fy and power P that is observed in white noise with variance o then P n o P o2 p where p is the number of autocorrelation lags used note that the noise power contribution has been reduced by a factor of p The resolving power of this estimator is reported to be 27 1 pVSNR where SNR is the signal to noise ratio ratio of signal variance to noise variance and p is the length of the autocorrelation sequence This estimator is implemented in routines har mest and doa Af gt 1 73 1 Tutorial 1 74 The maximum entropy spectral estimator of Burg is based on the availability of exactly known autocorrelation lags Ry m5 2p in which casea Gaussian pdf and an AR p model are obtained for process x n Routine trench may beusedtoestimatethe AR model parameters Theresolving power is reported to be 1 For SNR 1 the maxi mum entropy estimator has better resolving power
83. E Cxyz K D Cyyz K 1 Usually we set N3 to N and obtain estimates that are asymptotically unbiased Autocumulants are obtained when w x y z These estimates are known to be consistent provided the process x n satisfies some weak mixing conditions 5 For example for large N the variance of the sample estimate of the third order cross cumulant is given by AN var Cxyz k b C N where cis a finite constant that depends upon the auto and cross moments cumulants of orders 1 through 6 of the processes x n y n and z n These definitions assume that the processes are zero mean in practice the sample mean is removed first Routinescum2x cum3x and cum4x may be used to estimate cross cumulants of orders 2 3 and 4 cumest may be used to estimate the autocumulants 1 13 1 Tutorial Examples We will simulate a non Gaussian ARMA process and then estimate its cumulants rand seed 0 randn seed 0 u rpiid 1024 exp n 25 ysfilter 1 2 1 1 5 0 8 u for k n n cmat k n 1 cumest y 3 n 128 0 biased k end subplot 121 mesh n n n n cmat subplot 122 contour n n n n cmat 8 Time series y is segmented into records of 128 samples each with no overlap biased estimates of thethird order cumulants are obtained from each segment and then averaged the i j element of c mat will contain the estimate of Czyli n Lj n 1 for ij 21 2 n 1 You can usethefunction cumt r ue to com
84. Estimated Angular Spectra doa 1 78 Harmonic Processes and DO A Summary The estimation of the frequencies of harmonics observed in noise or that of determining the bearings of sources in the DOA problem can be treated from either a parametric or a nonparametric view point In the latter case we obtain nonparametric estimates of spectra and angular spectra n theformer case we have several high resolution algorithms Routinehar mest can beused for the harmonic retrieval problem the user has the choice of using second or fourth order cumulants This routine estimates power spectra using the MUSIC Eigenvector Pisarenko ML Capon AR and Min Norm methods based on either the diagonal slice of fourth order cumulants or on the covariance it also estimates the conventional periodogram The number of harmonics can usually be determined by examining the singular value plot generated by har mest Routinedoa can be used to solve the direction of arrival problem DOA for a uniformly spaced linear array It estimates angular spectra using the Beamformer ML Capon AR MUSIC Pisarenko s method eigenvector Min Norm and ESPRIT algorithms based on the spatial cross correlation or cross cumulant Routinepi ck peak may be used to pick peaks in the estimated angular spectra 1 79 1 Tutorial Nonlinear Processes 1 80 The simplest nonlinear system is the second order Volterra system whose input output relationship is de
85. F 0 and B 0 Its default valueis 1 where the sign is the sign of the cross correlation between the time series y and the corresponding instrumental variable z corresponding to order mor der and unity value for ambda z is computed via z ivcal y morder 1 nsmut h is the window length for estimating the steady state AR parameters Its default value is min length y 4 50 arvec istheAR parameter vector corresponding to the steady state final weight vector and has length ar order 1 Thelast ns mut h samples of the time varying weights wt are averaged and then converted to the AR parameters f pe is the final prediction error wt isansamp xarorder array of theestimated weights of the adaptivefilter as a function of time herensamp is thelength of thetime series y Details of the algorithm are given in the Tutorial i vcal ri vdl rivtr References 1 Swami A and M Mendel Lattice Algorithms for Recursive Instrumental Variable Methods USC SIPI Report 117 University of Southern California Los Angeles Dec 1987 2 Swami A System Identification Using Cumulants Ph D Dissertation University of Southern California pp 107 8 1988 2 75 rpiid Purpose Syntax Description 2 76 Generates i i d random sequence u rpiid nsamp u rpiid nsamp in type pspike rpiid generates a sequence of nsamp i i d random variables of the type described byin type in_t
86. Higher Order Spectral Analysis T ool box For Use with MATLAB Computation Visualization EXE Programming User s Guide Vesi on 2 x OD e E How to Contact The MathWorks 508 647 7000 Phone 508 647 7001 Fax The MathWorks Inc Mail 24 Prime Park Way Natick MA 01760 1500 http www mathworks com Web ftp mathworks com Anonymous FTP server comp soft sys matlab Newsgroup supportQmathworks com Technical support suggestQmathworks com Product enhancement suggestions bugsQmathworks com Bug reports docQmathworks com Documentation error reports subscribe mathworks com Subscribing user registration servicegmathworks com Order status license renewals passcodes info mat hworks com Sales pricing and general information Higher Order Spectral Analysis Toolbox User s Guide COPYRIGHT 1984 1998 by The MathWorks Inc All Rights Reserved The software described in this document is furnished under a license agreement The software may be used or copied only under the terms of the license agreement No part of this manual may be photocopied or repro duced in any form without prior written consent from The MathWorks Inc U S GOVERNMENT If Licensee is acquiring the Programs on behalf of any unit or agency of the U S Government the following shall apply a For units of the Department of Defense the Government shall have only the rights specified in the license under which the commercial computer software or com
87. Hinich M J Testing for Gaussianity and linearity of a stationary time series J TimeSeries Analysis Vol 3 pp 169 76 1982 2 Patel J K and C B Read Handbook of the Normal Distribution Sec 7 9 4 M Dekker New York 1982 2 49 harmest Purpose Syntax Description Estimation of frequencies of harmonics in colored Gaussian noise and power spectra Pxx arl ar2 harmest y Pxx arl ar2 harmest y maxlag p order flag nfft norder har mest estimates the frequencies of real harmonics in noise and power spectra using the MUSIC Eigenvector Pisarenko ML Capon AR and minimum norm methods based either on the diagonal slice of fourth order cumulants or on the covariance it also estimates the conventional periodogram y isthedata matrix each of its columns is assumed to correspond to a different realization maxlag specifies the number of cumulant lags to be computed max ag should be greater than twice the maximum number of harmonics expected The default value of max ag is arbitrarily set tonsamp 12 where nsamp is the row dimension of the data matrix p order specifiestheorder must be greater than or equal totwicethe number of harmonics if this parameter is not specified or is not positive you are prompted for the value of p order The order can be inferred from the display of thethe singular values of the cumulant matrix as follows If the singular values o k are more or less constant fo
88. Kay S M Modern Spectral Esti mation Theory and Applications New J ersey Prentice Hall 1988 28 Kay S M and S L Marple Spectrum Analysis A Modern Perspective Proc IEEE Vol 69 pp 1380 1418 1981 29 Krauss T J N Little and L Shure MATLAB Signal Processing T ool box User s Guide The MathWorks Inc 1994 30 Kumaresan R and D W Tufts Estimatingtheangles of arrival of multiple plane waves IEEE Trans AES Vol 19 pp 134 39 1983 31 Kung S Y et al State space and SVD approximation methods for the harmonic retrieval problem J Opt Soc Amer Vol 73 pp 1799 1811 1983 32 Makhoul J Stable and efficient lattice methods for linear prediction IEEE Trans ASSP Vol 25 pp 423 28 Oct 1977 33 Mallat S A theory for multi resolution signal representation the wavelet transform IEEE Trans PAMI Vol 11 pp 674 93 J uly 1989 34 Marple S L J r Digital Spectral Analysis with Applications pp 229 31 New J ersey Prentice Hall 1987 35 Marple J r A new autoregressive spectrum analysis algorithm IEEE Trans ASSP Vol 28 pp 441 50 Aug 1990 36 Matsuoka T and T J Ulrych Phase estimation using the bispectrum Proc IEEE Vol 72 pp 1403 11 1984 37 Mendel J M Tutorial on higher order statistics spectra in signal processing and system theory Theoretical results and some applications Proc IEEE Vol 79 pp 278 305 1991 1 141 1 Tutor
89. L L 1 31 m 4000 3000 2000 1000 0 1000 2000 3000 4000 Figure 1 38 Lynx Data Differenced 1 117 1 Tutorial 1 118 x 10 sval cum 2 pxx 2 0 1 d i 0 50 0 10 20 30 40 0 5 0 0 5 0 0 re D 2 20 50 40 0 5 0 0 5 0 5 0 0 5 0 0 8 MM 20 LM a 50 40 0 5 0 0 5 0 5 0 0 5 0 0 E i o 50 ft 20 amp E 100 40 0 5 0 0 5 0 5 0 0 5 Figure 1 39 Lynx Data Differenced Power Spectra Bispectrum estimated via the indirect method T T T T T T T T 0 37 Q HMM QW 01r CD SIS f2 o T 0 1 SERE e e 0 5 0 4 0 3 0 2 0 1 0 0 1 0 2 0 3 0 4 Figure 1 40 Lynx Data Differenced bispectrum C ase Studies Table 1 3 Listing of eda m function eda sp E DA Exploratory Data Analysis eda sp gt Sp is the time series data amp Author A Swami Jan 1995 disp 1 data and histogram figure 1 subplot 211 plot 1 length sp sp grid subplot 212 hist sp grid disp 2 Summary stats cl mean sp c2 cumest sp 2 c3 cumest sp 3 c2 31 2 c4 cumest sp 4 c2 2 fprintf Mean g n cl fprintf Variance g n c2 fprintf Skewness normalized g n c3 fprintf Kurtosis normalized g n c4 disp 3 power spectra and harmonic models figure 2 p2 a2 b2 harmest sp 30 0 u 256 2 r cpl xpair roots
90. P 37 m sg thestatisticfor the Gaussianity test is a three element vector sg 1 is the estimated statistics sg 2 is the number of degrees of freedom df p sg 3 is the probability of false alarm Pf a More specifically it is the probability that a x random variable with p degrees of freedom could have a value larger than the estimateds in sg 1 If this probability is small say 0 05 then we may reject the hypothesis of zero skewness at a Pf a or significance level of 0 05 In other words if you decide to accept the hypothesis that the data have nonzero skewness then the probability that the data may actually have zero skewness is given bys g 3 If Pfa islarge then the hypothesis of zero skewness cannot be easily rejected 2 47 gistat Algorithm 2 48 sI the statistic for the linearity test is a three element vector see below for more details s 1 is the estimated statisticR sI 2 is the estimated parameter A this parameter is called As in Hinich s paper sI 3 is the theoretical value of R Thelinearity hypothesis should berejected if the esti mated statistic R is much larger or much smaller than the interquartile range of X50 the x distributed random variable with 2 degrees of freedom and noncentrality parameter 1 In practice for a nonlinear process the estimated R value may be expected to be much larger than the theoretical R value 1 The number of samples available to es
91. SINQGFTS eeeeeeee re 1 82 EXAMPLES 22 wavy ares nthe adds bese eee ea tevee Yun 1 83 Quadratic Phase Coupling 00 c cece eee 1 84 Examples 5 s wisi uero TE RE oe ES 1 87 S rntnidEy sn eet um A Aiea we Ge ME Aen e are ARR n dt RR A 1 88 Time Frequency Distributions 1 89 Wigner Spectrum 0 00 nen 1 90 Examples esr EE EE bees Wee ERE 1 93 Examples seek ote eve ea E edo Gnd 1 94 Wigner Bispectrum ssesere eh 1 94 Exampl6s stes et eR Kern ate Aen Ron ge ee 1 96 Examples so teste eee Rade Soles NE ER UE T RE ERA ERE 1 97 Wigner Trispectrum lessen Ih 1 98 Examples 5S2 e eee eed VENUES Se EE au E Qt EE 1 99 ExampleS 23s 15 eee AER ESAE e ta heal enfe Ae 1 100 SUMMARY 23 inis i rp nete Sa en Da fre inf d ev DE ded Ts o s 1 100 Time Delay Estimation ss 1 101 A Cross Correlation Based Method 2 05 1 101 Examples cse ied c eR ec e Rep ee 1 103 A Cross Cumulant Based Method lsuslsusu 1 103 Examples 2222 ri AAT NEA ES a a ER EU A E ue 1 105 A Hologram Based Method sselsslsesesss 1 105 Examples 222 vvesle E Ives Paved I IIEfCYREEUEPPUYS 1 107 SUMMARY uere ee Ee WR e ERU T qx Ra e 1 107 iv Contents Case Studies oues PRIN SA 1 108 S nspotData oix cee did E X ONE DA 1 108 Canadian Lynx Data 0 00 es 1 114 Examples 2 a heh et R R ie pee E RRISMPS 1 114 A
92. a2 disp Estimated cycles a angle r yest a 0 ind find a 0 yest ind 2 pi ones size ind angle r ind yest disp 4 gaussianity and linearity tests glstat sp 51 length sp disp 5 the bispectrum fi gure 3 Bspec w bispeci sp 25 25 lags Wpcolor w w abs Bspec shading interp nice slow 1 119 1 Tutorial 1 120 A Classification Example The data consist of two underwater acoustic signals the sampling rate is fairly high and we have lots of data 131 072 samples The primary interest here is to classify the two signals Means and standard deviations were estimated over 262 nonoverlapped segments each of length 500 samples and are shown in Figure 1 41 the data appear to be highly nonstationary Because of the nonstationarity overall estimates of the power spectrum or the bispectrum or tests of Gaussianity and or linearity will not be meaningful mean A mean B 1 I 1 0 5 0 5 0 0 0 5 ai 1 1 0 100 200 300 0 100 200 300 Std A std B 1 1 0 8 0 8 0 6 0 6 0 4 0 4 0 2 l 0 2 MIU UI 0 100 200 300 0 100 200 300 Figure 1 41 Acoustic Data Means and Standard Deviations We segmented the data into lengths of 1024 samples power spectra esti mated from each segment indicate that the data are essentially low pass we expe
93. alizations mal mt An MA 3 synthetic the MA parameters are 1 0 9 0 385 0 771 the time series length is 1024 samples input distribution is single sided exponential additive white Gaussian noise was added so as to obtain a SNR of 20 dB The vector y contains the MA 3 synthetic The model is mixed phase nl1 mat This file contains data to test functions n pow andnlti ck A white Gaussian process x is passed through a second order Volterra system to obtain the output process y Matrices x andy are 64 x 64 with columns corresponding to realizations The input output relationship for a second order Volterra system is given by 1 135 1 Tutorial 1 136 The filter h was chosen to bethe IR of an FIR 12 filter with a cutoff of 0 2Hz and was generated with h fir1 11 0 4 The filter q was chosen such that co co y n h k x n K Y X q k Dx n k x n I k 0 k 01 0 q k I h1 k h1 I wherehlisthelR of an FIR 12 filter with a cutoff of 0 1Hz and was generated byh1 fir1 11 0 2 Theoutput data were generated using function nl gen y nl gen x h q Theimpulse responses IR and transfer functions TF of the linear part h and the quadratic part q are shown in Figure 1 51 IR trueh IRtrue q 0 4 40 0 3 30 0 2 20 0 1 10 0 Hd 0 5 10 15 10 20 30 40 TF trueh TF trueq TO NTSECUU ESI TERETE SERIEN uui Bonn 04 10 0 2 10 0 10 0 2
94. atial covariance or fourth order cumulant matrix spec theta dvec doa ymat spec theta dvec doa ymat dspace dtheta nsource order delta doa estimates the angular spectra of source bearings for a uniformly spaced linear array using the Eigenvector Music Pisarenko ML Capon AR minimum norm beamformer and ESPRIT methods based either on fourth order cumulants or the spatial covariance matrix The peaks in the angular spectra nominally indicate the direction of arrival DOA various algorithms based on the diagonal slice of the fourth order cumulant or the spatial cross correlation ymat is thedata array each column corresponds to a different sensor the rows correspond to snapshots dspace istheelement spacing in wavelength units the default valueis 0 5 half wavelength spacing dtheta is the angular spacing in degrees at which the spectra are to be computed the default value is 2 degrees nsource is the number of sources the default value is 0 If nsource is not positive you are prompted to choose the number of sources This value can be inferred from the display of the singular values of the covariance or the fourth order cumulant matrix as follows If the singular values o k are more or less constant for k gt p then a useful rule of thumb is to choose p as the number of sources With finite data records you can expect a slow decrease in the smaller singular values Usually there will
95. attice algorithms for recursive instrumental variable methods to appear in Thelnt l J ournal of Adaptive Control and Signal Processing 1996 See also USC SIPI Report 117 University of Southern California Los Angeles Dec 1987 Proc CASSP 88 New York pp 2248 51 Apr 1988 Proc Amer Control Conf Atlanta GA pp 2114 19 J une 1988 68 Swami A and J M Mendel ARMA parameter estimation using only output cumulants IEEE Trans ASSP Vol 38 pp 1257 65 J uly 1990 Swami A and J M Mendel Cumulant based approach to the harmonic retrieval and related problems IEEE Trans ASSP Vol 39 pp 1099 1109 May 1991 References 70 Swami A and J M Mendel I dentifiability of the AR parameters of an ARMA process using cumulants IEEE Trans on Automatic Control Vol 37 pp 268 73 F eb 1992 71 Swami A and P Tichavsky Then th order moment cumulant of multiple sinusoids strong ergodicity and related issues Proc IEEE SP ATHOS Workshop on HOS Girona Spain J une 1995 72 Tekalp and A T Erdem Higher order spectrum factorization in one and two dimensions with applications in signal modeling and nonminimum phase system identification IEEE Trans on Acoustics Speech and Signal Processing Vol 37 10 pp 1537 49 1989 73 Tick L J The estimation of transfer functions of quadraticsystems Technometrics 3 563 67 Nov 1961 74 Tugnait J K Identification of linear stochastic
96. ay also be defined by conjugating one or more of the terms x y and z This is readily accomplished by using the MATLAB function conj x y andz are segmented into possibly overlapping records the mean is removed from each record and the FF T computed The cross bispectrum of the kth record is computed as Byyz k M n X MY Zi n n where X Yk and Z are the FFT s of the kth segments of x y and z The bispectral estimates are averaged across records and an optional frequency domain smoother specified by parameter wi nd is applied bispecd 1 Subba Rao T and M Gabr An Introduction to Bispectral Analysis and Bilinear TimeSeries Models pp 42 43 New York Springer Verlag 1984 2 Nikias C L and A Petropulu Higher Order Spectra Analysis A Nonlinear Signal Processing Framework New J ersey Prentice Hall 1993 bispeci Purpose Syntax Description Bispectrum estimation using the indirect method bspec waxis bispeci y nlag bspec waxis bispeci y nlag samp seg overlap flag nfft wind The bispectrum of the process y is estimated via the indirect method y is the data vector or matrix nl ag specifies the number of cumulant lags to be computed the third order cumulants of y C3y m n will be estimated for nl ag lt m n nl ag This parameter must be specified A useful rule of thumb is to set nl ag tons amp 10 wherens amp is the length of the time series y samp_seg specifies
97. be a significant drop in the singular value from o p to o p 1 order specifies the cumulant order to use it should be either 2 for cross correlation based estimates or 4 for estimates based on the diagonal slice of the fourth order cross cumulant The default value is 4 delta is the displacement between the two subarrays for ESPRIT here we assume that the two subarrays are obtained by partitioning a single array If there are m sensors the column dimension of array y mat then the first array doa Algorithm See Also References will consist of sensors1 2 m delta andthe second array will consist of sensors delta deltatl m Thedefault value of delta is 1 spec isthearray of estimated spectra the columns correspond to esti mates based on the Eigenvector Music Pisarenko ML Capon AR minimum norm and beamformer methods the rows correspond to bearing angles which are returned in the vector t het a All estimated spectra are normalized to a maxi mum value of unity theta is the vector of bearings corresponding to the rows of spec dvec is the vector of bearings estimated by ESPRIT Details of the algorithm are given in the Tutorial har mest 1 J ohnson D H The application of spectrum estimation methods to bearing estimation problems Proc IEEE Vol 70 pp 975 89 1982 2 Pan R and C L Nikias Harmonic decomposition methods in cumulant domains Proc I CASS P 88 pp 2
98. be wrong In other words in a Monte Carlo simulation of 1000 trials one should expect the test results to be wrong 1000 pfa times Polyspectra and Linear Processes Examples load mal q maorder y 0 6 The following table will be displayed q var cqk cqk thres result 0 71 17383e 003 71 31541e 001 5 24957e 001 0 1 6 60864e 002 1 37221e 001 1 59332e 001 1 2 3 76124e 002 4 11813e 001 3 80114e 001 0 3 1 04832e 003 5 06864e 003 1 64547e 001 1 4 3 11151e 003 3 31784e 002 1 20463e 001 1 5 1 31938e 003 6 87545e 002 71 11923e 002 1 6 4 41359e 003 1 43717e 003 1 31092e 001 1 Estimated MA order is 3 The columns in the table correspond to the estimated variance of c3 q 0 the estimated value of c3 q 0 the corresponding threshold and whether or not the absolute value of the estimated c3 q 0 exceeded the threshold The time series y corresponds to an MA 3 process contaminated by AWGN with SNR of 20 dB Linear Processes Impulse Response Estimation The basic model here is x n Y h k u n k em 1 53 y n x n w n 1 54 where additive noise w n is assumed to be symmetric distributed not necessarily Gaussian Process u n is i i d non Gaussian independent of the noise and satisfies 0 C3 0 0 ee 1 37 1 Tutorial 1 38 Wewill discuss two algorithms both based on the notion of the logarithm of the bispectrum The Polycepstral Methods It is assumed that H z has no zeros on t
99. c Am Vol 91 2 pp 975 88 Feb 1992 48 Pisarenko V F The retrieval of harmonics from a covariance function Geophysical J Royal Astron Soc Vol 33 pp 347 66 1973 References 49 Porat B Digital Processing of Random Signals New J ersey Prentice Hall 1994 50 Porat B B Friedlander and M Morf Square root covariance ladder algorithms in Proc ICASSP 81 Atlanta GA pp 877 880 March 1981 51 Powers E J Ch P Ritz C K An S B Kim R W Miksad and S W Nam Applications of digital polyspectral analysis to nonlinear systems modeling and nonlinear wave phenomena Proc Workshop on Higher Order Spectral Analysis pp 73 77 Vail Colorado J une 1989 52 Priestley M B Spectral Analysisand TimeSeries Academic Press London 1981 53 Priestley M B Non linear and non stationary time series analysis Academic Press London 1988 54 Raghuveer M R and C L Nikias Bispectrum Estimation A Parametric Approach IEEE Trans ASSP Vol 33 pp 1213 30 1985 55 Rangoussi M and G B Giannakis FIR modeling using log bispectra Weighted least squares algorithms and performanceanalysis IEEE Trans Cir Sys Vol 38 pp 281 96 1991 56 Rioul O and M Vetterli Wavelets and signal processing IEEE Signal Processing M agazine pp 14 38 Oct 1991 57 Roy R and T Kailath ESPRIT Estimation of signal parameters via rotational invariance techniques IEEE Trans ASSP Vol 37
100. ce algorithm avec bvec armaqs y p q avec bvec armaqs y p q norder maxlag samp_seg overlap flag armaqs estimates the ARMA parameters of the ARMA p q process y usingthe q slice algorithm The AR parameters are estimated via the function arrcest then the impulse response is estimated and finally the MA parameters are estimated TheAR order p must be greater than zero function maest should beused in the p 0 case the MA order q must also be specified Variablenorder specifies the cumulant order to be used and should be 3 or 4 the default value is 3 max ag specifies the maximum number of cumulant lags to be used its default valueiSptq Variablessamp seg overlap and lag control the manner in which sample cumulants are esti mated samp seg specifies the number of samples per segment the default valueis the length of the time series overlap specifies the percentage overlap between segments the allowed range is 0 99 and the default value is O Ifflag is biased then biased sample estimates are computed default if the first letter is not b unbiased estimates are computed If y isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset to zero andsamp seg is set tothe row dimension Cumulants are estimated from each realization and then averaged The estimated AR and MA parameters are returned in the vectors avec and bvec respectively armaq
101. ces can be estimated viabi coher andbi coherx The theoretical bispectrum of a linear process can be computed via bi spect trispect computes slices of the theoretical trispectrum of a linear process M fileg stat implements tests for non Gaussianity actually nonzero bispectrum and for linearity 1 45 1 Tutorial 1 46 Three routines are available for nonparametric estimation of the impulse response of a linear non Gaussian process mat ul implements the M atsuoka U Irych method given the bispectrum of the data it estimates the amplitude and phase of thetransfer function separately bi ceps usesthe bicepstral method it uses sample estimates of the third order cumulant bi cepsf isa frequency domain implementation of bi ceps and is useful if the cepstra do not decay fast enough Several algorithms are available for estimating the ARM A parameters of non Gaussian linear processes The ARMA orders can be estimated via arorder and maor der Routinearrcest can be used to estimate the AR parameters of AR or ARMA processes based on second third and fourth order cumulants The esti mates based on cumulants of different orders may not bethe same if the process is nonlinear or has inherent all pass factors or is noisy The parameters of an MA process contaminated by white noise may be estimated via routinemaest ARMA parameters can be estimated via the residual time series algorithm in armarts theMA estimation part of this thr
102. ct that a second order AR model may be adequate to dassify the data We used C ase Studies arrcest Withpz2 g20 norder 22 max ag 24tofit AR 2 models to successive segments of the data nsamp 1024 indz 1 nsamp for k l length x nsamp arx k arrcest x ind 2 0 2 24 ary k arrcest y ind 2 0 2 24 ind ind nsamp end plot arx 2 arx 3 x ary 2 ary 3 0 grid In Figure 1 42 the AR coefficient a 2 is plotted against the coefficient a 1 circles and crosses correspond to the two time series note that the data appear to be well separated in the AR coefficient domain The relationship between a 2 and a 1 appears to be linear we usedt s to estimate the slope and intercept from which we derived the line separating the two classes Since the autocorrelation is adequate for classifying the data higher order statistics are not useful in this example In other examples perhaps when the data are contaminated by lot of Gaussian noise AR estimates based on third or fourth order cumulants may be more useful 0 9 x T T T 0 8 0 7 F Class A 0 6 0 5 0 4 Class B 0 3 0 2r 0 1r 4 a me 0 1 1 1 1 L 1 2 1 8 1 6 1 4 1 2 1 0 8 Figure 1 42 Acoustic Data AR 2 Coefficients 1 121 1 Tutorial 1 122 Laughter Data Speech signals are highly nonstationary and are known to exhibit phase and frequency coupling phenomena Here we will use some of the H
103. ction problem 1 47 beamformer 1 64 1 65 1 77 bibliography 1 139 bi ceps 2 19 bi cepsf 2 21 bicepstrum 1 38 bi coher 2 23 bicoherence 1 4 auto 2 23 cross 2 25 estimation 1 20 2 23 bicoherx 2 25 bispecd 2 27 bispecdx 2 29 bispeci 2 31 bispect 2 33 bispectrum 1 8 cross 2 29 direct estimate 1 19 direct method 2 27 estimation 1 17 indirect method 1 18 2 31 theoretical 2 33 Wigner 2 90 2 92 Burg s maximum entropy estimator 1 74 C Canadian lynx data 1 114 Capon ML 1 74 Capon s maxi mum likelihood estimator 1 73 Choi Williams distributions 1 89 filter 1 92 1 95 1 100 smoothing 2 88 cross bicoherence 1 10 cross bi periodogram 1 17 cross bispectra Volterra systems 1 80 cross bispectrum 1 9 direct esti mate 1 18 direct method 2 29 estimation 1 13 indirect estimate 1 16 cross cumulant 1 81 Index cum2x 2 34 cum3x 2 36 cum4x 2 38 cumest 2 40 cumtrue 2 42 cumulants 1 4 auto 2 40 definitions 1 6 fourth order 2 40 2 42 sample esti mates 1 12 second order 2 34 2 40 third order 2 36 2 41 true 2 42 D demos 2 54 DOA 1 62 1 64 AR 1 74 beamformer 1 74 Capon ML 1 74 cumulant based estimators 1 74 eigenvector 1 74 ESPRIT 1 74 fourth order cumulants 2 44 minimum norm 1 74 MUSIC 1 74 Pisarenko 1 74 spatial covariance matrix 2 44 doa 2 44 doal mat 1 135 doagen 2 46 E eda 1 112 1 114 1 120 eigenvector method DOA 2 44 harmonic retrieval 2 50 eigenvector methods 1 67 ESPRIT 1 70 1 74
104. d and theoretical interquartile ranges are close to one another glstat z 0 51 256 Test statistic for Gaussianity is 12640 0657 with df 48 Pfa 0 Linearity test R estimated 606 9323 lambda 492 5759 R theory 59 9088 N 14 z was obtained by passing x through a nonlinearity z n 2x n hence z is non Gaussian and nonsymmetric We reject the Gaussianity assumption since the Pfa is small we cannot accept the linearity hypothesis since the estimated interquartile ranges is much larger than the theoretical value glstat l 0 51 256 Test statistic for Gaussianity is 49 931 with df 48 Pfa 0 3965 Linearity test R estimated 2 6047 lambda 1 8124 R theory 4 0038 N 14 1 25 1 Tutorial 1 26 The data in I arei i d and Laplace distributed symmetric Since the Pfa is high we cannot reject the Gaussian hypothesis but if the Gaussian hypothesis holds the bispectrum must be zero and we cannot conclude on the basis of the bispectrum alone whether or not the process is linear hence the results of the linearity test should be ignored in this case Parametric Estimators ARMA Models So far we have looked at nonparametric estimators Parametric estimators are often useful either because they lead to parsimonious estimates or because the underlying physics of the problem suggest a parametric model The basic idea is that if x n depends upon a finite set of parameters 6 then al
105. d as the location of the peak of Ry m A three point interpolation is used to improve the delay estimate 1 Ryy d 1 Ry d 1 d d __ gt __ 2Ryy d 1 2 Ry Cd Ryy d 1 tde tdeb 1 Krauss T J N Little and L Shure MATLAB Signal Processing Toolbox User s Guide The MathWorks Inc 1994 2 Nikias C L and R Pan Time delay estimation in unknown Gaussian spatially correlated noise IEEE Trans ASSP Vol 36 pp 1706 14 Nov 1988 2 83 tis Purpose Syntax Description Algorithm Reference Total Least Squares solution to a system of linear equations x flag tls A b tI s obtains the Total Least Squares TL S solution to the linear system of equations Ax b IfA is m by n and b is m by k it is required that m gt n k that is in the usual case whereb is a vector we have k 1 so that the system of equations should be overdetermined If the problem does not havea unique solution f ag will beset to unity TLS is useful when both A andb have errors F or example sample esti mates of cumulants have errors in them due to finite record lengths The least squares solution x to the set of linear equations Ax b minimizes the norm of the residual vector r Ax b subject to the condition that b r is in the range space of A the implicit assumption is that A is clean but b is noisy The TL Squares solution 1 assumes that both 4 and b may be n
106. default value is the power of 2 equal to or just greater than twice the signal length sigma is the parameter in the Choi Williams window and controls the amount of cross term suppression very large values result in no suppression whereas very small values result in loss of signal terms The default value is 0 05 flag ifflag is nonzero and x is real valued the analytic form of x is used instead of x the default valueis 1 If x t is a real valued signal with Hilbert transform y t then theanalyticsignal is defined as the complex valued signal z t x t jy t wx isthe computed WS Therows correspond totime samples and the columns to frequencies waxis isa vector whose entries are the frequencies associated with the columns of the WS Let the instantaneous autocorrelation be defined by r m n x n m x n m where n is identified with time and m with lag The ambiguity function AF is then given by AF m 0 r m m expc j 229 m The AF is multiplied by the Choi Williams window function w m 0 7 exp m6 i g ma 2 wig2c See Also References The2 D FT of AF m 0 w m 0 0 ton and mto f yields the filtered WS W f n The original signal must be sampled at twice the Nyquist rate or faster Wi 92 wi g3 wi g3c wi g4 wi g4c 1 Hlaswatch F and G F Boudreaux Bartels Linear and Quadratic Time F requency Representations IEEE Signal Processing Magazine pp 21 67 Apr 1992 2 Cohen
107. der statistics we should test the data for Gaussianity and if applicable for linearity as well Notethat the symmetry of the univariate pdf does not imply that the bispectrum is zero Model parameters e g ARMA harmonics estimated by algorithms based on cumulants of different orders need not agree with one another F or example in the case of thelaughter data we saw that the second and fourth order spectra were different a possible reason could be that some of the harmonic components can be well modeled as narrow band Gaussian Similarly in the case of ARMA modeling algorithms based on cumulants of different orders may lead to different results Consistent results may be expected only in the case of noise free data and long data lengths Consider the following three component model 3 y n Y x n i l x n Y hitloui n IK i l Assume that uj n isi i d Gaussian u n isi i d single sided exponential and u3 n is i i d and double sided exponential Laplace further assume that the ui s are independent of one another Let o denote the variance and y3 and Y4 the skewness and kurtosis of ui n Then 2 2 2 2 2 2 S gt 03 H 1 o5 H 6 o3lH 3 S3x 1 05 3 2H 2 1 H 2 2 H 2 2 Ya 3H 3 4 H 3 2 H 3 3 H 4 4 gt 03 1 131 1 Tutorial 1 132 In this example the power spectrum is the sum of the power spectra of the three components since u4 n and u3 n are both
108. domain input output relationship of the second order Volterra system under consideration is Y f H f X f Y Q fy f XO X3 fitfoa f It is assumed that Q f1 f2 Q f2 f1 Q f1 f2 Transfer functions H f and Q f4 f2 are obtained as the least squares solution to the set of equations Y m 2 A b m where bim tm of 32 2 033 23 m S and a e FS rne T Here M isthe FFT lengthnfft andm 0 M 2 ni tick 2 65 nlpow Reference 1 Powers E J Ch P Ritz C K An S B Kim R W Miksad and S W Nam Applications of digital polyspectral analysis to nonlinear systems modeling and nonlinear wave phenomena Proc Workshop on Higher Order Spectral Analysis pp 73 77 Vail Colorado J une 1989 nitick Purpose Syntax Description Nonparametric second order Volterra System Identification M ethod for Gaussian input signals h q nltick x y nfft wind samp seg overlap Implements Tick s nonparametric method for the identification of a second order Volterra system given both inputs and outputs The inputs are assumed to be Gaussian x istheinput tothe nonlinear system y isthe output process it must have the same dimensions asx Ify isa matrix columns are assumed to correspond to independent realizations nf ft isthe FFT length for computing power spectra and cross bispectra the default value is the power of 2 greater than the length of the time series row dime
109. e default settings are used If the function is invoked ast degen def aul t wherethe variabledef aul t may take on any value s then the default settings are used The colored noise sequence gi n is obtained by passing an i i d sequence Laplace uniform or normally distributed through an ARMA filter The colored noise sequence g n is obtained by passing the noise gj n through another ARMA filter Notethat the two noise sequences are correlated with each other Noise free signals can be generated by specifying a noise variance of zero when prompted Vectorss1 ands2 contain the simulated signals at the two sensors The vectorss1 ands 2 in thefiletde1 mat can be regenerated via s1 52 7t degen 1 2 81 tder Purpose Syntax Description Algorithm 2 82 Time delay estimation TDE using the ML windowed cross correlation method delay rxy tder x y max delay samp seg overlap nfft The time delay between two signals possibly corrupted by colored Gaussian noise is estimated using the Maximum Likelihood ML windowed cross correlation method x and y should both be vectors or matrices with identical dimensions If they are matrices each column is assumed to be a different realization These are the signals at the two sensors max delay istheabsolute value of the maximum expected delay samp seg specifies the number of samples per segment the default valueis the power of 2 just greater
110. e is 3 nlags is the maximum number of cumulant lags to be computed The default value of nl ags iSq p If fourth order cumulants are requested that is nor der is 4 then k refers to the third lag of the fourth order cumulant Ca i j k Its default value is 0 cmat is the vector or matrix of theoretical cumulants If norder iS 2 cmat isa column vector of length 2 nl ags 1 and consists of C2 m m nlags nlags If norder is30r4 cmat isa2 nlags 1 by2 nlags 1 matrix If norder is 3 the i j element of the matrixis Ca i nl ags 1 j nl ags 1 and ifnor der is 4 the i j element of the matrix is C4 i nl ags 1 j nlags Lk Note that the nl ags 1 nl ags 1 element of the matrix c mat contains C3 0 0 or C4 0 0 k cumtrue Algorithm See Also Reference Let h n denote the impulse response of a linear system excited by an i i d process u n with kth order cumulant yku Then the kth order cumulant of the output of the linear system is given by the Bartlett Brillinger Rosenblatt formula CC gt Tk 1 Yku L h n h n e c h n y 4 n 0 In this module Yku is assumed to be unity bispect trispect 1 Mendel J M Tutorial on higher order statistics spectra in signal processing and system theory Theoretical results and some applications Proc IEEE Vol 79 pp 278 305 1991 2 43 doa Purpose Syntax Description 2 44 Direction of arrival estimation based on sp
111. e constant Here y n is the observed process and z n is the IV when z n y n the RIV reduces to RLS Examples load riv ar 1 rivtr y 2 ar 2 rivtr y 3 hi ar 3 srivtr y 4 ar 4 rivtr zw 4 ar 5 rivtr zc 4 Now let us look at the five steady state AR parameter estimates disp ar 1 000 1 000 1 000 1 000 1 000 1 4789 1 5075 1 4931 1 5354 1 4686 0 7600 0 7910 0 7867 0 8248 0 7684 The data used in this example correspond to an AR 2 model with AR parameter vector 1 1 5 0 8 The vector y contains the noiseless signal Additive white Gaussian noise was added to y to obtain the noisy signal zw with a SNR of 10 dB Colored Gaussian noise generated by passing a white Gaussian sequencethrough the AR filter 1 0 0 49 was added to y to obtain zc also at a SNR of 10 dB 1 57 1 Tutorial 1 58 RIV Algorithm Double Lattice Form An order and time recursive solution of 1 94 and 1 95 leads to the double lattice algorithm and is given below Starting with n 1 evaluate 1 106 through 1 117 for m 21 2 M where M is the final desired order then repeat with n 2 and soon b i ES B du cda qup ime uS if n 1 Ym 1 1 1 106 b n 1 fm 1 n aP Nn Maman p mA ma Ym 1 D 1 107 f Am 100 Pen f m Bry Q n 1 m 1 1 108 b Am 100 Pim E TW m 1 1 109 Fh Q0 u n T aes m 1 n 1 110 fmn fg 4 0 Tg mN bm n 1 1 111 bm N bm 1 N 1 amp T
112. e final estimate of the power spectrum bispecd bi speci 1 Subba Rao T and M Gabr An Introduction to Bispectral Analysis and Bilinear TimeSeries Models pp 42 43 New York Springer Verlag 1984 2 Nikias C L and M R Raghuveer Bispectrum estimation A digital signal processing framework Proc IEEE Vol 75 pp 869 91 J uly 1987 bicoherx Purpose Syntax Description Cross bicoherence esti mation using the direct FF T based method bicx waxis bicoherx w x y bicx waxis bicoherx w x y nfft wind samp seg overlap The cross bicoherence of the three processes w x and y is estimated via the direct F FT based method w x y should all have the same dimensions nf ft specifies the FFT length to be used for computing the cross bi coherence the nominal default value is 128 ifnfft is smaller than samp seg the power of 2 just larger than samp seg will beused wi nd specifies the time domain window to be used it should be a vector of length samp seg By default thehanni ng window is used Data segments are multiplied by the time domain window then F ourier transformed to compute the frequency domain double and triple products The F ourier transform of the window function should be real and nonnegative samp seg spedifies thenumber of samples per segment Thedefault valueis set such that eight possibly overlapped records are obtained overlap Specifies the percentage overlap between s
113. e know that these harmonics are at integer multiples of a fundamental In this case it suffices to examine the diagonal slice of the bispectrum that is shown in Figure 1 49 We used pi ckpeak to locate the dominant peaks which were found to be at 0 0664 0 1289 and 0 1953 Hz given the FFT length of 256 we can conclude that the these three harmonics are quadratically frequency coupled the frequency coupling is also evident but less obvious from the power spectral estimates and the spectrogram bispecd reshape x 100 14 256 1 T T 0 4r 1 03r 4 9 0 27 ed J e e 0 1r J o g e S Op 4 e 5 2 amp 0 1r 4 2 2 N 0 24 e 4 LJ 03r J 0 4 4 0 5 l l 1 i 1 L i l i 0 5 0 4 0 3 0 2 0 1 0 0 1 0 2 0 3 0 4 Figure 1 48 bispectrum nonparametric 1 127 1 Tutorial 1 128 diagonal slice 0 014 T T T 0 012r 4 0 008 4 0 006 F 4 0 004 F 4 Js 0 4 0 3 0 2 01 0 0 1 0 2 0 3 0 4 0 5 Figure 1 49 Diagonal Slice of bispectrum We used routineqpct or to confirm the quadratic frequency coupling we let nlags 25 segsamp 100 overlap 0 the singular values shown in Figure 1 50 appear to fall in four clusters we used p 10 to estimate the parametric bispectrum shown in the bottom panel The dominant peak is at 0 0625 0 0625 Hz indicating the presence of a second harmonic at 0 1250 Hz also note that the contour plot indicates peaks at fy kf k 1 2 3
114. e moon laws of musical harmony N ewton s spectral decomposition of light 1677 Bernouilli 1738 and Euler s 1755 analysis of vibrating membranes and Prony s approximation for vibrating mechanisms 1793 M odern Fourier Analysis as weknow it today received its foundations in the work of Fourier 1807 although the roots of the Fast Fourier Transform FFT can betraced back to Gauss s work on orbital mechanics 1805 We will assume without loss of generality that the processes or signals of interest to us are zero mean We will also assume that the processes are discrete time with a sampling interval of T 1 corresponding to a normalized sampling frequency of 1 Hz sothat the N yquist frequency is 0 5 Hz The power spectrum is the primary tool of signal processing and algorithms for estimating the power spectrum have found applications in areas such as radar sonar seismic biomedical communications and speech signal processing Analog equipment to estimate the spectrum namely the spectrum analyzer has been around for more than five decades and may be found in almost any lab Our toolbox not only offers a substitute for that equipment it expands the analyst s toolkit to include algorithms more sophisticated than the simple conventional spectral analysis techniques Polyspectra and Linear Processes The usefulness of the power spectrum arises from an important theorem known as Wold s decompositi on which states that any
115. e the slice f4 f2 Note that the original signal should be sampled at twice the Nyquist rate in order to avoid aliasing Also note that the frequency axes are scaled by the factor of 2 3 in this routine weundothe scaling sothat you can find the peaks at the expected frequencies wi g2 wi g2c wi g3c wi g4 wi g4c wig3 References 1 Fonollosa J R and C L Nikias Wigner Higher Order Moment Spectra Definitions Properties Computation and Applications to Transient Signal Detection IEEE Trans SP J an 1993 2 Gerr N L Introducing a third order Wigner distribution Proc IEEE pp 290 92 Mar 1988 3 Swami A Third order Wigner distributions Proc CASSP 91 pp 3081 84 Toronto Canada M ay 1991 2 91 wig3c Purpose Syntax Description Algorithm 2 92 Computes the diagonal slice of the Wigner bispectrum with Choi Williams filtering wx waxis wig3c x nfft sigma flag Computes the diagonal slice of the Wigner bispectrum WB the Choi Williams filter is applied in order to suppress the cross terms in the WB of multicomponent signals x is the signal and must be a vector nfft specifies the FFT length and hence the frequency resolution this parameter must be greater than threetimes thelength of the signal x in order to avoid aliasing The default value is the power of 2 equal to or just greater than thrice the signal length sigma isthe parameter in the Choi Will
116. e two peaks in the power spectrum The resolution of a power spectral density estimator is the smallest valueof f4 f2 that leads to two discernible peaks The resolution of the classical power spectrum estimators e g Welch and Blackman Tukey is of order 1 M where M isthe effective window length n contrast the resolution of the periodogram is 1 N and since M lt N the periodogram exhibits higher resolution than the Blackman Tukey and Welch estimators The choice of the window function dictates the resolution variance tradeoff if W f has a broad main lobe the power spectrum estimate will be smoother the variance of the estimate will be smaller and the resolution will be lower It should be noted that zero padding the data or the correlation sequence does not improve the resolution it only makes the spectral estimate denser this may lead to improved accuracy in estimating the locations of well separated peaks 1 65 1 Tutorial 1 66 AR and ARMA Models The basic idea in the parametric estimators is that the noiseless harmonic process in 1 124 obeys a self driving AR p model that is 28 p Y a k x n k 0 k 0 1 131 2nfk where a 0 1 and the roots of the A z polynomial are given by z d k 1 p Hence we can use any standard algorithm to estimate the AR parameters the frequencies can be determi ned from the roots of the esti mated AR polynomial The determination of p is now equivalent
117. ecord The default value is set such that eight possibly overlapped records are obtained overlap specifies the percentage overlap between segments The default value is 0 The allowed range is 0 99 Ify isa matrix the columns are assumed to correspond to independent realizations or records in this caseover lap is set to zero and samp seg is set tothe row dimension bi c is the estimated bicoherence it is an nf f t by nf ft array with origin at the center and axes pointing down and tothe right waxi s istheset of frequencies associated with the bicoherenceinbi c Thus the ith row or column of bi c corresponds to the frequency waxis i i 1 nfft Frequencies are normalized that is the sampling frequeng is assumed to be unity A contour plot of the magnitude of the esti mated bicoherence is displayed The data y are segmented into possibly overlapping records the mean is removed from each record the time domain window is applied and the FFT computed the bispectrum of the kth record is computed as 2 23 bicoher See Also References 2 24 B m n X m X n X m n where X denotes the FFT of the kth record and the spectrum is computed as P m X amp m 2 The spectral and bispectral estimates are averaged across records and the bicoherence is then estimated as BG s fD bi A SS a ON ICCf 4 f 2 PO PGOPGSU f where B f4 f2 is the final estimate of the bispectrum and P f is th
118. ed polyspectra of x n are zero but the even ordered polyspectra such as thetrispectrum are not identically equal to zero 1 23 1 Tutorial 1 24 Note also that univariate symmetry or zero skewness does not mean zero bispectrum In the above example suppose that y3 0 and yk h gt k 0 then C3 0 0 Y3u _k h gt k 0 so that process x n has zero skewness This does not imply that C3 t1 t2 0 or S3 1 2 0 In particular let h 0 h 1 1 h k 0 Vk gt 1 then x n u n u n 1 is easily seen to be symmetrically distributed around zero but the bispectrum of x n is not identically zero In general we use notions of skew and symmetry in the context of random variables or univariate pdfs bispectra and trispectra relate to random processes multivariate pdfs Examples load gldat glstat g 0 51 256 Test statistic for Gaussianity is 22 179 with df 48 Pfa 0 9995 Linearity test R estimated 0 88819 lambda 0 6982 R theory 2 9288 N 14 The data in g are Gaussian distributed Since the Pfa is high we cannot reject the Gaussian hypothesis but if the Gaussian hypothesis holds the bispectrum must be zero and we cannot conclude on the basis of the bispectrum alone whether or not the process is linear hence the results of the linearity test should be ignored in this case In this case the data are Gaussian hence also linear glstat u 0 51 256 Test statistic for Gaussia
119. ee step algorithm works well only under good SNR conditions Routinear maqs implements the q slice algorithm it estimates both AR and MA parameters Inthis section we have given a quick overview of the area of higher order statistics for more extended tutorial expositions see 41 37 38 linear Prediction M odels Linear Prediction Models TheAR model is also obtained if one considers the problem of linear prediction LP Depending upon whether one chooses forward backward or forward backward prediction criteria different solutions are obtained The Levinson algorithm enables efficient esti mation of the parameters of an AR process from its autocorrelation AC sequence In the deterministic context AR modeling leads to the forward backwards least squares F BLS solution The RLS algorithm is a time recursive solution to the LP problem whereas the lattice algorithm provides a time and order recursive solution Levinson Recursion Consider the linear prediction problem for a stationary process x n In the forward prediction problem we want to choose fap k I _ to minimize the forward prediction error variance Pp Po lain 1 63 p e m Y ajdxn k a0 1 k 0 1 64 where the subscript p denotes the order of the predictor Similarly in the backward prediction problem we want to choose Cp Ku _ 1 tominimize the backward prediction error variance Pp Pp e 1 65 p e Y cy K oX n K cy p
120. egments The default value is 50 Theallowed range is 0 99 If w x y are matrices the columns are assumed to correspond to independent realizations in this caseoverlap isset tozero andsamp seg is set totherow dimension bi cx is the estimated cross bicoherence it is an nf f t by nf f t array with origin at the center and axes pointing down and tothe right waxis is the set of frequencies associated with the cross bicoherence in bi cx Thus the ith row or column of bi cx corresponds tothe frequency waxi s i izl nfft Frequencies are normalized that is the sampling frequency is assumed to be unity A contour plot of the magnitude of the esti mated cross bicoherence is displayed 2 25 bicoherx Algorithm See Also References 2 26 The data w x y are segmented into possibly overlapping records the mean is removed from each record the ti me domain window is applied and theFFT computed the cross bispectrum of the kth record is computed as B wxy k m n Wyi m Xy n Y M n where Wy Xy and Yp denote the FFT of the kth segments of w x and y The spectra are computed as P y m Wy m a Pxk m XI m 2 and Py k m Yy m The spectral and cross bispectral estimates are averaged across records and the cross bicoherence is then estimated as 2 B wxyCf 1 f2 bi QUE EE T E EEUU EEUU KEEN EDS IC f4 f2 PG DPXRGOPVG fa where Bwxy f1 f2 is the averaged estimate of the cross bi
121. ence 1 61 represents a set of linear equations in thelog magnitudeln H f Itis shown in 36 that a full rank set of linear equations can be obtained by appropriate choice of f and fo The phase relationship in 1 62 holds only modulo 2z If the unwrapped 2 D phase were available then 1 62 yields a linear set of equations in the phase of the transfer function 0 f This algorithm is implemented in routine mat ul where we also incorporate the phase unwrapping algorithms of 54 Because of the phase ambiguity this routine is not recommended for routine use Examples cmatezcumtrue 1 3 5 1 5 1 3 5 bsp fft2 flipud cmat 64 64 hest matul bsp You should see the display in Figure 1 9 Here we computed the true bispectrum of a MA 2 process by evaluating the FT of its true cumulant sequence We then used the M atsuoka Ulrych algorithm to estimate the impulse response Note the scale ambiguity including sign Estimated impulse response 0 6 T T T T 1 L L L L L L L 40 30 20 10 0 10 20 30 40 sample number Figure 1 9 IR Estimated via the Matsuoka Ulrych Algorithm matul 1 42 Polyspectra and Linear Processes Linear Processes Theoretical Cumulants and Polyspectra If x n is an ARMA process that is p q x n Y x n k Y b k u n k k l k 0 then its impulse response can be calculated via h n p h n I9 b n where h n 0 for n 0 h 0 1 and b n 0 i
122. epstrum is nonzero only along thek 1 axes and the main diagonal k 1 k 1 bm m 9 2 Y Amp m l 1 nzl nzl1 k 1 x hemp 3m mp 1 2 k 1 InYku 90m 1 2 Tekalp and Erdem 69 showed a process whose polycepstrum exists is linear if and only if its polycepstrum has the above k line region of support Based on this they have also proposed measures of linearity 1 40 Polyspectra and Linear Processes Examples load mal hest ceps bicepsf y 8 Y ou should see the display in Figure 1 8 Note that the complex cepstrum decays rapidly to zero ThetrueMA parameters are 1 0 9 0 385 0 771 note the scale ambiguity including sign in the estimated impulse response complex cepstrum T 1 5 ah 0 5 4 or M 0 5 L L L L i 1 80 60 40 20 20 40 60 80 impulse response 3 T 2r zl 1E J or 4 1L 2 L L L L L L 1 80 60 40 20 o 20 40 60 80 Figure 1 8 Cepstrum and IR Estimated by bicepsf The Matsuoka Ulrych Algorithm We can write the bispectrum in terms of its magnitude and phase as S3x fa f2 2 M fif expGe f1 f2 Let H f H f exp j6 so that INM fp f2 Inlrsy InlIHG In H G2 In B G3 F 1 61 D f 0 0 f5 0 d 2 fy f2 9 9 2 f f2 mod 27 1 62 1 41 1 Tutorial where f4 f take on discretized values on a grid The assumption that y3y 1 will only introduce an overall scalar ambiguity h
123. ero if m gt q When the true cumulants are replaced by sample estimates the estimated values of C3y q i 0 i 20 will not beidentically zero a statistical test is used to determine whether the estimated values are close to zero This test is based on estimating the theoretical variance of the sample estimates of C3 m 0 Sample estimates C3y m 0 and their variances are estimated for m ranging from qmin to qmax which reflect our a priori knowledge of the bounds on the true order q For an MA q process the asymptotic variance of the sample esti mate of C3y q 1 0 can be estimated via 20 2041 2 y 2 l y DyG q 1 C3 q 1 0 j q x Iy G Dy j q 1 3 q 1 0 A N c q 1 o where N is thelength of the time series The sample esti mates are asymptotically Gaussian and unbiased hence the threshold t in Pr e5 m 1 0 lt t m 1 1 pfa is given by tm 1 erfinv 1 pfa 42o m 1 whereerfinv isthe MATLAB inverse error function in practice we use the sample estimate o Let m denote the largest value of m in the range qmin to qmax for which c3 m 1 0 gt t m 1 so that the hypothesis of MA m model fails then the estimated order is q 2 mg 1 if such an m does not exist the MA order is declared to be qmax 1 This algorithm is implemented in routine maor der This is a statistical test and pfa specifies the fraction of the time that the test results will
124. es are scaled by the factor of 1 2 which can be easily fixed The WT also suffers from the presence cross terms As in the case of WS we can attempt to suppress the cross terms by appropriate filtering Define the smoothing kernel 2 2 2 2 P O 14 T2 T3 exp 0 t4 t5 13 0 1 186 The filtered WT is then given by Won fy f f3 f dede dc dcsdu gr mo fate fata This algorithm is implemented in routine wi g4c where the slice f1 f2 f3 is computed the resulting sliced WT is real valued As in the case of the WS cross terms can usually be suppressed by this approach but with a concomitant loss of resolution In contrast to the WB the filtering does not unduly distort the auto terms 14 15 Time Frequency Distributions Examples clf load wi gdat subplot 221 wig4 s2 0 subplot 222 wig4 s2 Subpl ot 223 wig4 s3 Subpl ot 224 wig4 s4 You should see the display in Figure 1 24 Signal s2 isa harmonic with a nominal frequency of 0 05 Hz which has been multiplied by a Gaussian window and is centered at n 20 Note that the WT of the signal is concentrated around its nominal center frequency note also the interaction between the energies at the positive and negative frequencies giving risetothe interference terms around D C By usingtheanalytic version of the signal the negative frequency terms are suppressed this also suppresses the distortion terms around D C Signal s3 is a harmo
125. es each where K N M x O 100 M The kth record or segment of x consists of the samples y i 2y i k 1 M1 i 1 M k 1 K The sample mean is removed from the kth record and sample estimates of the cumulants are computed For example sample estimates of third order cumulants are obtained as M max 0 m n 1 C3 m n y YkDYk m y i N M m n i 1 MaxX 0 m n where M m n M for biased estimates and M max 0 m n min 0 m n for unbiased estimates The final estimate is given by Second and fourth order K 1 C3 m n Y Cs m n k 1 cumulant estimates are obtained similarly see also the algorithm descriptions for M files cum2x and andcum4x cumtrue cum2x cum3x cum4x 1 Nikias C L and M R Raghuveer Bispectrum estimation A digital signal processing framework Proc IEEE Vol 75 pp 869 91 J uly 1987 2 41 cumtrue Purpose Syntax Description 2 42 Computes theoretical that is true cumulants of a linear process cmat cumtrue ma cmat cumtrue ma ar norder nlags k cumt rue computes the theoretical that is true cumulants of a linear ARMA process ma is the MA parameter vector and must be specified I fq isthe MA order then ma Will be of length q 1 ar isthe AR parameter vector its default value is 1 Ifp isthe AR order then ar will be of length p 1 norder isthe cumulant order allowed values are 2 3 and 4 the default valu
126. es pointing down and to theright waxis is the set of frequencies associated with the bispectrum in bs pec thus the ith row or column of bs pec corresponds to the frequency waxis i The sampling frequency is assumed to be unity Let H f B f A f denote the transfer function of the ARMA filter then the bispectrum is given by B fif2 H CF H G2H f f2 cumtrue trispect 1 Subba Rao T and M Gabr An Introduction to Bispectral Analysis and Bilinear Time Series M odds pp 42 43 New York Springer Verlag 1984 2 Nikias C L and M R Raghuveer Bispectrum estimation A digital signal processing framework Proc IEEE Vol 75 pp 869 91 J uly 1987 2 33 cum2x Purpose Syntax Description 2 34 Computes the cross cumulant covariance of two signals cvec cum2x Xx y maxlag sSamp_seg overlap flag Computes the second order cross cumulant covariance of the two signals x andy x y should have identical dimensions maxlag specifies the maximum lag of the cumulant to be computed its default valueis O samp seg Specifies the number of samples per segment Its default valueisthe row dimension of y if y is a row vector the column dimension is used as the default value overlap specifies the percentage overlap between segments the allowed range is 0 99 the default valueis O Ifflag is biased then biased sample estimates are computed default if the first letter is not b unbiased e
127. escribe the temporal evolution of the spectrum or the polyspectrum moment spectrum are useful tools A time frequency distribution TF D is a transform that maps a 1 D signal into a 2 D time frequency map which describes how the spectral content of the data evolves with time As such TFD s are the natural tools for the analysis synthesis interpretation and processing of nonstationary signals Among the more well known TF D s are the short time Fourier transform STFT 1 the Gabor representation 17 and the Wavelet transform 33 55 Linearity is a desirable property of algorithms used in analyzing linear systems however quadratictime frequency distributions have been proposed analyzed and interpreted as time varying power spectra 13 25 The so called Cohen class of shift invariant distributions includes special cases such as the Spectrogram Rihaczek Page Wigner Ville distribution WD and Choi Williams distributions 13 25 The WD has become quite popular because it possesses a number of useful properties and has been used in the analysis of phase modulated signals which are common in radar and sonar 3 An important property of the WD is that every member of Cohen s class can beinterpreted as a 2 D linearly filtered version of the WD 25 A third order WD was introduced in 18 it was 1 89 1 Tutorial generalized and its properties were studied in 14 15 60 61 62 Higher order WD s describe the evolu
128. f ne O q For alinear process x n x n yk h lgu n K where u n is i i d cumulants and polyspectra of orders 2 3 4 are given by 1 16 through 1 18 and 1 19 through 1 21 We can estimate the true cumulants using routinecumest the theoretical bispectrum via routine bi spect and slices of the theoretical trispectrum via routinetrispect Examples First we will compute the theoretical third order cumulants of an ARMA 2 1 process with AR parameters 1 1 5 0 8 and MA parameters 1 2 wewill then display the estimates using MATLAB s functions mesh or contour cmat cumtrue 1 2 1 1 5 0 8 3 25 clf subplot 121 mesh 25 25 25 25 cmat grid on subplot 122 contour 25 25 25 25 cmat 8 grid on You should seethe display in Figure 1 10 1 Tutorial 1 44 SES SSN OSS COLD OC Oo 2 V Ss SS SAV SAG o N o lj N o Q 50 40 20 10 0 10 20 Figure 1 10 True Third Order Cumulants of an ARMA 2 1 Process cumtrue Let us now compute the theoretical bispectrum of another ARM A 2 1 process ma 1 2 ar 1 0 8 0 65 bisp bispect ma ar 128 Y ou should seethe display on Figure 1 11 The 12 dotted lines emanating from the origin divide the bifrequency domain of the bispectrum into 12 regions of symmetry 0 4 0 3 02 oal 2 0 o ab Soal 0 3 0 4 F 0 5 Ll f i 0 5 0 4 0 3 0 2 0 1 Figure 1 11 Bispectrum
129. f non Gaussian processes from high order moments IEEE Trans on Automatic Control Vol 35 1 pp 27 37 1990 17 Gabor D Theory of communication J i e E Vol 93 pp 429 59 1946 18 Gerr N L Introducing a third order Wigner distribution Proc IEEE pp 290 92 Mar 1988 19 Giannakis G B and J M Mendel Identification of non minimum phase systems using higher order statistics IEEE Trans ASSP Vol 37 pp 360 77 Mar 1989 20 Giannakis G B and M Mendel Cumulant based order determination of non Gaussian ARMA models IEEE Trans ASSP pp 1411 21 Aug 1990 21 Golub G H and C F Van Loan Matrix Computations pp 420 5 Baltimore TheJ ohns Hopkins University Press 1983 22 Hasan T Nonlinear time series regression for a class of amplitude modulated sinusoids J Time Series Analysis Vol 3 pp 2 109 22 1982 23 Haykin S Adaptive Filter Theory New J ersey Prentice Hall pp 391 94 198 205 464 70 234 36 2nd ed 1991 24 Hinich M J Testing for Gaussianity and linearity of a stationary time series J Time Series Analysis Vol 3 pp 169 76 1982 References 25 Hlaswatch F and G F Boudreaux Bartels Linear and Quadratic Time Frequency Representations IEEE Signal Processing Magazine pp 21 67 Apr 1992 26 ohnson D H The application of spectrum estimation methods to bearing estimation problems Proc IEEE Vol 70 pp 975 89 1982 27
130. f the function is invoked ashar mgen default where the variable def aul t may take on any value s then the default settings are used see Examples below Each column of z mat corresponds to a different realization The matrixz mat in thefilenar m mat can be regenerated via zmat harmgen 1 2 53 hosademo Purpose Syntax Description 2 54 A guided tour of the Higher Order Spectral Analysis Toolbox hosademo hosademo takes you on a guided tour of the Higher Order Spectral Analysis Toolbox hosademo also offers a brief introduction to the area of higher order statistics hosahelp Purpose Gives a one line synopsis for all M files in the Higher Order Spectral Analysis Toolbox Syntax hosahelp Description hosahel p gives a one line synopsis for all the documented M files in the Higher Order Spectral Analysis Toolbox 2 55 hprony Purpose Syntax Description Algorithm Reference 2 56 Estimates the parameters of a complex transient signal modeled as the sum of complex exponentials with decaying amplitudes a theta alpha fr hprony x p Models a complex transient signal as the sum of complex exponentials with decaying amplitudes p xm Y adoe P expincoco j2rf k k 1 where p is the order the a k s are amplitudes 6 k s are the initial phases a k s are the damping factors and f k s are the frequencies x is the time series it must be a vector p is the
131. fined by y n Y h k x n k Y ack Dx n k x n 1 k 0 k 0l 0 1 157 The corresponding frequency domain representation is Y f H f X f Y Q fy XCEOXG fitfo f 1 158 and is obtained by Fourier Transforming both sides of 1 157 The time domain product term x n k x n I leads to convolution in the frequency domain which is represented by the condition f f f gt It is usually assumed that the quadratic kernel is real and symmetric that is k l q I k or equivalently Q fzf2 Q f2 f1 Q fy f2 It is readily verified that Q f4 f2 in the region f f1 0 lt f4 lt 1 4 specifies Q f1 f2 everywhere The basic problem is given x n and y n n 1 N we want to estimate the linear part h k and the quadratic part q k I Solution Using Cross Bispectra In 70 an algorithm to estimate the parameters of the model in 1 157 was developed under the assumption that x n is stationary Gaussian From 1 157 we can see that the cross spectrum is given by Syx f ECY X f H PIELX A X f H A S AA 1 159 The term involving the triple product of the X f s disappears since x n is symmetrically distributed Since S f and Sxy f can be estimated we can esti mate H f the linear part N onlinear Processes Consider the cross cumulant Cyxy p ELx n x n x y n p 1 160 E amp x n x n 7 Y h k x n p k a k 20 Y alk Dxin p xin p
132. he cross bicoherence Examples load qpc dbspec bicoher zmat 128 Bicoherence estimated via the direct FFT method T T T T T T T 0 4 4 0 3 4 AS ose 4 O O 0 1 E N or J 0 1 0 2r 4 O 0 3 E 4 0 4 4 oE i i i i i 0 5 0 4 0 3 0 2 0 1 o 0 1 0 2 0 3 0 4 fi Figure 1 4 Bicoherence Estimate bicoher A contour plot of the estimated bicoherence is shown in Figure 1 4 The data consists of quadratically phase coupled harmonics with frequencies at 0 10 Polyspectra and Linear Processes 0 15 and 0 25 Hz and an uncoupled harmonicat 0 40 Hz The maximum value of the bicoherenceis 0 6442 the valueislessthan unity because of the additive noise which affects the power spectrum estimate load nll bicx bicoherx x x y You should see the display in Figure 1 5 The cross bicoherence is significantly nonzero and nonconstant indicating a nonlinear relationship between x and y The nonsharpness of the peaks as well as the presence of structure around the origin indicates that the nonlinear relationship is not purely of the form y n x n and that x y are not narrow band processes From the description of nl 1 mat weseethaty isthe output of a second order Volterra system whose input x is Gaussian Cross Bicoherence T O A g g j a E o 3r m o2 f2 o T j 1 1 1 1 1 1 1 1 0 5 0 4 0 3 0 2 0 1 o 0 1 0 2 0 3 0 4
133. he unit circle In this case the cepstrum of H z is well defined and is given by 42 ho df exp j2afk INH f 1 55 The bicepstrum is also well defined 43 and is given by A j2nfym j2nxf bonny adie e Ins gura 1 56 h m 8 m h n 8 n h m 8 m n Iny4 6 m n 1 57 The complex cepstrum h k istheinverse FT IFT of thelog of the FT of h n hence h n isthe IFT of the exponential of the FT of the complex cepstrum other techniques to obtain h n from h k are described in 42 Notice that for a linear process the bicepstrum is nonzero only along the axes m 0 n O and the diagonal line m n Direct implementation of 1 56 demands 2 D phase unwrapping Pan and Nikias 43 have developed an alternative method based on the relationship Y Y mim n Cay k m I n kC3y k 1 aS 1 58 which taking into account 1 57 reduces to Y mh m C3y k m n C3y k m 1 m kKCay k 1 m eo The complex cepstrum is known to be exponentially bounded hence one may replace theinfinite summation over m by a finite summation m q tom p Polyspectra and Linear Processes where p and q may be based on some a priori knowledge notethat the p and q have nothing todo with the orders of an ARM A p q process Now we have a set of linear equations for estimating a finite set of p q parameters p y mh m C3y k m n C3 k m 1 m kC3y k D mosd 1 59 Note that the above equation does no
134. her assumes that u n has unit variance or that b 0 1 Instead of esti mating P f Vfe 1 2 1 2 as in the nonparametric approach we have to estimate only p q 1 parameters namely a k 1 b l9 d j and o2 Let h n denote the impulse response of the model in 1 31 hence H z B z A z TheAR and MA parameters arerelated to the impulse response IR via n L a k h n k b n ne 0 q k 0 1 35 0 n O q 1 36 n practice the observed process is noisy that is y n x n w n 1 37 where process w n is additive colored Gaussian noise the color of the noise is usually not known Given the noisy observed data y n we want to estimate the a k s and the b k s in 1 31 We assume that the model orders p and q are known The determination of the model orders p and qis an important issue and is discussed later routinesarorder and maorder implement AR and MA model order determination techniques In general it is better to overestimate model 1 27 1 Tutorial 1 28 orders rather than to underestimate them however specific implementations or algorithms may suffer from zero divided by zero type problems Motivation to use cumulants for the ARMA parameter estimation problem is as follows Correlation based methods can be used successfully only if q 20 pure AR and the additive noise is white Even in the noiseless case correlation based methods cannot identify inherent al
135. hm section for details is computed and you are prompted to choose an order for the SVD low rank approximation Thesingular values of the cumulant matrix associated with the TDE problem rarely show a clean break in contrast with thoseassociated with theharmonicretrieval problem solved byhar mest ortheDOA problem sol ved by doa it is recommended that the order chosen for the SVD low rank approximation be greater than the value of max del ay del ay is the estimated delay of signal y with respect to signal x avec is the estimated parameter vector it is a vector of length 2 by max delay 1 and corresponds totime samples max delay max delay The estimated delay is the time sample at which avec attains its maximum absolute value 2 77 tde Algorithm See Also Reference 2 78 Let x n s n wy4 n and y n 2As n D w2 n denote the two observed signals where A is the amplitude gain D is the signal delay between the two sensors and w1 n and w2 n are assumed to be jointly Gaussian Let P be the maximum expected delay Then P y n Y avi x n i w n i P where a n 0 n z D and a D A Let Cy p E y n x n 7 x n p and Cyxx t P E fx n x n 7 x n p We obtain the third order recursion P Cyxx T p Ye a i Cd Gi pi i P Using this equation for various values of p and t we get a system of linear equations in the a i s Cy ya Cy The estimated delay is the index
136. ial 38 Nikias C L and J M Mendel Signal processing with higher order spectra IEEE Signal Processing Magazine Vol 10 No 3 pp 10 37 J uly 1993 39 Nikias C L and R Pan Time delay estimation in unknown Gaussian spatially correlated noise IEEE Trans ASSP Vol 36 pp 1706 14 Nov 1988 40 Nikias C L and A Petropulu Higher Order Spectra Analysis A Nonlinear Signal Processing Framework New J ersey Prentice Hall 1993 41 Nikias C L and M R Raghuveer Bispectrum estimation A digital signal processing framework Proc IEEE Vol 75 pp 869 91 J uly 1987 42 Oppenheim A V and R W Schafer Digital Signal Processing New J ersey Prentice Hall 1989 43 Pan R and C L Nikias The complex cepstrum of higher order cumulants and non minimum phase system identification IEEE Trans ASSP Vol 36 pp 186 205 F eb 1988 44 Pan R and C L Nikias Harmonic decomposition methods in cumulant domains Proc CASSP 88 pp 2356 59 New York 1988 45 Patel J K and C B Read Handbook of the Normal Distribution Sec 7 9 4 M Dekker New York 1982 46 Perryin F and R Frost A unified definition for the discrete time discrete frequency and discrete ti me frequency Wigner distributions IEEE Trans ASSP Vol 34 pp 858 67 Aug 1986 47 Pflug L A G E loup J W loupandR L Field Properties of higher order correlations and spectra for bandlimited transients J Acoust So
137. ial correlation matrix Periodogram Beamformer The autocorrelation matrix Ry with k l entry Ry I k E amp n x n k I which appears in the harmonic retrieval problem will be replaced by the Harmonic Processes and DO A spatial correlation matrix R E amp n x n wherex n is thenth snapshot in the DOA problem The FT of the spatial temporal correlation sequence is called the beamformer periodogram The autocorrelation sequence of the signal in 1 125 is given by p Ryy t Y jou expj 2ntf Ry t kaal 1 130 which follows because the harmonics are uncorrelated phases are mutually independent Note that the autocorrelation sequence is not absolutely summable however the FT of Ryy t exists provided we are willing to use the Dirac delta functions defined by g f 3p f 9 0 5 p 17 In principle one can estimate the frequencies by picking peaks in the esti mated power spectrum the variance of this frequency estimator is of order 1 N 22 where N is the length of the time series The power spectrum is limited by its poor frequency resolution that is frequency separation must be larger than 1 N In practice this severely limits the applicability of the power spectrum as a direct estimator of the frequency Synthetics for the two problems can be generated via har mgen and doagen Resolution and Variance If x n consists of two complex harmonics at frequencies f4 and fz we expect to se
138. iams window and controls the amount of cross term suppression very large values result in no suppression whereas very small values result in loss of signal terms The default valueis 0 5 flag ifflag is nonzero and x is real valued the analytic form of x is used instead of x the default valueis 1 If x t is a real valued signal with Hilbert transform y t then theanalyticsignal is defined as the complex valued signal z t x t jy t wx is the computed WB The rows correspond toti me samples and the columns to frequencies waxis isa vector whose entries are the frequencies associated with the columns of the WB Define the instantaneous triple product r3 t tz t2 X t at at2 x t Bt at2 x t at Bt2 where a 1 3 and B 2 3 Define the smoothing kernel 0 11 T2 exp e r r5 sigma wig3c See Also References The filtered WB is then given by WON fy f deds dz du e t 13 airt nue r3 t T1 Tz P 9 14 T2 In this routine we compute the slice f1 f2 Note that the original signal should be sampled at twice the N yquist rate in order to avoid aliasing Also note that the frequency axes are scaled by the factor of 2 3 In this routine we undo the scaling sothat you can find the peaks at the expected frequencies NOTE The application of the Choi Williams kernel to the Wigner bispectrum does not guarantee preservation of the auto terms consequently this routine should
139. igher Order Spectral Analysis Toolbox M files to illustrate these features We will use the initial 1400 sample segment of the data in thefile aught er mat which is included in the standard MATLAB distribution package The commands that weused to generatethe various figures and to estimate various quantities are listed in Table 1 4 Figure 1 43 shows the data and the histogram theunivariate distribution does not appear to be symmetrical The mean standard deviation skewness and kurtosis were esti mated to be 0 5621 536 69 0 1681 and 1 3277 indicating that the data are non Gaussian and that the univariate pdf is not symmetric distributed laughter data za L L 1 L L L L 0 20 40 60 80 100 120 140 160 180 time in ms 100r E 0 1 f dears 4 3 2 1 0 1 2 3 4 5 Figure 1 43 Laughter Data C ase Studies Figure 1 44 shows the spectrogram or the STFT computed viaspecgram we used an FFT length of 512 Hanning window of length 256 and an overlap of 240 samples Three dominant frequency tracks can be seen in the figure approximately around 550 Hz 1100 Hz and 1550 Hz thelast formant begins around 30 ms additional fragments are visible around 1800 Hz and 2100 Hz specgram x 512 8192 256 240 4000 J 3500 F 4 3000 F 4 m a Q o T i
140. inting down and to the right waxis isthe set of frequencies associated with the bispectrum in bs pec Thus the ith row or column of bs pec corresponds to the frequency waxi s i Frequencies are normalized that is the sampling frequeng is assumed to be unity The data y are segmented into possibly overlapping records biased or unbiased sample estimates of third order cumulants are computed for each record and then averaged across records a lag window is applied tothe esti mated cumulants and the bispectrum is obtained as the 2 D FFT of the windowed cumulant function bispecd 1 Subba Rao T and M Gabr An Introduction to Bispectral Analysis and Bilinear TimeSeries Models pp 42 43 New York Springer Verlag 1984 2 Nikias C L and M R Raghuveer Bispectrum estimation A digital signal processing framework Proc IEEE Vol 75 pp 869 91 J uly 1987 bispect Purpose Syntax Description Algorithm See Also References Theoretical bispectrum of an ARMA process bspec waxis bispect ma ar nfft The theoretical bispectrum corresponding to an ARMA process is computed ma is the MA parameter vector and must be specified ar is the AR parameter vector its default value is 1 0 nfft specifies the FFT length to be used for computing the bispectrum the default value is 512 bspec is the bispectrum corresponding to the ARMA model It is an nf ft by nfft array with origin at the center and ax
141. introducing the ideas of bias and consistency The bias of an estimator is defined as E Sy s the estimate is said to be unbiased if the bias is zero that is E Sy s Often this holds true only as N ee in which case the estimate is said to be asymptoti cally unbiased The bias by itself does not completely characterize the estimate If the esti mate is good we expect that Su will take on values around the true quantity s The natural measure of the spread is the squared deviation around thetrue quantity s E Sn s k The estimate is said to be asymptotically consistent if the squared deviation goes to zero as N e This condition is sometimes called mean square consistency A consistent esti mate is necessarily asymptotically unbiased Estimating Cumulants In practice we have a finite amount of data xn _ gt and we must obtain consistent estimates of cumulants The sample estimates are given by N xy k T y x myin k n N 1 22 M Myy k N Y xmyn k amp ionis 1 23 EA Cx D qz Y xam k z n 1 n N 1 24 Polyspectra and Linear Processes N Cwxyz k l m M y w n x n k y n Dz n m n N ii wx k Cyz l m we cess m M wz m M xy 1 25 where N4 and N are chosen such that the summations involve only x n s with ne 0 N 1 unbiased estimates are obtained if N3 is set equal tothe actual number of terms which are averaged for example AN
142. ith fixed f3 f3 Itisannfft by nfft array with origin at the center and axes pointing down and tothe right waxis is the set of frequencies associated with thet spec thus theith row or column of t spec corresponds to the frequency waxi s i Frequencies are normalized that is the sampling frequency is assumed to be unity Let H f B f A f denote the transfer function of the ARMA filter then the trispectrum is given by S3 f if 2 f3 HOH G2 H 3H f3 f2 f3 In this routine frequency f3 is fixed at the specified value of f 3 bispect cumtrue 1 Subba Rao T and M Gabr An Introduction to Bispectral Analysis and Bilinear Time Series Models pp 42 43 New York Springer Verlag 1984 2 Nikias C L and M R Raghuveer Bispectrum estimation a digital signal processing framework Proc IEEE Vol 75 pp 869 91 J uly 1987 Purpose Syntax Description Algorithm See Also References Computes the Wigner spectrum wx waxis wig2 x nfft flag Estimates the Wigner spectrum WS of a signal using the time domain approach x is the signal and must be a vector nfft specifies the FFT length and hence the frequency resolution this parameter must be greater than twice the length of the signal x in order to avoid aliasing The default valueis the power of 2 equal to or just greater than twice the signal length flag ifflag is nonzero and x is real valued the analytic form of x is used i
143. its nominal center frequency In the last panel we compute the WB of the sum signal s 452 3 note the cross terms in the WS some of these can be eliminated by using routinewi g3c A direct consequence of the definition of the WB isthat the frequency axes are scaled by the factor 2 3 4 5 n this routine this scaling has been undone so that the peaks appear at the expected frequencies Examples load wi gdat wi g3c 54 256 0 2 1 Comparethedisplay in Figure 1 23 with the 2 2 panel in Figure 1 20 and with Figure 1 21 notice the suppression of cross terms CWB f1 f2 s 0 2 60r 4 time in samples L L L 1 1 0 25 0 2 0 15 0 1 0 05 0 0 05 0 1 0 15 0 2 frequency Figure 1 23 Wigner bispectrum With Choi Williams Smoothing wig3c 1 97 1 Tutorial 1 98 Wigner Trispectrum Define the instantaneous fourth order product r t t4 To 13 X t T X t T T1 X t T T2 X t T 13 1 183 where t T1 T 13 4 The WT is then given by i2 Win fq f gt f3 Joridrodeze WM Tek V Ta 03 1 184 Notice that two terms are conjugated also note the sign of the f 315 term This algorithm is implemented in routinewi g4 where the slice f1 f2 f3 f is computed that is we compute j2 W n f dci dzodsse cL e Ti T gt 13 1 185 Note that the original signal should be sampled at twice the Nyquist rate in order to avoid aliasing Also note that the frequency ax
144. iversity of New York at Buffalo in 1980 and 1982 Currently heis a Professor of Electrical Engineering at USC in Los Angeles where heis also Director of CRASP and Associate Dean of Academic Research Dr Nikias is a Fellow of the IEEE Heis the author of over 150 journal and conference papers two textbooks a monograph and four patents Heis co author of the textbook Higher Order Spectral Analysis A Nonlinear Signal Processing Framework Prentice Hall Inc 1993 Tutorial 1 Tutorial Introduction This section of the User s Guidedescribes how to begin using the Higher Order Spectral Analysis Toolbox for your signal processing applications t assumes familiarity with basic MATLAB gt as well as a basic understanding of signals and systems Thereis much more information in a stochastic non Gaussian or deterministic signal than is conveyed by its autocorrelation or power spectrum Higher order spectra which are defined in terms of the higher order moments or cumulants of a signal contain this additional information The Higher Order Spectral Analysis HOSA Toolbox provides comprehensive higher order spectral analysis capabilities for signal processing applications The toolbox is an excellent resource for the advanced researcher and the practicing engineer as well as the novice student who wants to learn about concepts and algorithms in statistical signal processing The Higher Order Spectral Analysis Toolbox is a collectio
145. l Equation Approach Proc IEEE Vol 70 pp 907 38 1982 Capon J High resolution frequency wavenumber analysis Proc IEEE Vol 57 pp 1408 18 1969 Chandran V and S Elgar A general procedure for the derivation of principal domains of higher order spectra IEEE Trans on Signal Processing SP 42 pp 229 33 J an 1994 10 Choi H and WJ Williams Improved Time F requency Representation of Multicomponent Signals Using Exponential Kernels IEEE Trans ASSP pp 862 71 J une 1989 11 Classen T A C M and W F G Mecklenbra ker TheWigner distribution a tool for time frequency signal analysis Phillips Res Vol 35 pp 217 50 276 300 372 89 1067 72 1980 1 139 1 Tutorial 22 Classen T A C M and W F G Mecklenbratker The aliasing problem in discrete time Wigner distributions IEEE Trans ASSP pp 1067 72 Oct 1983 13 Cohen L Time Frequency Distributions A Review Proc IEEE pp 941 81 J uly 1989 14 Fonollosa J R and C L Nikias Wigner Higher Order Moment Specra Definitions Properties Computation and Applications to Transient Signal Detection IEEE Trans SP J an 1993 15 Fonollosa J R and C L Nikias Analysis of finite energy signals using higher order moments and spectra based ti me frequency distributions Signal Processing Vol 36 pp 315 28 1994 16 Friedlander B and B Porat Asymptotically optimal esti mation of MA and ARMA parameters o
146. l of its statistics can be expressed in terms of 0 For example we obtain parametric estimates of the power spectrum by first estimating 6 and then evaluating P f 0 The specific form we postulate for the relationship between 0 and the sequence x n constitutes a model A popular model in time series analysis is the Auto Regressive Moving Average ARMA model p q x n Y a k x n k Y b k u n k k 1 k 0 1 31 where u n is assumed to be an i i d sequence with variance c The Auto Regressive AR polynomial is defined by p A z Y a k z k 0 1 32 where a 0 1 A z is assumed to have all its roots inside the unit circle that is A Z 0 gt Z lt 1 this condition is also referred to as the minimum phase or causal and stable condition In general no restrictions need to be placed on the zeros of the Moving Average MA polynomial d k B z Y b k z k 0 1 33 Polyspectra and Linear Processes however B z is usually assumed to be minimum phase The minimum phase assumption is usually not true for discrete time processes which are obtained by sampling a continuous time process algorithms based on H OS do not require this assumption The power spectrum of the ARMA process is given by ICON c d P9 IA cr 1 34 Note that the power spectrum does not retain any phase information about the transfer function H z B z A z Since we do not have access to the sequence u n one eit
147. l pass factors and cannot resolve the phase of the system Non Gaussian processes are not completely characterized by their Second order statistics by using higher order statistics we are exploiting more of the information contained in the data A note of caution for odd ordered cumulants yy z 0 does not imply that the kth order polyspectrum can be used to reconstruct H f For example even if Y3u 0 S3x f 4 f 2 in 1 20 may be identically equal to zero This can happen when H f isa relatively narrow band bandpass signal explicit conditions are given in 63 An even more trivial example is y4 E u n lt 0 but H 0 0 in which case y x 0 For polyspectra of even order say 2k Si EP RO lS a cannot be identically zero if yz 0 63 Thetransfer function H z B z A z is said to have an inherent all pass factor if a root of B z lies at 1 2 where z is a root of A z The power spectrum is blind toall pass factors If the data are noise free and if the ARM A model does not have any inherent all pass factors then techniques based on the autocorrelation power spectrum may be used to estimate the model parameters OncetheARMA parameters have been esti mated M ATLAB s routinef r eqz can be used to estimate the transfer function and theoretical spectrum Higher Order Spectral Analysis Toolbox routinesbi spect andtrispect can be used to compute the theoretical bispectrum and slices of the theoretical trispectrum c
148. le p Y ap k C3 m k p 0 m q 1 q p nt 1 75 In both cases the resulting matrix is Toeplitz but not symmetric An efficient order recursive solution to the system of equations is given by the Trench recursion 74 Let R denote the nonsymmetric Toeplitz matrix whose i j entry is p i j The recursion is given by m 1 Al i12 Y aq Qd0p m I9 Kee 1 76 m 1 Am 1 Y Sm_ 1 k p 1 k k 0 1 77 1 49 1 Tutorial 1 50 f rf Am 1 m7 P me 1 78 b b Am 1 Im 5 e 1 79 f wb Pm Pode En dose a_ m rf m p 1 81 b c m Tr Qucm 1 82 f a k 2 a 4 K I c q k 1 kz1l m 1 m m 1 m m 1 1 83 b C k 2c k 1 rma k k 1 m 1 m m 1f m m 1f 1 84 with initial conditions Pg p 0 ag 0 1 and co 0 1 The inversion algorithm requires that the matrix 9t be strongly nonsingular that is R as well as all of its principal minors are nonsingular The condition of nonsingularity is weaker than the usual assumption of positi ve definiteness In the symmetric positive definite case the T s are bounded by unity and the Pas are nonincreasing and Cm k a m k k Z0 m In the non Hermitian case theT s are not bounded and the P s may increase or decrease This algorithm is implemented in trench Examples Let us compute the autocorrelation sequence for an AR 4 model with AR parameters 1 0 3 0 1 0 39 0 72 and then apply the Trench recur
149. lem and nl gen can beused to compute the output of a second order Volterra system Time Frequency Distributions Time Frequency Distributions The popular and natural tool in linear system analysis is the Fourier transform which decomposes a signal into its frequency components The power spectrum or energy spectrum gives us information about thefrequency components in the data but not about the temporal localization of these components Alternatively the original time signal itself gives us good time resolution obviously but does not give us any idea about frequency localization except in trivial cases So far we have defined spectra bispectra and trispectra only for stationary processes as such they are inappropriate for the analysis of nonstationary processes or transient signals F or nonstationary processes we can define time varying cumulants and polyspectra in the same way that one talks about time varying covariances and power spectra Consider the signal y t s t g t where s t is a deterministic signal and g t is stationary zero mean noise perhaps Gaussian The process y t is nonstationary its mean is s t and for k gt 1 Cky ag thus only the first order cumulant carries information about the signal s t In contrast the higher order moments of y t depend both upon s t and g t thus in this case it is preferable to use higher order moments 63 For such signals time frequency distributions that d
150. lihood algorithm The conventional estimators are easy to understand and easy to implement but are limited by their resolving power the ability to separate two closely spaced harmonics particularly when the number of samples is small F or random signals these estimators typically require long observation intervals in order to achieve acceptably low values for the variances of the estimate AN The natural estimator of the power spectrum is the FT of Rxx m N po j2 E POPE MGE Rime 5 y xckye 24 N m N 1 k 0 This estimator also known as the periodogram can be computed as the squared magnitude of an N point FFT of the observed time series Since E Rxx m R xxm the periodogram is an unbiased estimator of P f However the periodogram is not a consistent estimator because var Py4 f Py f that is its variance does not goto zero as N gt 56 1 15 1 Tutorial 1 16 If x n is Gaussian white noise with variance o2 its power spectrum should be flat Py f o the variance of the periodogram esti mate P f is given by A N 4 var Px f 6 Ee sin 2xf Note that the variance does not gotozero as N ee that is the estimateis not consistent The covariance between esti mates at frequencies f k N and f2 m N is zero thus on the one hand as N increases the variance at any f does not go to zero however the spacing between estimates that are uncorrelated decreases as 1 N as a con
151. lled Choi Williams filter which can help suppress the cross terms in the WS WB and WT It should be noted that this filtering may in fact destroy signal terms in the WB but not in WS or WT Wigner spectra of different orders offer different perspectives on the data and are useful for exploratory data analysis Time Delay Estimation Time Delay Estimation Thetime delay esti mation problem occurs in various applications for example in the determination of range and bearing in radar and sonar It also has esoteric applications such as the measurement of the temperature of a molten alloy by measuring the passage time of a signal Other applications include analysis of EEG data The basic model is as follows two sensors record delayed replicas of a signal in the presence of noise x t s t w t 1 187 y t As t D w t 1 188 where D is the delay of the signal at the y sensor relative to the signal at the x sensor A is the relative amplitude gain and w t and wy t are sensor noises Given x t y t t 0 N 1 we want to estimate the delay D The basic idea is to shift the signal y t and compare the shifted waveform with x t the best match occurs when the shift equals the delay D This notion is made more precise by using the cross correlation between the two signals We assume that s t is a stationary process and that the noises are zero mean A Cross Correlation Based Method The cross correla
152. logram third order 1 106 hosademo 2 54 hosahelp 2 55 hprony 2 56 instrumental variables 2 57 ivcal 2 57 K kurtosis 1 7 L Levinson Durbin recursion 1 48 2 85 linear models frequency domain bicepstral method 2 21 lag domain bicepstral method 2 19 linear prediction 1 31 1 47 adaptive 1 54 linear processes impulse response esti mation 1 37 theoretical cumulants 1 43 theoretical pol yspectra 1 43 linearity test 1 24 linearity tests 2 47 M MA models 1 29 order estimation 2 61 parameters estimation 2 58 mal mat 1 135 maest 2 58 maorder 2 61 Matsuoka Ulrych algorithm 2 63 matul 2 63 minimum phase 1 11 minimum norm method DOA 2 44 harmonic retrieval 2 50 1 3 Index 1 4 mixed phase 1 11 1 135 ML Capon 1 74 MUSIC 1 68 DOA 2 44 harmonic retrieval 2 50 N nl 1 mat 1 135 nl2 mat 1 136 nlgen 2 64 nl gen 2 64 nl pow 2 65 nltick 2 67 noise subspace 1 69 nonredundant region 1 133 normal equations cumulant based 2 17 deterministic 1 53 P peak picking 2 69 periodogram 1 15 1 64 phase coupling 1 84 1 132 pi ckpeak 2 69 Pisarenko s method 1 67 DOA 2 44 harmonic retrieval 2 50 pitfalls 1 131 polycepstra methods 1 38 polycepstral methods 1 38 polycepstrum 1 40 polyspectra linear processes 1 4 windows 1 16 polyspectrum definitions 1 6 power spectrum Wigner 2 87 power spectrum estimation 1 4 conventional methods 1 15 criterion based estimators 1 72 criterion based methods 1
153. lues and these are identical with the eigenvalues of B A The singular B case is discussed in 21 It is shown in 30 57 that the generalized eigenvalues of Cy Cyz are given by Ak exp j2nAdsin 0 A k 1 p Once A has been esti mated we can readily obtain the 6 s This algorithm for a uniformly spaced linear array is implemented in routinedoa the generalized eigenvalues are computed via the MATLAB built in function ei g Algorithms based on fourth order cumulants are discussed next Criterion Based Estimators The FT at frequency wo can be viewed as a filter that passes only the components with frequency og while suppressing the rest With finite data lengths however the spectral estimate suffers from the problems of sidelobe leakage with details dictated by the choice of window function Consider the length p finite impulse response FIR filter p y n J a k x n k a 1 151 where the input x n is assumed to be zero mean The variance of the filter output is given by Harmonic Processes and DO A 2 2 H Oo E n a R a y ly n Xx 1 152 where a a 0 a p and Rxx is the p 1 x p 1 autocorrelation matrix We would like the output variance o to be a measure of the power spectral density of x n at frequency fs Alternatively the FIR filter should pass a sinusoid at frequency fo with unity gain hence we have the constraint p Y a k exp j2nkf 1 E
154. lues of the AR and MA orders Let 1 50 i ca m C m p Ca m q ay SE mM 1 51 c4y m C4y m p 0 Cas y m d 01 1 52 Then the singular values of the matrix where k 2 k 3 or k 24 are computed Polyspectra and Linear Processes Let s m denotethe singular values If the cumulant estimates are perfect we expect that exactly p of the singular values will be nonzero with sample estimates we expect p dominant singular values The AR order pis C Cky q p 1 c y q 2p 1 then given by the value of n which maximizes s n s n 1 that is it corresponds to the index at which the singular values show the maximum drop This is a nonstatistical test and is based on the fact that for an ARMA p q model only p of the singular values should be nonzero This algorithm is implemented in routinearorder Examples load armal p arorder y 3 You should seethe display in Figure 1 6 The estimated AR order is 2 The time series y corresponds to a non Gaussian ARMA 2 1 process contaminated by AWGN with SNR of 20 dB The order determination is based on third order cumulants difference in singular values cumulant order 3 T T T T T 0 6 Figure 1 6 Estimate of AR Order Using SVD Method arorder 1 35 1 Tutorial 1 36 MA Order Determination The basic idea is that for an MA q process the true values of the cumulant C3y m 0 will be identically z
155. mercial software documentation was obtained as set forth in subparagraph a of the Rights in Commercial Computer Software or Commercial Software Documentation Clause at DFARS 227 7202 3 therefore the rights set forth herein shall apply and b For any other unit or agency NOTICE Notwithstanding any other lease or license agreement that may pertain to or accompany the delivery of the computer software and accompanying documentation the rights of the Government regarding its use reproduction and disdo sure are as set forth in Clause 52 227 19 c 2 of the FAR MATLAB Simulink Handle Graphics and Real Time Workshop are registered trademarks and Statefl ow and Target Language Compiler are trademarks of The Math Works I nc Other product or brand names are trademarks or registered trademarks of their respective holders Printing History May 1993 First printing September 1995 Second printing J anuary 1998 Third printing About the Authors Tutorial l Introduction ssseses esses 1 2 Polyspectra and Linear Processes 05 1 4 tr oductlon s Seated Re pcm EEUU EE 1 4 Definitions Xlsl4desssteREIdeSRUN AS RE een ade dalek a 1 6 Why Do We Need Higher Order Statistics 1 10 Bias and Variance of an Estimator 1 11 Estimating Cumulants 00 000 1 12 Exatmplesz ccn EE ee Ce AR 1 14 Estimating Polyspectra and Cross polyspectra 1 15 Estima
156. n m C5y k C5 1 m C44 I C54 K m M m M Kk 1 2x 2x 1 7 where M 54 m E amp n x n m and equals C2 m for a real valued process The first order cumulant is the mean of the process and the second order cumulant is the autocovariance sequence Note that for complex processes there are several ways of defining cumulants depending upon which terms are conjugated The zero lag cumulants have special names C2 0 is the variance and is usually denoted by 65 C3 0 0 and C4 0 0 0 are usually denoted by y3 and Yax We will refer to the normalized quantities y3 65 as the skewness and Yax 02x as the kurtosis These normalized quantities are both shift and scale invariant If x n is symmetric distributed its skewness is necessarily zero but not vice versa if x n is Gaussian distributed its kurtosis is necessarily zero but not vice versa Often theterms skewness and kurtosis are used to refer to the unnormalized quantities y3 and yay If x n isan i i d process its cumulants are nonzero only at the origin If x n is statistically independent of y n and z n x n y n then C42 k l m C4 kl m Cay k l m with similar relationships holding for cumulants of all orders This additivity property simplifies cumulant based analysis 1 Tutorial 1 8 The cumulants of a stationary real valued process are symmetric in their arguments that is C5 k C k C3 k 1 C3x l K
157. n which maximizes a n A low rank approximation of the cumulant matrix C44 may be used tdeb tder 1 Nikias C L and R Pan Time delay estimation in unknown Gaussian spatially correlated noise IEEE Trans ASSP Vol 36 pp 1706 14 N ov 1988 tdeb Purpose Syntax Description Time delay estimation TDE using the conventional bispectrum method delay h tdeb x y max delay nfft wind samp seg overlap The time delay between two signals possibly corrupted by colored Gaussian noise is estimated using the conventional bispectrum method This involves computing the third order hologram described below that exhibits a peak at thetrue delay x and y should both be vectors or matrices with identical dimensions If they are matrices each column is assumed to be a different realization These are the signals at the two sensors The signals are assumed to have nonzero third order cumulants for example they cannot be symmetrically distributed max delay istheabsolute value of the maximum expected delay nf ft specifies the FFT size to be used the nominal default value is 128 the actual FFT size used will bemax samp seg nfft wi nd specifies the frequency domain smoothing window If wi nd is a scalar the Rao Gabr window 2 ry W m n Apne m n e G T N of length wi nd will be used here N is half the FFT length nf ft and G is the set of points m n satisfying _ 42 m n74 mn lt ne nfft 2
158. n of M files that implement a variety of advanced signal processing algorithms for spectral estimation polyspectral estimation and computation of time frequency distributions with applications such as parametric and nonparametric blind system identification time delay estimation harmonic retrieval direction of arrival estimation parameter estimation of Volterra nonlinear models and adaptive linear prediction Other potential applications include acoustics biomedicine econometrics exploration seismology nondestructive testing oceanography plasma physics radar sonar speech etc For the newcomer to the field of higher order statistics spectra some excellent starting places are T2 Mendel J M Tutorial on higher order statistics spectra in signal processing and system theory Theoretical results and some applications Proc IEEE Vol 79 pp 278 305 1991 T3 Nikias C L andJ M Mendel Signal processing with higher order spectra IEEE Signal Processing Magazine Vol 10 No3 pp 10 37 J uly 1993 T4 Nikias C L and A P Petropulu Higher Order Spectra Analysis A Nonlinear Signal Processing Framework New J ersey Prentice Hall 1993 T1 Nikias C L and M R Raghuveer Bispectrum estimation A digital signal processing framework Proc IEEE Vol 75 pp 869 91 J uly 1987 Introduction The field of higher order statistics spectra and its applications to various signal processing problems
159. ng frequency of 1 Hz Symmetry properties of the trispectrum include Sax fr fo f3 Sax fa f f Sax fo f f3 Say fa faf fafa Sax fy fof The literature is somewhat confusing both in the derivation as well as in the description of the nonredundant regions The nonredundant region for a continuous time band limited process with a Nyquist frequency of 0 5 Hz is the triangle with vertices 0 0 1 4 1 4 1 2 0 A tutorial treatment of the differences between the continuous time and the discrete time cases is given in 47 related discussions may be found in 63 The nonredundant region of the trispectrum is discussed in 6 9 47 Similar to the cross correlation we can also define cross cumulants for example C k l E x n y n k z n 1 xyz l x N y Z dua The cross bispectrum is defined by a ie j2xf k d2nfil S Up ey Uy Gdkhe Ue SS MESS 1 13 Note that the bispectrum S3 f1 f2 is a special case of the cross bispectrum obtained when x y z The cross bi coherenceis another useful statistic which is defined as 1 Tutorial 1 10 Sxyzf r f2 Diy Ty f2 0 M xyz r 7 2 dS2x f3 f295 f5240 2 1 14 The autobicoherenceis obtained when x y z M filebi coherx can be used to estimate the cross bicoherence and bi coher can beused to estimate the bicoherence The cross bicepstrum of three processes is defined by bm m f Inu fp fe g ap ar Te and is well defined only if Sy
160. nic with a nominal frequency of 0 15 Hz which has been multiplied by a Gaussian window and is centered at n 50 Notethat the WT of the signal is concentrated around its nominal center frequency In the last panel we compute the WS of the sum signal s 4 2 453 note the cross terms in the WT some of these can be eliminated by using routinewi g4c WT WT 60 60 50 50 o E 40 g40 5 E e 30 307 220 Ez 10 10 0 2 0 1 0 0 1 0 2 0 2 0 1 0 0 1 0 2 frequency frequency WT WT 60 60r g 50 g50 amp 40 amp 40 d 5 8 iE 30 x 2 E30p c FEE 2p 10 10 0 2 0 1 0 0 1 0 2 0 2 0 1 o 0 1 0 2 frequency frequency Figure 1 24 Wigner trispectra wig4 1 99 1 Tutorial 1 100 Examples load wigdat wig4c s4 Compare the display in Figure 1 25 with the 2 2 panel of Figure 1 24 and Figure 1 21 and Figure 1 23 note that the cross terms in Figure 1 24 have been suppressed CWT s 0 05 T T T T T T T 60 4 50r om 40r 4 o a E oO o amp 30 J o E 20r 4 10r z L f L L L L L L L 0 25 0 2 0 15 0 1 0 05 0 0 05 0 1 0 15 0 2 frequency Figure 1 25 Wigner Trispectrum With Choi Williams Filtering wig4c Summary The Higher Order Spectral Analysis Toolbox offers routines to compute the WS and slices of the WB and the WT and areimplemented in wi g2 wi g3 and wi g4 Routines wi g2c wi g3c and wi g4c implement the so ca
161. nity is 17 4885 with df 48 Pfa 1 Linearity test R estimated 0 72383 lambda 0 51704 R theory 2 7453 N 14 The data in u areuniformly distributed Since the Pfa is high we cannot reject the Gaussian hypothesis but if the Gaussian hypothesis holds the bispectrum must be zero and we cannot conclude on the basis of the bispectrum alone whether or not the process is linear hence the results of the linearity test should beignored in this case glstat e 0 51 256 Polyspectra and Linear Processes Test statistic for Gaussianity is 253 3529 with df 48 Pfa O Linearity test R estimated 7 8894 lambda 9 4555 R theory 8 4655 N 14 The data in e obey the single sided exponential distribution The Pfa in rejecting the Gaussian hypothesis is very small hence we are comfortable in accepting the hypothesis of non Gaussianity The estimated and theoretical interquartile ranges the R s are fairly close to one another hence we accept the linearity hypothesis The estimate of the interquartile range was based on N 14 samples glstat x 0 51 256 Test statistic for Gaussianity is 277 5194 with df 48 Pfa 0 Linearity test R estimated 6 7513 lambda 10 6519 R theory 8 968 N 14 x was obtained by passing through a linear filter hence it is linear and non Gaussian Wereject the Gaussianity assumption sincethe Pfa is small we accept the linearity hypothesis since the estimate
162. ns The Min Norm solution is not useful in this case A modification developed by Tugnait 72 appends the following sets of equations to the preceding set of 2q equations q amp 3C5y n b k C3y k n q C3y n q Kat 1 40 q 4C5y n b I9C4y k n q 0 C4y n q 0 kz1 1 41 1 29 1 Tutorial 1 30 where n q q N 0 Now we have 4q equations in 2q 1 unknowns hence we can obtain the least squares solution We solve either 1 38 and 1 40 or 1 39 and 1 41 Since both b k and b2 k or b3 k are estimated we need some method for combining the estimates Let b4 k and b k denote the estimates of b k and b2 k If all of the estimated b2 k s are nonnegative then the final MA parameter estimate is obtained as k signfb k 0 5 b3 k b5 k otherwise b k b k When fourth order cumulants are used the method estimates both b k and b3 k Let b Kk and b3 k denote the estimates of b k and b k If all the esti mated b3 k s have thesamesign as the corresponding b k s then the final MA parameter estimate is obtained as b signf b k t lo po 7172 otherwise b k b k This algorithm is implemented in routine maest If the observed process is z n y n g n where g n is additive noise independent of y n then Coz n C5y n C5g n Vn If g n is white C55 n 0 n 0 hence C5z n 2 Coy n n 0 In 1 40 and 1
163. ns at time n are given by see 1 87 and 1 88 linear Prediction M odels de cs Snan Pp 0 x a 1 94 n c n EH c i P N 7 7 1 95 where the deterministic correlation matrix n is given by n em Y xoxoa k p 1 1 96 where a 0 lt A lt 1 is the forgetting factor which is introduced so that the algorithm will be sensitive to the nonstationary environment We could solve 1 94 or 1 95 at each n but the computational load is very high Let the order p be fixed and let ap n 1 w n Then a time recursive solution is given by X P n 1 x n H0 7 13A x my Pc DX 1 97 a n x n amp 2 w n 1 x n 1 98 w n w n 1 k n a n 1 99 P n X IPi 15 X ki x n P n 1 1 100 with initial values w 0 20 and P 0 8 41 wheres is a small positive constant 1 55 1 Tutorial 1 56 k n is the gain vector a n is thea priori prediction error and P n l n This algorithm is known as the Recursive Least Squares Algorithm RLS and is a special case of the Recursive Instrumental Variable Algorithm RIV which is discussed next RIV Algorithm Transversal Form Consider the noisy AR process p y n x n w n x n Y a k u n k k 1 where u n is assumed to be zero mean i i d and independent of additive noise w n which may be colored The normal equations with m 1 p p Y a k E y n z n m k 0 m gt 0 k 0 1 101 whe
164. nsion for matrices wi nd specifies the frequency domain smoothing window If wi nd is a scalar the Rao Gabr window of length wi nd will be used This window is defined by 2 2 2 wim n Gpr manjeg T N where N is half the FFT length nf f t and G is the set of points m n satisfying wind 2 2 m n mn cz 5 nfft 2 e A unity value for wi nd results in no windowing e Ifwind lt 0 the default value of 5 will be used e f wi nd is a vector it is assumed to specify a 1 D window from which a 2 D window is computed W m n 7 w m w n w m n 1H2 e f wind is a 2 D matrix it is assumed to specify the 2 D smoother directly The bispectrum estimate averaged across records is smoothed by convolving with the 2 D window function The window function should be real and nonnegative 2 67 nitick Algorithm See Also Reference 2 68 samp seg specifies the number of samples per segment The default value is set such that eight possibly overlapped records are obtained overlap specifies the percentage overlap between segments The default value is 50 The allowed range is 0 99 h is the estimated impulse response of the linear part q is the estimated impulse response of the quadratic part The output of a second order Volterra system is given by co co y n Y h oxin o Y Y ak Dxin oxi 1 k 0 k 0l 0 The linear part h k is estimated in the frequency domain from Syx H Sy
165. nstead of x the default valueis 1 If x t is a real valued signal with Hilbert transform y t then the analytic signal is defined as the complex valued signal z t x t jy t wx isthe computed WS Therows correspond to time samples and the columns to frequencies waxis isa vector whose entries are the frequencies associated with the columns of the WS Let the instantaneous autocorrelation be defined by r m n 2x n m x n m wheren is identified with time and m with lag The WS is then given by W f n Y rim n exp jnfm m The original signal must be sampled at twice the Nyquist rate or faster in order to avoid aliasing wi g2c wi g3 wi g3c wi g4 wi g4c 1 Hlaswatch F and G F Boudreaux Bartels Linear and Quadratic Time F requency Representations IEEE Signal Processing Magazine pp 21 67 Apr 1992 2 Cohen L Time F requency Distributions A Review Proc IEEE pp 941 81 J uly 1989 2 87 wig2c Purpose Syntax Description Algorithm Computes the Wigner spectrum with Choi Williams filtering wx waxis wig2c x nfft sigma flag Estimates the Wigner spectrum WS of a signal eliminates cross terms in the WS of multi component signals by applying the Choi Williams filter x is the signal and must be a vector nfft specifies the FFT length and hence the frequency resolution this parameter must be greater than twice the length of the signal x in order to avoid aliasing The
166. o Reference 2 16 Let s m denote the singular values The AR order p is then given by the value of n which maximizes s n s n 1 that is it correspondstotheindex at which the singular values show the maximum drop maorder 1 Giannakis G B andJ M Mendel Cumulant based order determination of non Gaussian ARMA models IEEE Trans ASSP pp 1411 21 Aug 1990 arrcest Purpose Syntax Description Estimates AR parameters using the normal equations based on autocorrelations and or cumulants avec arrcest y p avec arrcest y p minlag norder maxlag samp seg overlap flag arrcest estimates the AR parameters of the ARMA p q process y via the normal equations based on autocorrelations and or cumulants p isthe AR order and must be specifi ed mi nl ag is the minimum value of m to be used in the normal equations If you want to estimate the AR parameters of an ARMA p q process mi nI ag should be greater than q Seethe Algorithm subsection for more details Theabsolute valueofnorder specifies the cumulant order to be used if nor der is negative least squares solutions based on the simultaneous solution of both autocorrelation and cumulant based normal equations are obtained The allowed values of nor der are 2 3 4 The default value is 2 maxlag specifies the maximum number of cumulant lags to be used its default valueis p mi nl ag Variablessamp seg overlap andflag control the manner in
167. o be computed its default valueis O samp seg Specifies the number of samples per segment Its default valueisthe row dimension of y if y is a row vector the column dimension is used as the default value overlap Spedifies the percentage overlap between segments the allowed range is 0 99 and the default value is O Iff I ag is biased then biased sample estimates are computed default if the first letter is not b unbiased estimates are computed Parametersk1 and k2 control which 1 D slice of the cross cumulant is computed the default value for both the parameters is 0 By varyingk1 andk2 we can obtain the entire fourth order cross cumulant cvec Will contain the sample estimates of cum w n x n m y n k1 Z n k2 m maxl ag maxlag where the superscript denotes complex conjugation and cum a b c d E abcd E ab E cd E ac E bd E ad E bc and it is assumed that a b c d are zero mean random variables Ifw x y z are matrices columns are assumed to correspond to different realizations in this caseoverlap isset to zero andsamp_seg is set tothe row dimension the cross cumulant is estimated from each realization and then averaged across the suite of realizations cum4x Algorithm See Also Let w n x n y n and z n n 1 N denote the three time series Let M denote the number of samples per segment the value of the variable samp_seg Let O denote the percentage overlap be
168. ociated process z n I m n and Th m Nn arethereflection coefficients for theforward backward prediction of y n and may beinterpreted as normalized cross correlations between the forward prediction error in one lattice and the backward prediction error in the other u n is the scale factor for converting the reflection coefficients associated with y n to those associated with z n F m N is the correlation between the mth order forward prediction errors fm Nn and fm N Bin is the correlation between the mth order backward prediction errors b n and ban Ym N is the scale factor for converting the a priori prediction errors totheir a posteriori values its absolute value is bounded by unity Aq n is the cross correlation between the a posteriori backward prediction error bj 4 n 1 and the a priori tilded forward prediction error ty 4 n and Aq 4 n is the cross correlation between fma n and oq n 1 Time update equations for F Nn and B n are gi ven by 1 115 and 1 116 This algorithm is implemented in routiner i vdl Examples load riv ar 1 2zri vdl y 2 ar 2 rivdl y 3 ar 3 rivdl y 4 ar 4 rivdl zw 4 ar 5 rivdl zc 4 Now let s see the five AR parameter estimates disp ar 1 000 1 000 1 000 1 000 1 000 1 4814 1 5106 1 4941 1 5358 1 4692 0 7623 0 7941 0 7876 0 8252 0 7690 linear Prediction M odels The data used in this example correspond to an AR 2 model wi
169. oes not guarantee preservation of the auto terms consequently routine wi g3c should be used with great caution 1 95 1 Tutorial 1 96 Examples cl f load wigdat Subpl ot 221 Subpl ot 222 subpl ot 223 Subpl ot 224 You should see the wig3 s2 0 wi g3 52 wi g3 53 wi g3 54 display in Figure 1 22 WB f1 f2 WB f1 f2 60 60 a 50 v 50 2 2 40 40 5 30 ES 30 B20 CE 20 10 10 0 2 0 1 0 0 1 0 2 0 2 0 1 o 0 1 0 2 frequency frequency WB f1 f2 WB f1 f2 60 60 v 50 v 50 2 o amp 40 amp 40 5 5 c 30 30 20 20 10 10 0 2 0 1 Q 0 1 0 2 0 2 0 1 0 0 1 0 2 frequency frequency Figure 1 22 Wigner bispectra wig3 Signal s2 is a harmonic with a nominal frequency of 0 05 Hz which has been multiplied by a Gaussian window and is centered at n 20 Note that the WB of the signal is concentrated around its nominal center frequency note also the interaction between the energies at the positive and negative frequencies giving risetothe interference terms around D C By usingtheanalytic version of thesignal the negativefrequency terms are suppressed this also suppresses the distortion terms around D C Time Frequency Distributions Signal s3 is a harmonic with a nominal frequency of 0 15 Hz which has been multiplied by a Gaussian window and is centered at n 50 Notethat the WB of the signal is concentrated around
170. of an ARMA 2 1 Process bispect Polyspectra and Linear Processes We will compute several slices of thetrispectrum and then display them as a movi e clf clear ma 1 2 ar 1 0 8 0 65 nfft 64 n 10 M movi ei n 2 n41 for k nin trispect ma ar nfft k 2 n M k n 1 getframe end movie M clear We computed several slices of thetheoretical trispectrum of an ARMA process theslices correspond to f3 0 5 0 45 0 45 0 5 Then we used MATLAB s movi e command to display them movi ei n andgetframe are also MATLAB commands It is instructive to change the MA parameters to 1 1 and to compare the f3 0 slice with the theoretical bispectrum created by bi spect NOTE if you encounter OUT OF MEMORY problems in running the example dear the workspace and then try again If the problem persists decreasen to say 5 or install more memory on your machine Summary We have discussed various algorithms for estimating bispectra and trispectra as well as algorithms for the blind system identification problem Sample estimates of cross cumulants of orders 2 3 and 4 may be obtained via cum2x cum3x andcum4x autocumulants can be estimated via cumest Nonparametric estimates of the bispectrum can be obtained via routines bispeci andbispecd which implement theindirect and direct estimators The direct estimate of the cross bispectrum can be obtained via routinebi specdx Auto and cross bicoheren
171. oisy the TLS solution if one exists minimizes the Frobenius norm of the matrix E r subject to the condition that b r isin therangespace of A E Here E anda have the same dimensions Leta bemxn bbemxk and m2 n k Consider the SVD of the composite matrix Ab UXV where is the diagonal matrix of singular values o t t 1 n k In 1 itis shown that if s n gt s n 1 then the TLS solution exists is unique and is given by 1 x ViVa where the n by n matrix and the k by k matrix V22 are defined by Vii V y 2 Voi V22 Weighted solutions and other details are discussed in 1 1 Golub G H and C F Van Loan Matrix Computations pp 420 5 Baltimore TheJ ohns Hopkins University Press 1983 trench Purpose Syntax Description Algorithm References Estimates the AR parameters and the reflection coefficients given a slice of a cumulant sequence amat cmat pf gamf gamb trench c r Given a nonsymmetric Toeplitz matrix A described by the first column vector c and the first row vector r obtains the solutions to the set of equations AkKk k Pk0 7 k 1 M whereM 1isthelength of vectors c andr and A K sStheleading k x k submatrix of A Ifr is not specified it is set toc the conjugate transpose ofc The typical situation is the estimation of AR parameters based on correlation or cumulant sequenes amat is the matrix of estimated forward AR prediction vectors for
172. order default value is n 10 where n is the length of x a is the vector of estimated amplitudes theta is the vector of estimated initial phases al pha is the vector of estimated damping factors fr isthe vector of estimated frequencies The input time series and the time series corresponding to the estimated parameters are plotted The Signal Processing Toolbox function pr ony is used to fit an ARMA p p 1 model tothetransient signal the MATLAB functionresi due is used to convert the ARMA parameters to the pole residue form these arethen converted tothe amplitude phase frequency and damping factor terms Note that estimates of the amplitude and starting phase are sensitive to the presence of additive noise 1 Krauss T J N Little and L Shure MATLAB Signal Processing Tool box User s Guide The MathWorks Inc 1994 ivcal Purpose Syntax Description Reference Computes instrumental variables Used byri vdl andrivtr Z Z ivcal y ivcal y morder lambda i vcal computes the instrumental variable corresponding to the correlation or to the diagonal slice of third or fourth order cumulants y is the time series If y is a matrix each column is treated independently mor der is the cumulant order it should be less than 5 default value is 3 lambda is the forgetting factor O lt ambda lt 1 the default value is 1 z is the computed instrumental variable z n is obtained from y n as follows
173. orresponding to the ARMA model Synthetic ARMA processes can be generated via routinear masyn routine rpiid can be used to generate i i d sequences with various probability density functions pdf s Polyspectra and Linear Processes MA Models Let us consider the pure MA case that is p O in 1 31 through 1 37 It was established in 19 that the process y n satisfies the set of equations q q Coin Y esb kyC3y n k n k E bikcn k k 0 k 1 1 38 q q 3 Coy n Y eqgb k Cyy n k n Kin Y b 9Cs K k 0 k 1 1 39 wheren q 2q 3 3yb q yoy and 4 Y 4y q yoy Equations 1 38 and 1 39 can be readily verified by using the impulse response summation formulas in 1 16 through 1 18 These equations are sometimes referred to as the GM equations in the literature 19 Equations 1 38 and 1 39 represent a set of 3q 1 linear equations in the 2q 1 unknowns y 44b K k 20 q and b k k 21 q where m 220r m 3 In 19 MA parameters were esti mated from 1 38 or 1 40 by simultaneously solving for the em b k s and b2 k s or b3 k s A disadvantage of the method is that it is overparameterized and does not take into account the relationship between the e 1b k s and the b2 k s or b k s Additive white noise may be permitted provided we eliminate the equations involving C5 0 this eliminates q 1 equations leaving us only 2q equations in the 2q 1 unknow
174. p fq D 1 112 fm n fm 1 n T s n bm n 1 1 113 1 1 bm n bm 1 n 1 p n E f m n fm 1 n 1 114 linear Prediction M odels f b Am 3 0 Aq 1 0 F n Pane ee B n 1 m 1n 9 1 115 b n A n Bn B n 1 a 109 ee Fm 1 1 116 b n 1 b n 1 1 n 1 b _ 4 n 1 Ym n 1 e EET PU ERE m 1 D 1 117 Note that 1 106 through 1 117 arein the proper order required for successful evaluation The initial values for the algorithm are given by AL SOFAS 0s 1 118 Fm 1 0 Bg 0 8 1 1 119 fon bg n y n 1 120 fo n bo n z n 1 121 Fo n Bg n AF g n 1 z m y n 1 122 yo n 1 1 123 The following discussion is adapted from 64 which also gives a derivation of the algorithm Equations 1 111 and 1 112 with theinitial condition in 1 120 definea lattice structure excited by the observed process y n Equations 1 113 and 1 114 with the initial condition in 1 121 define a lattice structure excited by the 1 59 1 Tutorial 1 60 associated process z n The two lattices are coupled through the time update equations for the A s 1 106 and 1 107 and the order update equation for the conversion factor 1 117 When y n z n the two lattices collapse into one In the above f n and b n are the mth order prediction errors for the observed process y n f n and bm n arethemth order prediction errors for the ass
175. pp 984 95 J uly 1989 58 Stoica P T Soderstrom and B Friedlander Optimal instrumental variable estimates of the AR parameters of an ARMA process IEEE Trans Auto Control Vol 30 pp 1066 74 Nov 1985 59 Subba Rao T and M Gabr An Introduction to Bispectral Analysis and Bilinear Time Series Models pp 42 43 New York Springer Verlag 1984 1 Tutorial 60 Swami A System dentification Using Cumulants Ph D Dissertation University of Southern California pp 107 8 1988 61 Swami A Higher Order Wigner Distributions in Proc SPIE 92 Session on Higher Order and Time Varying Spectral Analysis Vol 1770 290 301 San Diego CA J uly 19 24 1992 amp Swami A Third order Wigner distributions Proc ICASSP 91 pp 3081 84 Toronto Canada May 1991 Swami A Some new results in higher order statistics Proc Intl Signal Processing Workshop on Higher Order Statistics Chamrousse F rance pp 135 38 J ul y 1991 64 Swami A Higher order Wigner distributions in Proc SPIE 92 Session on Higher Order and Time Varying Spectral Analysis Vol 1770 pp 290 301 San Diego CA J uly 1992 65 Swami A Pitfalls in polyspectra Proc CASSP 93 Minneapolis Vol IV pp 97 100 Apr 1993 66 Swami A and J M Mendel Adaptive Cumulant Based Estimation of ARMA Parameters Proc Amer Control Conf ACC 88 Atlanta GA pp 2114 19 J une 1988 67 Swami A and J M Mendel L
176. problem it estimates the angular spectrum using the Eigenvector MUSIC Pisarenko ML Capon ML and AR algorithms based on the diagonal slice of the fourth order cumulant Routinedoa implements the corresponding algorithms based on second order statistics it also implements ESPRIT Examples load harm Pxx harmest zmat 12 4 unbiased 256 4 You should see the display in Figure 1 12 sval cum 4 pxx 4 o i EN a ES o 40 o 5 10 15 0 5 0 0 5 o o 2 g E 20 50 40 0 5 o 0 5 0 5 o 0 5 o o 50 LMA 20 m EE a 100 40 0 5 o 0 5 0 5 o 0 5 o o E eo 50 ft 20 zi 100 40 0 5 o 0 5 0 5 o 0 5 Figure 1 12 Singular Values of Cumulant Matrix and Estimated Spectra harmest The data consist of two unit amplitude real harmonics at frequencies 0 1 Hz and 0 2 Hz and are contaminated by colored Gaussian noise with a variance of 0 5 the noise spectrum has a pole at 0 15 Hz with a damping factor of 0 9 The SVD plot is indicative of two real harmonics and the spectrum appears to peak at the correct frequencies 0 1 and 0 2 Hz The fourth order cumulants do not see the Gaussian process at 0 15 Hz Note that the power spectra are displayed in dB scale that is 10log4o9 Pxx Pxx harmest zmat 12 6 unbi ased 256 2 1 75 1 Tutorial You should see the display in Figure 1 13 The SVD plot is indicative of three
177. pute and display the true cumulants The contour plot in Figure 1 1 reveals the basic symmetry of third order cumulants namely C3 11 15 C3y t2 11 Other symmetry properties may be verified by usingcumt rue to estimate the true cumulants of a linear process Qa Leb eS 50 _40 20 10 Oo 10 Figure 1 1 Estimated Third Order Cumulants of an ARMA 2 1 Process cumest Polyspectra and Linear Processes Estimating Polyspectra and Cross polyspectra Estimators of polyspectra are natural extensions of estimators of the power spectrum with some important differences in the smoothing requirements Hence it will be useful to review power spectrum esti mation techniques first Estimating the Power Spectrum Techniques for estimating the conventional power spectrum fall into three broad categories the nonparametric or conventional methods the parametric or model based methods and the criterion based methods Thefirst category includes two classes the direct methods which are based on theFT of the observed data and indirect methods which are based on computing the FT of the estimated autocorrelation sequence of the data The dass of parametric methods includes algorithms such as MA AR and ARMA modeling and eigen space based methods such as MUSIC Min Norm etc which are appropriate for harmonic models Criterion based methods include Burg s Maximum Entropy algorithm and Capon s Maximum Like
178. quadratic phase coupling is discussed further in the section on Nonlinear Processes where we will see that the presence of pronounced peaks in the bispectrum is indicative of nonlinear phenomena Let x n be a zero mean white Gaussian process with variance 62 let y n x2 n It is easy to show that the cross bispectrum of x x y should bea constant equal to 20 This follows by noting that if a b c d are zero mean and jointly Gaussian then E abcd cum a b cum c d cum a d cum b c cum a c cum b d Examples randn seed 0 x randn 64 64 yzx x dbiczbispecdx x x y 128 5 Notice that an apparent structure along the axes which is an artifact due to the removal of the mean consistent estimates along the axes and the 1 19 1 Tutorial 1 20 anti diagonal can be obtained only by sufficient smoothing a smoothing window of size 5 is inadequate The lines o4 0 0 and e 65 0 are called the principal submanifolds 4 6 This example also illustrates that passing a Gaussian process through a nonlinearity makes it non Gaussian Estimating Bicoherence Given estimates of the power spectra and the cross bispectrum we can esti mate the cross bicoherence as indicated in 1 14 It has been shown that consistent estimates of the power spectrum and the bispectrum lead to consistent esti mates of the bicoherence Routines bi coher and bi coherx may be used to estimate the autobicoherence and t
179. r M R and C L Nikias Bispectrum estimation a parametric approach IEEE Trans ASSP Vol 33 pp 1213 30 Oct 1985 2 71 rivdl Purpose Syntax Description Algorithm 2 72 Adaptive AR parameter estimation using the double lattice form of the recursive instrumental variable algorithm arvec fref bref fpe rivdl y arvec fref bref fpe rivdl y morder arorder lambda delta thres nsmuth The recursive instrumental variable double lattice algorithm is applied to estimatethe AR parameters of a possibly nonstationary time series T he model is assumed to be causal y isthetime series must be a vector mor der specifies the cumulant order to be used it should be 2 3 or 4 The default value is 4 arorder istheAR order the number of stages in the lattice The default value is 2 lambda istheforgetting parameter O lt I ambda lt 1 The default valueis 0 998 delta istheinitialization value for F 0 and B 0 Its default value is 0 01 thres is the threshold check for division by zero lf x lt t hr es 1 xis set to zero The default value is 0 0001 1 ns mut h is the window length for estimating the steady state AR parameters Its default value is min length y 4 50 arvec is the steady state AR parameter vector of length arorder 1 The last ns muth samples of the reflection coefficients t ref andbref are averaged and then converted to the AR parameters fref isannsamp xa
180. r k gt p then a useful rule of thumb is to choose p as the number of harmonics With finite data records one should expect a slow decrease in the singular values Usually there will bea significant drop in the singular value from o p to o p 1 Iff Iag is biased biased sample estimates of cumulants are computed default if the first letter is not b unbiased estimates are computed nfft specifies theFFT length its default value is 256 norder Specifies the cumulant order touse it should beeither 2 for covariance based estimates or 4 for estimates based on the diagonal slice of the fourth order cumulant The default value is 4 Pxx isannfft 2 x 7 matrix whose first five columns respectively contain estimates of power spectra obtained via the MUSIC Eigenvector Pisarenko ML Capon and AR methods based on the diagonal slice of the fourth order harmest Algorithm cumulant or the covariance function The sixth column contains the standard power spectrum esti mate obtained via the periodogram method The last column is the estimate obtained by applying the minimum norm algorithm AII estimated spectra are normalized to a maxi mum value of unity for purposes of display only The estimated power spectra are displayed on theMATLAB graphics window ar 1 is the estimated parameter vector for the AR method ar2 is the estimated parameter vector for the minimum norm method Let C denote the max ag by max ag matrix
181. radar is defined as the 2 D FT of the WS and can also be expressed as AF 0 1 fr t t 387 dt 1 172 Time Frequency Distributions Cohen s general class of TFD s is defined by j2nt0 j2 e T j2nf Welt f ec AF 6 2 dt de 1 173 where the kernel 0 8 t specifies the particular distribution that is obtained Let x n y n 4 z n it follows that WyyCf n WyC n WzzCf n Wy f n Wz Cf n 1 174 Sincethe WS is a quadratictransform the WS of the sum of two signals is not the sum of the individual WS s but also has cross terms These cross terms make it difficult to interpret the WS These cross terms may be suppressed by appropriate filtering in the ambiguity function domain that is by appropriate choice of the kernel 0 1 In order to reduce the effects of the cross terms Choi and Williams 10 proposed using the kernel 22 2 0 1 exp 0 Tt 8 7 p 0 t 7o 1 175 where the parameter o controls the amount of attenuation the amplitude of the cross terms is proportional to c Unfortunately increased suppression of cross terms invariably leads to smearing or loss of resolution of the auto terms in the time frequency plane n practice signals aresampled in ti me and we compute F Ts also on a sampled frequency grid Thediscretization in timeand frequency of the continuous ti me WS leadstothe nonaliasing requirement that the original signal be sampled at twice the Nyquist rate
182. re all equal to o With M p 1 it follows from 1 133 that 48 R c yy wa 1 135 Thevector a is the eigenvector corresponding tothe minimum eigenvalue of the p 1 x p 1 correlation matrix Ryy Once a has been estimated the frequencies can be estimated as the roots of the A z polynomial The minimum eigenvalue solution to 1 135 is known as Pisarenko s method and is implemented in routines har mest and doa In order to ensure positi ve definiteness of matrix 3iyy it is essential to use biased sample estimates of the autocorrelation For a signal consisting of p real harmonics the order should be taken as 2p The Pisarenko algorithm separates the eigenvectors and eigenvalues of the matrix Ryy into two classes those associated with the signal and those associated with the noise Extensions of the idea have resulted in the high resolution eigen subspace methods such as the Min Norm and MUSIC algorithms 1 67 1 Tutorial 1 68 Multiple Signal Classification MUSIC Let M gt p The autocorrelation matrix Ry has a linearly independent orthogonal set of eigenvectors Denote the eigenvalues and eigenvectors by A and vy Define the M element signal vector i E H22M f T e E d f qgi2m e Pry 1 136 The eigenvectors corresponding to the p largest eigenvalues are said to constitute the signal subspace the remaining M p eigenvectors constitute the noise subspace The signals that is the harmonics in 1
183. re z n y n lead to inconsistent estimates because of the noise w n If w n can be modeled as MA noise that is Qu w n hy n k g n k 0 where g n is i i d then consistent estimates can be obtained by using z n y n d where d gt p qw Process z n is called an instrumental variable IV if 1 it is uncorrelated with the additive noise w n that is E w n z n m 0 m gt 0 and 2 with m 1 M gt p we get a full rank set of linear equations from 1 01 The delayed process z n y n d is an IV when w n is MA qy noise and d gt p qw Appropriate choices of z n lead to estimates based on higher order cumulants for example z n y n y n p leads to estimates based on the 1 D slice C3y m p of the third order cumulant other examples are discussed in 64 For a discussion of optimal IV s see 58 When z n z y n theresulting matrix is no longer Hermitian symmetric A variation of RLS that is appropriate for this case is called the Recursivel nstrumental Variable RIV algorithm and is implemented inri vtr Thealgorithm is given by linear Prediction M odels Compute the weight vector w i recursively as X 1P n 1 z n koce II 1 A y n P n 1 z n 1 102 a n din w n y n 1 103 w n w n 1 k n a n 1 104 1 1 H P n A P n 1 A k n y n P n 1 1 105 where the desired signal d n y n 1 The initial values are w 0 20 and P 0 281 where 8 is a small positiv
184. rgodic that is consistent estimates of the third order cumulant C3 t1 t2 cannot be obtained from a single realization unless frequency coupling is always accompanied by phase coupling 68 Indeed the random phase assumption itself suggests that multiple realizations are required Further creating multiple realizations by record segmentation does not lead to consistent estimates Given a single realization C3 14 15 will show impulses if frequency coupling exists that is f3 f2 f1 Given multiple realizations C3 t1 t2 and C3 t t will be nonzero only if both frequency and phase coupling exist that is f3 fo f1 and 63 202 64 The parametric bispectrum in 1 168 shows a peak whenever A n4 A n2 and A n n2 have peaks Hence when Ng gt 1 spurious peaks might result from interactions between different triplets of the form f n fin In the case of Nq gt 1 interactions may occur between different QPC triplets leading to additional peaks in the parametric bispectrum given by 1 168 Summary The Higher Order Spectral Analysis Toolbox offers two algorithms to estimate the parameters of a second order Volterra model n t i ck which assumes that theinput is Gaussian and nl pow which is applicable to arbitrary inputs Note that both algorithms require access to inputs and outputs Quadratic phase coupling QPC can be detected using routinegpctor Routineqpcgen can be used to generate synthetics for the QPC prob
185. rial The third order hologram h t is then defined by 39 40 Byyx fr f2 1 194 i j2nD f fas df exo G2xf20e 1 195 8 t D ee 1 196 The hologram should display a strong peak at the location of the true delay Since third order statistics are used the method is insensitive in theory to both spatially and temporally colored noise which is symmetrically distributed e g Gaussian In practice the estimated hologram h t will not be an impulse hence we estimate D as the index t which maximizes h x This algorithm is implemented in routinet deb 1 106 Time Delay Estimation Examples load tdel delay tdeb s1 s2 30 Delay estimated by TDEB is 16 You should see the display on Figure 1 28 and the E sti mated delay E sti mated delay 16 The two signals are corrupted by spatially correlated noise the cross correlation between the noises at the two sensors has a strong peak at a delay of 5 samples the signal delay is 16 samples TDEB Hologram delay 16 70 T T T 60r 4 40 E 30 J 20 J 10 80 80 Figure 1 28 Time Delay Estimated Using Cross Bispectrum tdeb Summary The basic idea in time delay estimation is to locate the peak in the cross correlation between two signals routinet der implements the Maximum likelihood window of Hannan and Thompson and is useful if the noises at the two sensors are not spatially correlated Routinest de andtder u
186. rorder array containing the forward reflection coefficients of the upper lattice where ns amp is the length of the time series y bref isannsamp xarorder array containing the backward reflection coefficients of the upper lattice f pe is the final prediction error Details of the algorithm are given in the Tutorial rivdl See Also ivcal ri vtr References 1 Swami A and J M Mendel Adaptive Cumulant Based Estimation of ARMA Parameters Proc Amer Control Conf ACC 88 Atlanta GA pp 2114 19 J une 1988 2 Porat B B Friedlander and M Morf Square root covariance ladder algorithms in Proc ICASSP 81 Atlanta GA pp 877 880 March 1981 2 73 rivtr Purpose Syntax Description Algorithm See Also 2 74 Adaptive AR parameter estimation using the transversal form of the recursive instrumental variable algorithm arvec fpe wt rivtr y arvec fpe wt rivtr y morder arorder ambda delta nsmuth The transversal form of the recursive instrumental variable algorithm isapplied to estimate the AR parameters of a possibly nonstationary time series The model is assumed to be causal y isthetime series must be a vector mor der specifies the cumulant order to be used it should be 2 3 or 4 The default value is 4 arorder istheAR order the default valueis 2 lambda istheforgetting parameter O lt I ambda lt 1 The default valueis 0 998 delta istheinitialization value for
187. ross the suite of realizations cum3x Algorithm See Also Let x n y n and z n n 1 N denote the three time series Let M denote the number of samples per segment the value of the variablesamp seg Let O denote the percentage overlap between segments the value of the variable overlap Let My M M x O 100 Then thetimeseries x y and z are segmented into K records of M samples each whereK N M O 100 M The kth record or segment of x consists of the samples x i x i kK 1 M1 i 1 M k 1 K y i and z i are defined similarly The sample mean is removed from each record and sample estimates of the third order cross cumulants are obtained as 1 CK 7 nn yxy i mz m i where the summation over i extends from 1 max 0 m n toM max 0 m n The normalizing parameter M m n is identically equal to M for biased estimates and equals M max 0 m n max 0 m n for unbiased estimates The cumulants estimated from the K records are then averaged to obtain the final esti mate C m n Ale K Y C m n k 1 cum2x cum4x cumest 2 37 cum4x Purpose Syntax Description 2 38 Computes the fourth order cross cumulant of four signals cvec cum4x w x y z maxlag samp_seg overlap flag kl k2 Computes the fourth order cross cumulant of the four signals w x y and z which should have identical dimensions maxlag specifies the maximum lag of the cumulant t
188. rs of an ARMA process using cumulants IEEE Trans on Automatic Control Vol 37 pp 268 73 Feb 1992 2 10 armarts Purpose Syntax Description Estimates ARMA parameters using the residual time series method avec bvec armarts y p q avec bvec armarts y p q norder maxlag samp seg overlap flag ar marts estimates the ARMA parameters of the ARM A p q process y via a three step procedure First AR parameters are estimated via the function arrcest next the MA q residual time series y n um oa K y n K is created and finally its MA parameters are estimated via thefunction maest p and q are the AR and MA orders Variable nor der specifies the cumulant order to be used nor der should be 3 or 4 a negative value indicates that autocorrelations as well as cumulants of order nor der should be used to estimate the AR parameters The default valueis 3 maxlag specifies the maximum number of cumulant lags to be used its default value is qtp Variablessamp seg overlap andflag control the manner in which sample cumulants are esti mated samp seg specifies the number of samples per segment the default valueisthe length of thetime series overlap specifies the percentage overlap between segments the allowed range is 0 99 and the default valueis O Ifflag is biased then biased sample estimates are computed default if the first letter is not b unbiased estimates are computed
189. s Algorithm The signal is assumed to be described by p q 90 x n Y a k x n k Y b k u n k h k u n k k 1 k 0 k 0 y n x n g n where u n is i i d non Gaussian and g n is Gaussian The AR parameters are obtained as the least squares solution to the normal equations given by p y a K C4 m k p 0 m gt q k 0 p a k C m k p 0 0 m gt q k 0 where a 0 21 m q 1 maxl ag andp q p q Since one must have atleast m q 1 q p 1 2 it follows from the preceding normal equations that the cumulant lags involved must include qt 1 p q p The default value of max ag follows from this The impulse response is then estimated via p h n 2j o0 Cs a km n 1 q Ve oa K C3 q k 0 or via P a K C4 q k n 0 h n Dic o9 Ca Qq k n 0 n 1 q Veep atk Cala k 0 0 depending upon whether third or fourth order cumulants are used We assume b 0 h 0 1 2 9 armaqs In our implementation we estimate the AR and IR parameters simultaneously and we use the Total Least Squares TL S solution seet s TheMA parameters are then obtained via p b n Y a k h n k 0 n 1 q k 0 See Also arrcest maest cumest armarts tls References 1 Swami A and M Mendel ARMA parameter estimation using only output cumulants IEEE Trans ASSP Vol 38 pp 1257 65 J uly 1990 2 Swami A and J M Mendel I dentifiability of the AR paramete
190. s Table 1 1 Summary Statistics for Sunspot Data Original Differenced Mean u E amp t 48 434 0 085 Variance c E x t u 1548 781 547 773 Skewness E x t u 3yc 1 038 0 904 Kurtosis E x t u yo 3 0 657 1 502 Thedata in Figure 1 29 display an apparent periodicity we can usehar mest to verify this The 1 1 panel of Figure 1 30 displays the singular values of the covariance matrix the plot indicates that only three singular values are significant p 3 hence we used an order of 3 The corresponding power spectral estimates are shown in the remaining panels of Figure 1 30 All of the estimates havea strong peak at 0 Hz and at approximately 0 1 Hz We can use pi ckpeak to accurately locate the peaks in the power spectra we can also use roots to esti mate the locations of the roots of the AR polynomial esti mated by the AR method doing so we obtain an estimate of 10 64 years for the period x 10 sval cum 2 pxx eig a wl Oo music S Qoo pisar I a o o c ml S o o 5 o 0 5 0 5 0 0 5 o 0 D sp x 5 o 5 100 20 amp E 200 40 0 5 0 0 5 0 5 0 0 5 Figure 1 30 Sunspot Data Power Spectra 1 109 1 Tutorial 1 110 We used gl st at totest for Gaussianity and linearity the test results are Test statistic for Gaussianity is 357 4639 with df 60 Pfa 0 Linearity test R estimated 14 8592 lambda 11 0332 R
191. s data file contains the 4096 by 8 matrix y mat It is a synthetic for the DOA problem with two sources at bearings 01 2 15 and 05 25 the eight sensors are spaced half a wavelength apart The source signals are Laplace distributed and have unity variance Spatially correlated Gaussian noise with unity variance was added to obtain the noisy signals The spatial color is such that the angular noise spectrum of the noise shows sharp peaks at 30 gl dat mat This file contains the vectorse g u x andz each vector consists of 512 samples Sequencese g andu are realizations of i i d sequences with exponential Gaussian Laplace and uniform distributions Sequence x is obtained by passing e through the AR filter 1 1 5 0 8 and zin x n harm mt A harmonics in noise synthetic Two unity amplitude harmonics with frequencies 44 0 1 and Az 0 2 were corrupted with additive colored Gaussian noise with a variance of 0 5 Since each of the harmonics has power 0 5 the SNR is 3 dB Thirty two independent realizations each with 128 samples were generated The phases of the harmonics were chosen randomly for each realization The colored Gaussian noise was generated by passing white Gaussian noise through an AR filter with coefficients 1 1 058 0 81 The noise spectrum has a pole at 0 15 Hz with a damping factor of 0 9 The set of realizations is contained in the matrix z mat the columns correspond to independent re
192. s not b unbiased estimates are computed only the first letter of f ag is checked Cumulants are esti mated from each record and then averaged across the set of records h the default valueis 2 p q the estimated impulse response will range from samples I h tol h hest is the estimated impulse response h n n 2 1 h I1 h theimpulse response is normalized so that h 0 1 Samples of h n n lt 0 will have significant values if the original linear system is not minimum phase e g the original system is an ARMA model some of whose zeros or poles lie outside the unit cirde ceps is the estimated complex cepstrum Ri n n q p 2 19 biceps Algorithm See Also Reference 2 20 a is the vector of the minimum phase cepstral parameters AM n 1 p b is the vector of the maximum phase cepstral parameters B n 1 q mi nh is the minimum phase component of the I R maxh is the maximum phase component of the I R so that hest conv minh maxh The cepstral parameters and complex cepstral coefficients are related via n A 7n n hine 0 B n n lt 0 Any linear system can be approximated by an MA L model provided L is large enough If the MA model has Li zeros inside the unit circle and L L Li zeros outsidethe unit circle then the estimated impulse response will show an apparent shift of L osamples to the left of time zero Details of the algorithm are given in the Tutorial
193. scale and shift ambiguity Any linear System can be approximated by an MA L model provided L is large enough If the MA model has Li zeros insidethe unit circle and Lo L Lj zeros outside the unit circle then the estimated impulse response will show an apparent shift of L samples to the left of time zero Details of the algorithm are given in the Tutorial biceps 1 Pan R and C L Nikias The complex cepstrum of higher order cumulants and non minimum phase system identification IEEE Trans ASSP Vol 36 pp 186 205 F eb 1988 bicoher Purpose Syntax Description Algorithm Bicoherence estimation using the direct FFT based method bic waxis bicoher y bic waxis bicoher y nfft wind samp seg overlap The bicoherence of the process y is estimated via the direct FFT based method y isthe data vector or matrix nfft specifies the FFT length touse for computing the bicoherence the nominal default value is 128 if nf ft is smaller than samp seg the power of 2 just larger than samp seg will beused wi nd specifies the time domain window to be used it should be a vector of length seg samp By default thehanni ng window is used Data segments are multiplied by the time domain window then Fourier transformed to compute the frequency domain double and triple products The Fourier transform of the window function should be real and nonnegative samp seg specifies the number of samples per segment or r
194. scellaneous Miscellaneous hel p hosa will display thelist of HOSA Toolbox M files along with one line descriptions For additional help on a function typehel p functi on name For example to get help on function doa typehelp doa Most of the arguments in the Higher Order Sectral Analysis Toolbox M files have default settings Recall that arguments in a function have positional significance Thus if you want to specify the fourth argument of a function you must specify at least the first four arguments Function invocations of the formfunction name a b are invalid Vectors must be specified within square brackets e g 1 2 3 Text or string variables must be specified within single quotes that is biased For more details seethe Tutorial section of the MATLAB User s Guide For details on the algorithms and on the theory of cumulants and polyspectra see the Tutorial section of this manual Prompting Several Higher Order Sectral Analysis Toolbox routines such asar masyn and qpcgen areinteractive the routines ask the user to specify various parameters we refer to this as prompting Guided tour Invoke the function hosademo for a guided tour of the Higher Order Sectral Analysis Toolbox demo will take you through all the examples in the manual Addenda SeetheReadme mfile or typewhatsnewhosa attheMATLAB command prompt armaqs Purpose Syntax Description Estimates ARMA parameters using the q sli
195. sequence the fluctuations in the periodogram become more rapid as N increases As an example try semilogy abs fft rpiid n nor for increasing values of n Proper smoothing smooths out the fluctuations and yields consistent estimates It should be emphasized that the variance expressions are meaningful only in the context of random processes the FT itself is a very useful tool for analyzing deterministic signals The periodogram esti mate can be made consistent in several ways by e Smoothing filtering in the frequency domain Multiplying the autocorrelation sequence by a lag window function Multiplying the time domain data by a window function or Averaging several periodogram estimates These observations carry over to bispectral estimates as we see next Estimating Bispectra and Cross Bispectra The natural estimate of the cross bispectrum is the FT of the third order cumulant sequence that is N 1 N 1 eh Y Y Cxyz k 1 k N 11 N 1 j2nf k 2nfil e fika SRTA N 1 26 Polyspectra and Linear Processes where Xy f is the FT of x n _ 9 This estimate known as the cross bi peri odogram is not a consistent estimate As in the case of the power spectrum the estimate can be made consistent by suitable smoothing The bispectrum and the biperiodograms are special cases obtained when x y z Smoothing can be accomplished by multiplying the third order cumulant estimates by a lag window func
196. sethird order cross cumulants and cross bispectra respectively hence the sensor noises may be spatially correlated provided they are symmetric cross bispectrum of the noises is zero Routinet degen can be used to generate synthetics for the TDE problem 1 107 1 Tutorial Case Studies 1 108 In this section we will use the Higher Order Spectral Analysis Toolbox to process some real data in doing so we will learn some of the pitfalls and tricks of the trade Sunspot Data Schuster was perhaps the first to use the periodogram to analyze sunspot data in the nineteenth century subsequently these data have been virtually adopted as a standard by statisticians and signal analysts These data are known to have an approximate 11 year cycle Sunspot data for the years 1700 1987 are available in the files unspot dat which is included in the standard distribution package Figure 1 29 shows the data and the histogram the latter is indicative of an exponential distribution We used cumest to obtain the estimates of some summary statistics which are shown in the third column of Table 1 1 the last column is discussed later 200 T T 150r 4 100r 4 A 4 0 L L L 1 0 50 100 150 200 250 300 100 T T 807 n 607 n 40H J 20H 4 0 L L L L L 1 he 0 20 40 60 80 100 120 140 160 180 200 Figure 1 29 Sunspot Data C ase Studie
197. signal the Wigner trispectrum can be defined in two different ways here we implement the symmetric Wigner trispectrum 1 Definethe instantaneous fourth order product ra t11 12 73 X t t x t t v1 x t T v2 x t T 13 where t 14 445 13 4 Notice that two of the terms are conjugated The symmetric Wigner trispectrum WT is then given by j2 IST 242 3 3 Win fy fo f3 dxidzadsse u fitit fotat fst ra t t4 1 T3 In this routine we compute the slice f1 f2 f3 f Note that the original signal should be sampled at twice the Nyquist rate in order to avoid aliasing Also note that the frequency axes are scaled by the wig4 factor of 1 2 in this routine weundothe scaling sothat you can find the peaks at the expected frequencies See Also wi g2 wi g2c wi g3 wi g3c wi g4c References 1 Fonollosa J R and C L Nikias Analysis of finite energy signals using higher order moments and spectra based ti me frequency distributions Signal Processing Vol 36 pp 315 28 1994 2 Swami A Higher Order Wigner Distributions in Proc SPIE 92 Session on Higher Order and Time Varying Spectral Analysis Vol 1770 290 301 San Diego CA J uly 19 24 1992 2 95 wig4c Purpose Syntax Description Algorithm 2 96 Computes the Sliced Reduced I nterference Wigner trispectrum wx waxis wig4c x nfft sigma flag Computes the Sliced Reduced I nterference Wigner trispectrum WT 2
198. sion ar 1 0 3 0 1 0 39 0 72 ma 1 nlags 6 rvec cumtrue ma ar 2 nlags a2 c2 p2 gf2 gb2 trench rvec 7 13 rvec 7 1 1 linear Prediction M odels The AR matrix a2 should be 1 0000 1 0000 1 0000 1 0000 1 0000 1 0000 0 0764 0 0700 0 0399 0 0300 0 3000 0 3000 0 0 0834 0 0581 0 1000 0 1000 0 1000 0 0 0 3613 0 3900 0 3900 0 3900 0 0 0 0 7200 0 7200 0 7200 0 0 0 0 0 0000 0 0000 0 0 0 0 0 0 0000 Note that the AR vectors for orders greater than three are all essentially the same The prediction error variance p2 should be 2 4049 2 3881 2 0764 1 0000 1 0000 1 0000 The forward reflection coefficients gf 2 should be 0 0764 0 0834 0 3613 0 7200 0 0000 0 0000 Since the Toeplitz matrix in this example is symmetric 9b2 gb1 and c2 flipud a2 Since the autocorrelation sequence corresponds to an AR 4 model the reflection coefficients are zero for m gt 4 and the pfe s remain constant for m 2 4 Now let us repeat this example using third order cumulants rvec cumtrue ma ar 3 third order cumulants nlags rvec rvec nlags l the C m 0 slice c rvec nlagstl 2 nlags positive lags r rvec nlagstl 1 1 negative lags a4 c4 p4 gf4 gb4 trench c r 1 51 1 Tutorial The forward AR matrix a4 should be 1 0000 1 0000 1 0000 1 0000 1 0000 0 2218 0 2050 0 1834 0 3000 0 3000 0 0 1163 0 0673 0 1000 0 1000 0 0 0 2757 0 3900 0 3900 0 0 0 0 7200 0
199. sity and the University of Southern California respectively He has held positions with Unocal USC CS 3 and Malgudi Systems He is currently a Senior SCEEE Research Fellow at the Army Research Lab Adelphi MD Dr Swami has published over fifty journal and conference papers in the areas of modeling and parameter estimation of non Gaussian processes He is co organizer and co chair of the Eighth IEEE Signal Processing Workshop on Statistical Signal and Array Processing Corfu Greece J une 1996 Jerry M Mendel J erry Mendel received his B S degree in Mechanical Engineering in 1959 his M S in 1960 and his Ph D in 1963 in Electrical Engineering from the Polytechnic I nstitute of Brooklyn NY Currently he is Professor of Electrical Engineering at USC in Los Angeles Dr Mendel is a Fellow of the IEEE a Distinguished Member of the IEEE Control Systems Society a member of Tau Beta Pi Pi Tau Sigma and Sigma Xi and a registered Professional Control Systems Engineer in California He has authored more than 300 technical papers three textbooks and four other books related to his research in estimation theory deconvolution higher order statistics neural networks and fuzzy logic Chrysostomos L Max Nikias Chrysostomos L Max Nikias received his B S degree in Electrical and M echanical Engineering from the National Technical University of Athens Greeceand his M S and Ph D degrees in Electrical Engineering from the State Un
200. so that we may conclude that three of the four harmonics are quadratically phase coupled The singular value plot shows six significant singular values corresponding to one quadratically coupled triplet as in the case of the power spectrum harmonics in noise overestimating the number of harmonics usually leads to better results in this case we used an AR order of 12 Power Spectrum Singular values 0 0 2 0 4 0 6 0 5 10 15 20 frequency Estimated Parametric Bispectrum T T T T T T T 0 2r pi 5 a 0 15F gae ad e Pd O 1rF Er d J 0 05 L 1 4 0 i 1 L L L A L 1 1 1 0 0 05 0 1 0 15 0 2 0 25 0 3 0 35 0 4 0 45 Figure 1 19 Parametric Estimate of bispectrum qpctor Ideally A n contains impulses at f K 1 2 3 and S3x n1m2 is given by 3 3 3 S3x Np 12 x tfo Y notfp Y Sev emat k 1 k 1 m 1 1 169 1 87 1 Tutorial 1 88 From 1 166 we note that it suffices to examine the bispectrum in the region O lt n lt 1 lt 1 2 If fo 2f1 the parametric bispectrum will show a single impulse at f4 f2 corresponding to k 1 22 m 3 in 1 169 indicating the existence of the quadratically phase coupled triplet f4 f2 f4 f2 If fo 2f4 the bispectrum will have an additional peak at n4 n2 f1 f1 corresponding tok 1 21 m 2 in 1 169 which can also result from the phase coupled triplet f f 2f It is important to note that the process in 1 163 is not e
201. ssianity linearity tests har m mat Synthetic for the harmonics in noise problem mal mat MA 3 synthetic nl 1 mat Synthetic for testing the nonlinear routines nl 2 mat Synthetic for testing the nonlinear routines qpc mat Synthetic for the quadratic phase coupling problem riv mat Synthetic for the adaptive LP problem t del mat Synthetic for the time delay estimation problem tprony mat Synthetic for testing Prony s method wigdat mat Synthetic for testing the Wigner routines In the following signal to noise ratio SNR is defined as the ratio of the signal variance to the noise variance when expressed in dB it is given by 1010919 65 6 where o is the variance of the signal and oj isthe variance of the noise All frequencies arein Hz relative to a nominal sampling frequency of 1 Hz arl mat An AR 2 synthetic the AR parameters are 1 1 5 0 8 the time series length is 1024 samples input distribution is single sided exponential additive white Gaussian noise was added so as to obtain a SNR of 20 dB The vector y contains the AR 2 time series 1 134 Data Files armal mat An ARMA 2 1 synthetic the AR parameters are 1 0 8 0 65 the MA parameters are 1 2 thetime series length is 1024 samples input distribution is single sided exponential additive white Gaussian noise was added so as to obtain a SNR of 20 dB The vector y contains the ARMA 2 1 synthetic The model is mixed phase doal mat Thi
202. stimates are computed cvec Will contain the sample estimates of E x n i n m iiy m max ag s s m xl ag Here ux denotes the mean of process x and the superscript denotes complex conjugation If x y are matrices columns are assumed to correspond to different realizations in this case overlap is set toO and samp seg tothe row dimension the cross cumulant is estimated from each realization and then averaged across the suite of realizations cum2x Algorithm See Also Let x n and y n n 21 N denote the two time series Let M denote the number of samples per segment the value of the variablesamp seg Let O denote the percentage overlap between segments the value of the variable overlap Let My M M x O 100 Then the time series x andy are segmented into K records of M samples each where K N M O 100 M The kth record or segment of x consists of the samples X i 2x i kK 1 M1 i 1 M k 1 K y i is defined similarly The sample mean is removed from each record and sample estimates of the cross covariance are computed as 1 a Ck gg LO m where the summation over i extends from 1 max 0 m to M max 0 m The normalizing parameter M m is identically equal to M for biased estimates and equals M max 0 m max 0 m for unbiased estimates The cumulants estimated from the K records are then averaged to obtain the final estimate K C m
203. systems via second and fourth order cumulant matching IEEE Trans on Information Theory Vol 33 3 pp 393 407 May 1987 75 Tugnait J K Approaches to FIR system identification with noisy data using higher order statistics IEEE Trans ASSP Vol 38 pp 1307 17 J uly 1990 76 Zhou G and G B Giannakis Self coupled harmonics stationary and cyclostationary approaches Proc ICASSP pp 153 156 Adelaide Australia Apr 1994 77 Zohar S Toeplitz Matrix Inversion the Algorithm of W F Trench J Assoc Comp Mach Vol 16 pp 592 601 1968 1 Tutorial Reference 2 Reference Function Tables The following tables summarize the suite of M files in the Higher Order Spectral Analysis Toolbox by functionality Reference entries for the functions arethen listed in alphabetical order On linehelpis available via thehel p and hosahel p commands Higher Order Spectrum Estimation Conventional Methods Function Purpose cum2x Estimates second order cross cumulants cum3x Estimates third order cross cumulants cum4x Estimates fourth order cross cumulants cumes t Estimates auto cumulants orders 2 3 4 bi coher Estimates auto bicoherence bi coherx Estimates cross bicoherence bispecd Direct method for bispectrum estimation bispecdx Direct method for cross bispectrum estimation bispeci Indirect method for bispectrum estimation gl stat Detection statistics for Gaussianity and linearity tests 2 2
204. t involve h 0 the cepstrum at the origin is related to the overall gain of the system Because of the inherent scalar ambiguity in estimating H z from power spectra or cumulant spectra we let h 0 1 In practice sample esti mates of the third order cumulants are used this algorithm is implemented in routine bi ceps where k lt max p q and lt max p q 2 which are the recommended ranges in 43 Examples load mal hest ceps biceps y 8 8 128 You should see the display in Figure 1 7 Note that the complex cepstrum decays rapidly to zero The true MA parameters are 1 0 9 0 385 0 771 complex cepstrum T o sample number impulse response T sh 1 1 1 1 1 1 1 40 30 20 10 o 10 20 30 40 sample number Figure 1 7 Cepstrum and IR Estimated by Biceps 1 39 1 Tutorial An alternative approach is based on the FFT We may rewrite 1 58 as 43 7 FT kKC3 k 1 mb m n FT se FT C3 k where IFT denotes 2 D inverse Fourier transform Then we make use of 1 57 to obtain the A k s from mb m n The FFT approach is useful if the cepstrum h k has long support This algorithm is implemented in routine bicepsf More generally the polycepstrum is defined as the inverse FT of the log of the corresponding polyspectrum assuming that it exists b j2 1m i2 fy ama Dkx ms my 4 df dr nf m 6j2n Ins Goss For a linear process the polyc
205. ter than unsegmented estimates In the case of bispectra as we have seen earlier we can obtain consistent estimates by either estimating bispectra from segments and averaging them frequency resolution limited by the segment length or by applying an appropriate smoothing window in the frequency data or lag domain Be aware that the higher order moments and cumulants of complex processes can be defined in different ways and that their polyspectra do not possess all the symmetry properties of their real counterparts C ase Studies There has been lot of confusion about the nonredundant region NRR of the bispectrum we refer the reader to Brillinger 4 for an explanation of why the NRR is the triangle with vertices 0 0 0 2 and 21 3 21 3 and not the triangle with vertices 0 0 0 x and 1 2 1 2 Finally it is easy to be misled by contour plots where100 eps may stand out as a peak in a background of eps careful attention should be paid toscaling the data normalizing the data to unity variance is usually helpful 1 133 1 Tutorial Data Files The Higher Order Spectral Analysis Toolbox distribution diskette includes several mat files Thesefiles are used in the Tutorial to demonstrate various toolbox functions Descriptions of these mat files are given below arl mat AR 2 synthetic armal mat ARMA 2 1 synthetic doal mat Synthetic for the direction of arrival problem gl dat mat Data for Gau
206. th AR parameter vector 1 1 5 0 8 The vector y contains the noiseless signal Additive white Gaussian noise was added to y to obtain the noisy signal zw with a SNR of 10 dB Colored Gaussian noise generated by passing a white Gaussian sequence through the AR filter 1 O 0 49 was added to y to obtain zc also at a SNR of 10 dB Summary Wehave discussed several linear prediction problems TheTrench algorithmis implemented in t rench it computes the backward and forward AR prediction filters associated with a given cumulant slice Routinesrivtr andri vdl implement the transversal and double lattice version of the Recursive Instrumental Variable RIV method these algorithms may be based upon the correlation or a sliceof thethird or fourth order cumulant depending upon the Instrumental Variable IV used i vcal can be used to compute the IV 1 61 1 Tutorial Harmonic Processes and DOA 1 62 In the harmonic retrieval problem the observed data are of the form p y n Y o4 exp j2xn f jO W N x n w n k 1 1 124 where w n is additive noise a s are the amplitudes fy s are the frequencies and 9 s are the phases Additional assumptions are often made regarding the phases the typical assumption is that the phases are independent random variables uniformly distributed over 0 2x1 Such an assumption implies that multiple realizations are available that is p YMN Y apexp j2anty joy mm Wm M
207. than does Capon s maximum li kelihood estimator which in turn is better than the Blackman Tukey estimator Cumulant Based Estimators In 66 it was established that the fourth order cumulants of the process in 1 124 are given by p Cyy l m n Y k 1 1 155 d um ale C4w l m n In particular the diagonal slice is given by p Cyy m m m Y k 1 1 156 DAT TM C4 m m m If w n is Gaussian then Cay m m m 0 If w n is i i d non Gaussian then C4w m l n is a delta function at the origin Note that when w t is Gaussian or i i d non Gaussian 1 156 is similar to the autocorrelation sequence of the noiseless signal Consequently all of theanalyses based on the autocorrelation carry over to the fourth order cumulant In particular the development of the Eigenvector Pisarenko ML Capon AR MUSIC Min Norm and beamformer spectral estimates based on second order statistics can also be based on fourth order statistics 44 Similarly we can base ESPRIT on the fourth order cross cumulant matrix as well The fourth order statistics are most useful when the additive noise is narrow band Gaussian Routinehar mest implements algorithms for the harmonic retrieval problem it estimates the spectrum using the Eigenvector MUSIC Pisarenko ML and AR algorithms based on the diagonal slice of the fourth order cumulant Harmonic Processes and DO A Routinedoa implements algorithms for the direction of arrival
208. the estimates of b k and b k If all the 2 59 maest References 2 60 estimated b k s havethe samesign as the corresponding b3 K s then the final MA parameter estimate is obtained as B k sign by k by k b3 k 9 2 otherwise b k b k 1 Giannakis G B and J M Mendel Identification of non minimum phase systems using higher order statistics IEEE Trans ASSP Vol 37 pp 360 77 Mar 1989 2 Tugnait J K Approaches to FIR system identification with noisy data using higher order statistics IEEE Trans ASSP Vol 38 pp 1307 17 J uly 1990 maorder Purpose Syntax Description Algorithm Estimates the order of an MA process q maorder y qmin qmax pfa flag Estimates the order of an MA process using third order cumulants y is the observed MA process and must be a vector qmi n specifies the minimum expected MA order the default valueis O qmax specifies the maximum expected MA order the default value is 10 pf a specifies the probability of false alarm for the hypthesis testi ng procedure the default value is 0 05 flag ifflag is nonzero the sample estimate of C3 q k its estimated variance the threshold corresponding to the specified pf a and whether the absolute value of the esti mated C3 q k is less than the threshold are displayed on the command window The default value is 1 q is the estimated MA order The basic idea is that for an M A q process
209. ti on based estimates will be biased estimates based on third order fourth order cumulants will be unbiased if the additive noise is symmetric Gaussian ARMA Models As discussed in the previous section we can determine the AR parameters easily We will ignore the estimation errors and assume that k a k k 1 p this is justified if the data lengths are long enough to ensure good estimates of the cumulants In practice errors in estimating the a k s will show up as an additive non Gaussian noise term on the right hand side of 1 46 Consider the residual time series obtained via p zin Y a k y n k Keg 1 45 1 32 Polyspectra and Linear Processes q p Y b k u n k Y a k w n k k 0 k 0 q p Y b k u n k Y w4 n k 0 k 0 1 46 Routinear marts uses the residual time series method to estimate the ARMA parameters it estimates the AR parameters first via the normal equations it then computes the AR compensated time series via 1 45 and finally estimates the MA parameters via routine maest See 1 46 Since w n was assumed to be Gaussian w4 n is also Gaussian if w n is white w n is MA p noise Routinemaest assumes that the additive noiseis white hencethe results ofar marts are meaningful only at high SNR An alternative solution based on is implemented in routine ar mags whichis a q slice method The AR parameters are estimated via the normal equations as before The
210. timate the interquartile range is also printed note that the estimate cannot bereliable if the number of samples is small Let B3 0 02 denote the bispectrum and let P c denote the power spectrum The normalized bispectrum or bicoherence is defined as B3 1 02 bic 0 05 a P 93 yy 2 P yy y 03 7 Under the Gaussianity zero skewness assumption the expected value of the bicoherence is zero that is E 15ic 01 05 0 Thetest of Gaussianity is based on the mean bicoherence power S y bic o 65 where the summation is performed over the nonredundant region of the bispectrum details are given in 1 The statistics is x distributed with p degrees of freedom wherep is a function of the FFT length nf ft and the resolution parameter c par m 1 2 Under the linearity assumption B 2 is constant for all and c Let 1 IN 1 2 cparm X 04 05 B 01 02 gistat An estimate of 2N X Par y is obtained Note that X w gt is chi squared distributed with noncentrality parameter X The sample interquartilerangeR of theX m n sis estimated and should be compared with the theoretical interquartile range of a chi squared distribution with two degrees of freedom and noncentrality parameter i These statistics are defined in Hinich s paper 1 Sankaran s approxi mations in 2 are used to estimate the Pf a and the theoretical interquartile value of X20 References 1
211. ting the Power Spectrum ssssuss 1 15 Estimating Bispectra and Cross Bispectra 1 16 Examples vestre ec cee EVER YSTRUUE PRIX 1 18 EXAMP ES 4 occuro ESI E GENRE uu erate d 1 19 Estimating Bicoherence 00 cece eee 1 20 Examples ud x Seg Ravens See RENEW Fue RUE Re Aces 1 20 Testing for Linearity and Gaussianity 1 22 Examples nea her dlowew Dog ieee O derer 1 24 Contents Parametric Estimators ARMA Models 1 26 MA Models rela Re ue DADO Ee as 1 29 Examples i ede hed ve wed saeees YER DA ER wee ee YR 1 30 AR Models ranerne arenie niara kbi n ba eae toate 1 31 Examples si iminent e e 745 e SG Wr ature yee s 1 32 ARMA Models sese 1 32 Examples ice recwekkR EG eR ETE eed Ped teed 1 34 AR Order Determination 0 00 c eee eee eens 1 34 Examples 2 zu aces Sates his AREE i KEARE RAE 1 35 MA Order Determination 0 00 cece eee ee 1 36 Examples i2 za Re E Mt dla ee a a A a 1 37 Linear Processes Impulse Response Estimation 1 37 ThePolycepstral Methods 000 cee eee eee 1 38 Examples cy niia ia e a ey Mee We Et Rede 1 39 Examples 1 451 Aleta tle a RE ARA RT eet s 1 41 The Matsuoka Ulrych Algorithm suslsss 1 41 Examples esser bw EE Rated Went E pote 1 42 Linear Processes Theoretical Cumulants and Polyspectra 1 43 Examples ste ue eELE uT Y EID A EI Weis 1 43 SUMMARY vs szcdaendb
212. tion Let w t s be a 2 D window function whose 2 D FT is bounded and nonnegative further assume w 0 0 1 f w t sdt ds o EWG foaru drj riw fdf dr 0 The window function w t s must also satisfy the symmetry properties of third order cumulants For example 2 D lag windows may bederived from 1 D lag windows as follows w t s 2 w t w s wtt s which satisfies the symmetry conditions of C3 m n Consider the scaled parameter window wy t s 2 w t M s M and the smoothed esti mate N 1 N 1 Gee fo Y Y Cxyz k 1 Wyy K 1 k N 11 N 1 1 27 _j2nf k 2nfil e fi e Under the assumption that the cross bispectrum Sxyz f1 f2 is sufficiently smooth the smoothed estimate is known to be consistent with variance given by 2 var Sxyel Fy 2 SESS 4 Fa Say F Sox Fa WE dt ds 1 28 for 0 f1 f n Note the implied consistency condition is M e and M2 N gt e as N gt and J w t s dt ds lt The estimator in 1 27 for x y z is implemented in routinebi speci An alternative approach is to perform the smoothing in the frequency domain As in the case of power spectra we may segment the data into K records of 1 17 1 Tutorial 1 18 length L N K compute and average the biperiodograms and then perform the frequency smoothing using the frequency domain filter Ww f 1 f2 the FT of wm t s In this case 2 var Sxyz f 4 f2 Moon f3523 052542
213. tion and then averaged The estimated MA parameters are returned in theg 1 element column vector bvec The GM method obtains the least squares solution to the set of equations q q y eb loCs n k n 19 Y b gR ko Roy k 0 k 1 where n q 2q andnorder 3 Fornorder 4 the equations are maest q q y sb oC4 n k n kn I9 Y b oR o Ro k 0 k 1 wheren q 2q Here R C3 and C4 are the second third and fourth order cumulants respectively In order to handle additive white noise equations involving R 0 are eliminated sincethe variance of the noise is unknown The modified GM method incorporates a fix dueto Tugnait 2 which appends the following set of equations to the preceding set of equations for q e3R n b k C3 k n q C3 n q k 1 n q Q or q e R n L2 b k C4 k n q 0 C n q 0 k 1 depending upon the specified cumulant order The least squares solution tothe combined system of equations is obtained Here 4 y4 b q o and 4 y4b 9 o When third order cumulants are used the method estimates both b k and b k Let DIU and b k denote the estimates of b k and b k If all the estimated b k s are nonnegative thefinal MA parameter estimateis obtained as ck sign b k 0 5 b4 k b gt k otherwise b k b k When fourth order cumulants are used the method estimates both b k and b gt k Let b k and b3 k denote
214. tion of thetrue delay Sincethird order statistics are used the method is insensitive in theory to both spatially and temporally colored Gaussian noise Sample estimates of the auto and cross bispectra are obtained via routine bi specdx thedatais segmented each segment is Fourier transformed and the triple frequency product is computed these products are then averaged across the suite of segments to give the final cross and auto bispectral estimates In order to obtain good estimates it is critical that the data length be much larger than the specified maximum delay The estimated delay will range from nfft 2tonfft 2 1 tde tder 1 Nikias C L and R Pan Time delay estimation in unknown Gaussian spatially correlated noise IEEE Trans ASSP Vol 36 pp 1706 14 Nov 1988 tdegen Purpose Syntax Description Generates synthetic sequences for the time delay estimation TDE probl em 51 52 tdegen s1 s2 tdegen default tdegen generates the data sequences s n 2x n g1 n s2 n Ax n D gz n where x n is a zero mean i i d random sequence with single sided exponential density function g3 n and ga5 n are zero mean colored noise sequences A is the amplitude gain from sensor 1 to sensor 2 and D is the delay of the signal at the second sensor with respec to that at the first sensor If the function is invoked without any input arguments then you are prompted for all parameters otherwis
215. tion between the two signals x t and y t is given by Ryy t AR t D Rw w T d 1 189 where Ry 1 is the cross correlation between the two noise processes and Rss T is the autocorrelation of the signal If the noises are uncorrelated Ryy x will havea peak at x D the unknown delay In practice due to effects of finite length estimates and due to the presence of noise the cross correlation estimate may not have a sharp peak The data may be prefiltered in order to sharpen the peak equivalently we can multiply the estimated cross correlation by a window function Different choices of the window function lead to different estimates The most popular 1 101 1 Tutorial 1 102 window function is the maximum likelihood window of Hannan and Thompson which is described below Let S y f denote the cross spectrum between the two signals x and y and let Sxx f and Syy f denote the autospectra of x and y The squared coherence function is defined by Iscr C eat se x 7 SiS The optimal maximum likelihood window is then 1 Cy GP Wes E Hee DICES Ge COD and the windowed cross correlati on Ry m is the lFT of W f Sxy f Estimates of the auto and cross spectra and the coherence can be obtained via the MATLAB routines pect rum the segment length must be at least twice the expected maximum delay Since good estimates of the spectra demand a large number of segments it is critical that the lengths of
216. tion in time of the higher order moment spectra of the signal Wigner Spectrum The Wigner distribution WD was introduced in 1932 by Wigner in the context of quantum mechanics its usefulness to problems in communication theory was discovered by Ville in 1948 consequently it is often called the Wigner Ville distribution A series of papers by Classen and M ecklenbra ker 11 discussed the usefulness of the WD for time frequency analysis of continuous and discrete ti me signals and was devoted to applications in digital signal processing sampling issues are also discussed in 12 Tutorial overviews of the WD and its relationships with other time frequency distributions TF D may be found in 3 13 25 In order to differentiate the second order WD from the higher order WD s to be introduced later we will refer to the conventional WD as the Wigner spectrum WS In contrast tothe STFT which is resolution limited either in timeor in frequency dictated by the window function and suffers from smearing and sidelobe leakage the WS offers excellent resolution in both the frequency and time domains The Wigner cross spectrum WCS of two signals x t y t is defined via Wyy t o i xt v2yst v2e dt 1 170 ix E 2 Y E 2 e d 1 171 where X o and x t constitute an FT pair The auto WS is obtained when x t y t and is a bilinear function of the signal x t The ambiguity function which is widely used in
217. tween segments the value of thevariableover ap Let Mq M M x O 100 Then thetimeseries w x y and z aresegmented into K records of M samples each where K N M x O 100 M The kth record or segment of x consists of the samples x i x i kK 1 M1 i 1 M k 1 K Wei i and Z i are defined similarly The sample mean is removed from each record and sample estimates of the fourth order cross cumulants are obtained as C m n t Wang rt OR m y i nz t l yw 1 T E mem EM Odd e m Tin Etat et m 1 1 uu sis Event em a sos t 1 Ok mig Ern o a EON m where the summation over i is such that the indices of w x y and z are all in therange 1 M The normalizing parameter M m n t is identically equal to M for biased estimates and equals M max 0 m n t max 0 m n t for unbiased estimates and the normalizing parameter M m is identically equal to M for biased estimates and equals M max 0 m max 0 m for unbiased estimates The cumulants esti mated from the K records arethen averaged to obtain the final estimate K C m n t i Y C m n t k zd cum2x cum3x cumest 2 39 cumest Purpose Syntax Description 2 40 Computes sample estimates of cumulants using the overlapped segment method cvec cumest y cvec cumest y norder maxlag samp_seg overlap flag kl k2 cumest computes sample estimates of a 1 D slice of the cum
218. u SABI A eR A es Sets WE 1 45 Linear Prediction Models L uuuusss 1 47 Levinson Recursion s s s esaesa naear naar 1 47 Trench Recursion ccc eee 1 49 Examples neieecG eae ds Fak Pa eee EASES eee eed 1 50 Deterministic Formulation of FBLS 00 eee 1 53 Adaptive Linear Prediction 00c cee eee eee 1 54 RIV Algorithm Transversal Form lseuss 1 56 Examples sein vic ate te EF eed a EE SYETWYevE ew ua 1 57 RIV Algorithm Double LatticeForm 1 58 Examples ih nag I RIETI ARRA EM A Ie a 1 60 SUMMARY suere aut psa det ets peas aes dvd eda E Dort 1 61 Harmonic Processes and DOA 1 62 Resolution and Variance 0 cee eee 1 65 AR and ARMA Models 00 c eens 1 66 Pisarenkos Method 0 00 cece eee eee eee nee 1 67 Multiple Signal Classification MUSIC 000 1 68 Minimum Norm Method 0 0c cece een eee eee 1 69 ESPRIT cuz yw zi RR e E Y REERIBVWEGSTETS 1 70 Contents Criterion Based Estimators seen 1 72 Cumulant Based Estimators 00 00a eee eee 1 74 EXAMPpPICS ws anevada caw ik IRE eee eee ee weet ee 1 75 Examples acn 21 arot HD IQ ewe pete ee Ee Me 1 77 S mrmatry ezstvrev e ALET DRA WI EZERPSTRIICUSNPESTQT 1 79 Nonlinear Processes usulsslssls sisse 1 80 Solution Using Cross Bispectra lues 1 80 Examples 22 dea REEL Al A RETE a 1 82 Solution U
219. ulants of the process y y is the data matrix or vector norder specifies the cumulant order and should be 2 3 or 4 the default value is 2 maxl ag specifies the maximum lag of the cumulant to be computed its default value is 0 samp seg specifies the number of samples per segment the default value is the length of the time series overlap specifies the percentage overlap between segments maximum allowed value is 99 default valueis O Iff I ag is biased then biased sample estimates are computed default if the first letter is not b unbiased estimates are computed If y isa matrix the columns are assumed to correspond to independent realizations in this caseoverlap isset to zero andsamp seg is set to the row dimension Cumulants are estimated from each realization and then averaged Parametersk1 and k2 control which 1 D slice of the cumulant function is computed their default values are zero cvec Will contain C gt m C3 m k1 or C4 m k1 k2 m maxlag maxl ag depending upon the specified cumulant order Notethat cumest estimates a 1 D slice of the cumulant cumest Algorithm See Also Reference Le y n n 21 N denote the three time series Let M denote the number of samples per segment the value of the variablesamp seg Let O denotethe percentage overlap between segments variableover ap Let M4 M M O 100 Then thetime series y is segmented into K records of M sampl
220. und 0 1 0 1 and symmetric locations and is indicative of possible quadratic frequency coupling The data appear to be band pass in nature Bispectrum estimated via the indirect method 0 44 j i 4 0 17 f2 eo T Figure 1 37 Lynx Data bispectrum As with the sun spot data we repeated the exercise with the first differences corresponding plots are shown in Figure 1 38 through Figure 1 40 Noticethat the histogram appears to be somewhat symmetrical but the skewness is significantly nonzero the correlation matrix indicates two harmonics order 4 as evidenced by the singular value plot shown in the 1 1 panel of Figure 1 38 the estimated cycles are 4 8 years and 9 6 years The bispectrum confirms the coupling between these two periods C ase Studies In this example we used har mest and bi speci to estimate the periods of the harmonics and to detect quadratic coupling differencing the data helped to darify the esti mates of the power spectrum as well asthe bispectrum Our tests confirm that the data are non Gaussian and show evidence of a fundamental period of around 9 6 years as well as a harmonic at around 4 8 years this may be indicative of quadratic nonlinearities 4000 T T 2000 J eo T i 2000r J 4000 L i i 1 1 0 20 40 60 80 100 120 50 T 40 J 30F J 20r J 10 0 L L
221. us assume that we have estimated S and are confident that the data is non Gaussian Now if the data are also linear we expect the squared bicoherence to be constant for all f4 and fp In practice the estimated bicoherence will not be flat we can obtain an estimate of the constant value by computing the mean value of the bicoherence over the points in the nonredundant region let A denote this mean value The squared bicoherenceis chi squared distributed with two degrees of freedom and noncentrality parameter Thesampleinterquartile range R of the squared bicoherence can be estimated and compared with the theoretical interquartile range of a chi squared distribution with two degrees of freedom and noncentrality parameter i If the estimated interquartile range is much larger or much smaller than the theoretical value then we should reject the linearity hypothesis This test is implemented in routine gl st at Note that a zero bispectrum is not proof of Gaussianity since the higher order cumulants and polyspectra need not be identically zero The bispectrum of a time reversible process is identically equal to zero For example consider the linear process x n pe k where u n is an i i d process If u n is symmetrically distributed then the bispectrum of x n is zero If u n is Gaussian the bispectrum as well as all higher order polyspectra of x n are identically zero If u n is Laplace distributed the bispectrum and all odd order
222. w t s dt ds BIS 1 29 for 0 f4 f n Windowing is not required in this case provided K is large however this does not ensure that the bias will go to zero Rao and Gabr 56 Sec 2 4 have derived a bispectral window which is optimum in terms of bias and variance windows satisfying other optimality criteria are discussed in 56 40 The Rao Gabr window is available as an option in routinebi speci Routinebi specdx computes estimates of the cross bispectrum Examples load qpc bspec bispeci zmat 21 64 0 unbiased 128 1 dbspec bispecdx zmat zmat zmat 128 3 64 0 The contour plots of the two estimates of the bispectrum are shown in Figure 1 2 and Figure 1 3 Bispectrum estimated via the indirect method T T T T T T O 0 0 Figure 1 2 Indirect Estimate of the Bispectrum bispeci Polyspectra and Linear Processes Cross Bispectrum T T T O 4r 5 0 3 0 2 E 0 1 4 S 0r J 0 1 F I 0 2 0 3 f 4 0 4 J 0 5 i i i i i i ri f 0 5 0 4 0 3 0 2 0 1 o 0 1 0 2 0 3 0 4 Figure 1 3 Direct Estimate of the Bispectrum bispecdx Both the direct and the indirect estimates reveal peaks at 0 10 0 15 and the 11 other symmetric locations as indicated by 1 11 The data used in this example consist of quadratically phase coupled harmonics with frequencies at 0 10 0 15 and 0 25 Hz and an uncoupled harmonic at 0 40 Hz
223. with respect to its second order statistics that is given Ryy we can always find h k and an uncorrelated process u n such that Ryy m Ryy m where X N ye haoum k In other words the autocorrelation sequence cannot give any evidence of nonlinearity In contrast higher order cumulants can give evidence of nonlinearity Processes of the form x t a t exp j Y E gait whose phase is a polynomial in timet are called polynomial phase processes the FTs of such processes tend to be flat whereas suitably defined slices of higher order spectra reveal structure that permits esti mation of p and the a s To summarize cumulants are useful 1 if the additive noise is Gaussian and the signal is non Gaussian 2 the linear system is non minimum phase that is mixed phase or 3 the process is nonlinear Bias and Variance of an Estimator In practice we estimate cumulants and polyspectra from data These esti mates are themselves random and are characterized by their bias and variance Let x n denote a stationary process we assume that all relevant statistics exist and havefinite values Let s denote somestatistic defined on x n Let Su denote an estimate of the statistic based on N observations x N 9 Since x n is a random process the estimate y is also random clearly will not 1 11 1 Tutorial 1 12 equal s The estimate Sy is a good estimate if it is near s This notion is clarified by
224. y exp j2xAdsin 0 A matrix A has n k entry exp jo k n am 2 o4 m ow m is the source signal vector at the mth snapshot and wm and vm are the noise vectors The bearing information can be recovered from matrix d The auto and cross correlation matrices of y n and z n are given by H H 2 yy E vs ASA ol 1 147 H H H 2 Ryz i E vs AS A 6 1 148 where we have assumed that the additive noiseis white spatially uncorrelated from sensor to sensor and has variance o2 Matrix S Elaman is the signal correlation matrix and is nonsingular if the sources are incoherent The structure of matrix J depends upon the common elements of the two sub arrays If the two subarrays have no common elements then J is the zero matrix In the uniform linear array example considered earlier Zm N 2ym A n m 1 L A and J _ O xL On XA is l AO AxA In any case given the knowledge of the geometries of the two arrays we can compute matrix J 1 71 1 Tutorial 1 72 We can estimate c as the smallest eigenvalue of matrix Ry We can now eliminate the noise contribution and create the signal correlation matrices 2 H Cy 8 o l ASA LEE iat 1 149 2 HH Co R 0 l AS A Meus AMES 1 150 Given two n x n matrices A and B the generalized eigenvalues of A B are defined by 21 MA B z det A zB 0 If B is nonsingular there are exactly n generalized eigenva
225. ype density function exp Single sided exponential lap Double sided exponential Laplacian nor Normal Gaussian bga Bernoulli Gaussian uni Uniform The default distribution is nor p spi ke specifies theevent probability for the Bernoulli Gaussian distribution the default value is 0 1 The output sequence u is normalized to zero mean by subtracting the theoretical not the sample mean Note that the mean of any vector u can be set to a and its standard deviation to b via u a b x u mean u std u tde Purpose Syntax Description Time delay estimation TDE using the parametricthird order cross cumulant method delay avec tde x y max delay delay avec tde x y max delay samp seg svdflag The time delay between two signals possibly corrupted with spatially correlated colored Gaussian noise is estimated using the parametric third order cross cumulant method x andy arethesignals at thetwo sensors they must haveidentical dimensions if they are matrices it is assumed that the columns correspond to different realizations max delay istheabsolute value of the maximum expected delay samp seg specifies the number of samples per segment for estimating of cumulants it defaults to the length of the time series row dimension if x isa matrix svdflag is an optional parameter with a default value of zero if its valueis not zero the SVD of the cumulant matrix see the Algorit

Higher-Order Spectral Analysis Toolbox User's Guide

Contents

Download Pdf Manuals

Related Search

Related Contents