demqpc.m

高阶谱分析工具箱

%DEMQPC	HOSA Toolbox Demo:  Detection of Quadratic Phase Coupling 

echo off

% demos of qpctor (qpcgen)

% A. Swami April 15, 1993
% Copyright (c) 1991-2001 by United Signals & Systems, Inc. 
%       $Revision: 1.6 $

%     RESTRICTED RIGHTS LEGEND
% Use, duplication, or disclosure by the Government is subject to
% restrictions as set forth in subparagraph (c) (1) (ii) of the 
% Rights in Technical Data and Computer Software clause of DFARS
% 252.227-7013. 
% Manufacturer: United Signals & Systems, Inc., P.O. Box 2374, 
% Culver City, California 90231. 
%
%  This material may be reproduced by or for the U.S. Government pursuant 
%  to the copyright license under the clause at DFARS 252.227-7013. 

clear, clc, 
echo on
%       Non-Linear Processes: Quadratic Phase Coupling 
%
% When a signal consisting of two harmonics, say at frequencies f1 and f2,
% with starting phases, ph1(1) and phi(2), is passed through a square-law 
% device (a quadratic non-linearity), the output consists of harmonics at 
% frequencies f1, f2, f1+f2 and f1-f2, with starting phases phi(1), phi(2),
% phi(1)+phi(2) and phi(1)-phi(2).   These harmonics are said to be 
% quadratically phase-coupled.    This phenomenon is observed in several
% situations, such as plasma physics, oceanography, EEG signals, otoacoustic 
% emissions, signals generated by rotating machinery, etc. 

% The bispectrum estimation routines BISPECD and BISPECI may be used
% to obtain (non-parametric) estimates of the bispectrum.  If QPC is present,
% the bispectrum will exhibit sharp peaks.  (Choose the option 
% `Estimating Bispectra ...' from the main menu, to learn more about this.)
% Here we will look at parametric methods, which offer the promise of 
% high resolution.   

% Hit any key to continue
pause
clc 

% The HOSA Toolbox offers two routines:
% 1. QPCGEN may be used to generate synthetics consisting of 
%    sinusoids, whose phases may be independent or quadratically
%    phase coupled. 
% 
% 2. QPCTOR may be used to detect the presence of quadratic 
%    phase-coupling;  this function uses the "Third-Order Recursion"
%    method. 
% 
% We will test QPCTOR on some data generated by QPCGEN. 
% The data for this demo is a synthetic for the quadratic-phase coupling (QPC)
% problem.  The data consist of four unity amplitude harmonics with 
% frequencies, 0.1 Hz, 0.15 Hz, 0.25 Hz and 0.40 Hz.  64 independent 
% realizations, each consisting of 64 samples, were generated using routine 
% QPCGEN. The starting phases for the harmonics with frequencies, 0.1 Hz,
% 0.25 Hz and 0.40 Hz were chosen independently from an uniform distribution; 
% the starting phase of the third harmonic, at 0.25 Hz, was set equal to the 
% the sum of the starting phases of the harmonics at 0.1 Hz and 0.15 Hz. 
% White Gaussian noise with a variance of 1.5 was added to the signal. 
% Note that the harmonics at (0.1, 0.15, 0.25) are both frequency- and phase-
% coupled; whereas those at (0.15, 0.25, 0.40) are frequency-coupled, but
% not phase-coupled.  

% Hit any key to continue
pause
clc

load qpc

% We will use QPCTOR, with maxlag=18, and an assumed order of 12. 
% QPCTOR sets up a set of linear equations based on third-order cumulants.
% If the data has n quadratically phase-coupled frequency triplets,
% then, the rank of this matrix should be 6n. 
% We compute the SVD of this matrix, in order to determine its rank. 
% We also estimate the conventional power spectrum to get a rough idea of
% the number of harmonics present in the data. 
%
% Notice that the power spectrum shows peaks at roughly 0.10, 0.15, 0.25
% and 0.40 Hz.  We note that quadratic phase coupling may exist in the
% triplets (0.1,0.15,0.25) and (0.15,0.25,0.40). 
% You will notice a sharp drop in the singular value plot, after n=6
% Overfitting usually improves the estimate: so we will use an order of 12 

% Hit any key to continue
pause

clf
[arvec, Bspec] = qpctor(zmat,18,12);
set(gcf, 'Name', 'HOSA QPCTOR')

% QPCTOR estimates the conventional power spectrum, as well as the parametric 
%        bispectrum. 
%
% In the plots, you should note that: 
% The singular values drop off sharply after n=6 indicating the presence of
%     three quadratically phase-coupled harmonics. 
% The power spectrum has four peaks corresponding to the four harmonics,
%    (0.1, 0.15, 0.25 and 0.40 Hz). 
% The bispectrum shows a single sharp peak at around (f1,f2) = (0.10,0.15),
%   indicating qpc for frequencies (0.10, 0.15, 0.25). 
% Notice that only the 0 <= f1 <= f2 portion of the bispectrum is displayed.


% Hit any key to return to the main menu
pause 
echo off
clc