Mercurial > hg > ltpda

annotate m-toolbox/classes/+utils/@math/music.m @ 19:69e3d49b4b0c database-connection-manager

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update ltpda_uo.fromRepository
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Mon, 05 Dec 2011 16:20:06 +0100
parents f0afece42f48
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
1 function [music_data,msg] = music(x,p,varargin)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
2 %MUSIC Implements the heart of the MUSIC algorithm of line spectra estimation.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
3 % MUSIC is called by both PMUSIC and ROOTMUSIC.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
4 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
5 % Inputs:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
6 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
7 % x - vector or matrix. If vector it is a signal, if matrix it may be either a data
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
8 % matrix such that x'*x=R, or a correlation matrix R.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
9 % p - scalar or two element vector. If scalar, it indicates the dimension of the
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
10 % signal subspace. If vector, p(2) is a threshold used to determine the
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
11 % aforementioned dimension.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
12 % nfft - (optional) to be used only with PMUSIC. A scalar indicating the number of
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
13 % points used in the evaluation of the pseudospectrum.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
14 % Fs - (optional) a scalar specifying the sampling frequency. If omitted, we work
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
15 % in rad/sample; if empty it defaults to 1 Hz.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
16 % nw - (optional) a scalar or vector indicating either the order of the correlation
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
17 % matrix or (when a vector) a window whose length is the order of the matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
18 % and whose values are used to window each column of the data matrix.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
19 % noverlap - (optional) a integer indicating the number of samples to overlap from
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
20 % column to column.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
21 % strings - Optional input strings are: 'corr', 'EV' and range ('half' or 'whole').
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
22 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
23 % Outputs:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
24 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
25 % msg - a possible error message.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
26 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
27 % music_data - a structure with the following fields:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
28 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
29 % noise_eigenvects - a matrix whose columns are the noise subspace eigenvectors.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
30 % signal_eigenvects - a matrix whose columns are the signal subspace eigenvectors.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
31 % eigenvals - the eigenvalues of the correlation matrix.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
32 % p_eff - the effective dimension of the signal subspace.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
33 % nfft - number of points used to evaluate the pseudospectrum (only used in PMUSIC).
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
34 % Fs - sampling freq.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
35 % range - string indicating whether 'half' or the 'whole' pseudospectrum should be
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
36 % computed. (Only used in PMUSIC.)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
37 % EVFlag - flag, 0 = MUSIC method; 1 = EigenVector method.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
38
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
39 % Author(s): R. Losada
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
40 % Copyright 1988-2006 The MathWorks, Inc.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
41 % $Revision: 1.1 $ $Date: 2010/02/18 11:16:20 $
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
42
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
43 % References:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
44 % [1] Petre Stoica and Randolph Moses, Introduction To Spectral
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
45 % Analysis, Prentice-Hall, 1997, pg. 15
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
46 % [2] S. J. Orfanidis, Optimum Signal Processing. An Introduction.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
47 % 2nd Ed., Macmillan, 1988.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
48
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
49
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
50 xIsReal = isreal(x);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
51 msg = '';
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
52 music_data = [];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
53
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
54 if isempty(p),
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
55 msg = 'The signal subspace dimension cannot be empty.';
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
56 return
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
57 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
58
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
59 [opts,msg] = music_options(xIsReal,p,varargin{:});
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
60 if ~isempty(msg),
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
61 return
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
62 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
63
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
64 % Compute the eigenvalues and eigenvectors of the correlation matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
65 [eigenvals,eigenvects] = computeeig(x,opts.CorrFlag,opts.CorrMatrOrd,opts.nw,opts.noverlap,opts.window,opts.EVFlag);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
66
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
67 % Determine the effective dimension of the signal subspace
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
68 p_eff = determine_signal_space(p,eigenvals);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
69
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
70 % Separate the signal and noise eigenvectors
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
71 signal_eigenvects = eigenvects(:,1:p_eff);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
72 noise_eigenvects = eigenvects(:,p_eff+1:end);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
73
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
74 % Generate the output structure
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
75 music_data.noise_eigenvects = noise_eigenvects;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
76 music_data.signal_eigenvects = signal_eigenvects;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
77 music_data.eigenvals = eigenvals;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
78 music_data.p_eff = p_eff;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
79 music_data.nfft = opts.nfft;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
80 music_data.Fs = opts.Fs;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
81 music_data.EVFlag = opts.EVFlag;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
82 music_data.range = opts.range;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
83
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
84
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
85 %--------------------------------------------------------------------------------------
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
86 function [options,msg] = music_options(xIsReal,p,varargin)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
87 %MUSIC_OPTIONS Parse the optional inputs to the MUSIC function.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
88 % MUSIC_OPTIONS returns a structure, OPTIONS, with the following fields:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
89 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
90 % options.nfft - number of freq. points at which the psd is estimated
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
91 % options.Fs - sampling freq. if any
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
92 % options.range - 'onesided' or 'twosided' pseudospectrum (they correspond to
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
93 % 'half' and 'whole' respectively, but are returned as is by
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
94 % psdoptions.m
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
95 % options.nw - number of columns in the data matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
96 % options.noverlap - number of samples to overlap
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
97 % options.window - a vector with window coefficients
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
98 % options.CorrFlag - a flag indicating whether the input is a correlation matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
99 % options.EVFlag - flag, 0 = MUSIC method ; 1 = EigenVector method
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
100 % options.CorrMatrOrd - order of the correlation matrix to be used in computations
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
101
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
102
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
103
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
104 % Assign Defaults
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
105 msg = '';
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
106 options.nw = [];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
107 options.noverlap = [];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
108 options.window = [];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
109 options.nfft = 256;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
110 options.Fs = [];
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
111 options.CorrFlag = 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
112 options.EVFlag = 0;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
113 % Determine if frequency vector specified
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
114 freqVecSpec = false;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
115 if (length(varargin) > 0 && isnumeric(varargin{1}) && length(varargin{1}) > 1)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
116 freqVecSpec = true;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
117 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
118
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
119 if xIsReal && ~freqVecSpec,
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
120 options.range = 'onesided';
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
121 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
122 options.range = 'twosided';
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
123 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
124
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
125 [options,msg] = psdoptions(xIsReal,options,varargin{:});
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
126
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
127 if length(options.nfft) > 1,
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
128 if strcmpi(options.range,'onesided')
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
129 warning(generatemsgid('InconsistentRangeOption'),...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
130 'Ignoring the ''onesided'' option. When a frequency vector is specified, a ''twosided'' PSD is computed');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
131 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
132 options.range = 'twosided';
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
133 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
134
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
135 % psdoptions doesn't handle this field, assign it separetely
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
136 options.CorrMatrOrd = 2*p(1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
137
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
138 %-----------------------------------------------------------------------------------------
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
139 function [eigenvals,eigenvects] = computeeig(x,CorrFlag,CorrMatrOrd,nw,noverlap,window,EVFlag)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
140 %COMPUTEEIG Compute eigenvalues and eigenvectors of correlation matrix.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
141 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
142 % Inputs:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
143 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
144 % x - input vector or matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
145 % CorrFlag - (flag) indicates whether x is a correlation matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
146 % nw - (integer) length of the rows of the data matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
147 % (only used if x is vector)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
148 % noverlap - (integer) overlap between the rows of the data matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
149 % (used in conjunction with nw)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
150 % window - (vector) window to be applied to each column of data
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
151 % matrix (not used if x is a correlation matrix)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
152 % EVFlag - True if eigenvector method, false if MUSIC.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
153 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
154 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
155 % Outputs:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
156 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
157 % eigenvals
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
158 % eigenvects
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
159 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
160 % If x is a matrix,
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
161 % If CorrFlag = 1, input x is a correlation matrix, we compute the
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
162 % eigendecomposition and order the eigenvalues and eigenvectors.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
163 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
164 % If x is a vector,
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
165 % a data matrix is formed by calling corrmtx unless a custom nw
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
166 % and noverlap are specified. In that case, we use buffer to form
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
167 % the data matrix.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
168 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
169 % If window is not empty, each row of the data matrix will be
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
170 % multiplied by the window.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
171
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
172 % Determine if the input is a matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
173 xIsMatrix = ~any(size(x)==1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
174
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
175 if xIsMatrix && CorrFlag,
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
176 % Input is Correlation matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
177
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
178 % Compute the eigenvectors and eigenvalues
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
179 [E,D] = eig((x+x')/2); % Ensure hermitian
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
180 [eigenvals,indx] = sort(diag(D),'descend');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
181 eigenvects = E(:,indx);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
182
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
183 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
184 if xIsMatrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
185 % Input is already a data matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
186 [Mx,Nx] = size(x); % Determine size of data matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
187 if EVFlag && (Nx > Mx),
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
188 errmsg = 'The number of columns in the data matrix cannot exceed the number of rows.';
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
189 error(generatemsgid('invalidDataMatrix'),errmsg);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
190 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
191
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
192 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
193 % x is a vector
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
194 x = x(:); % Make it a column
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
195 if isempty(nw),
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
196 x = corrmtx(x,CorrMatrOrd-1,'cov');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
197 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
198 if EVFlag && nw > (ceil((length(x)-nw)/(nw-noverlap))+1),
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
199 errmsg = sprintf(['The segment length and overlap specified result in\n',...
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
200 'a data matrix with more columns than rows.']);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
201 error(generatemsgid('invalidDataMatrix'),errmsg);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
202 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
203 Lx = length(x);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
204 x = buffer(x,nw,noverlap,'nodelay');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
205 if Lx <= nw,
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
206 error(generatemsgid('invalidSegmentLength'),'The segment length, NW, must be smaller than the signal length.');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
207 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
208 x = x'./sqrt(Lx-nw); % Scale appropriately such that X'*X is a scaled estimate of R
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
209 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
210 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
211 if ~isempty(window),
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
212 % Apply window to each row of data matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
213 if length(window) ~= size(x,2),
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
214 error(generatemsgid('InvalidDimensions'),'Window length must equal the number of columns in the data matrix.');
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
215 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
216 window = repmat(window(:).',size(x,1),1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
217 x = x.*window;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
218 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
219
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
220 % Compute the eigenvectors and eigenvalues via the SVD
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
221 [U,S,eigenvects] = svd(x,0);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
222 eigenvals = diag(S).^2; % We need to square the singular values here
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
223 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
224
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
225
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
226 %--------------------------------------------------------------------------------------------
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
227 function p_eff = determine_signal_space(p,eigenvals)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
228 %DETERMINE_SIGNAL_SPACE Determines the effective dimension of the signal subspace.
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
229 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
230 % Inputs:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
231 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
232 % p - (scalar or vector) signal subspace dimension
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
233 % (but may contain a desired threshold).
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
234 % eigenvals - (vector) contains the eigenvalues (sorted in decreasing order)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
235 % of the correlation matrix
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
236 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
237 % Outputs:
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
238 %
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
239 % p_eff - The effective dimension of the signal subspace. If a threshold
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
240 % is given as p(2), the signal subspace will be equal to the number
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
241 % of eigenvalues, NEIG, greater than the threshold times the smallest
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
242 % eigenvalue. However, the dimension of the signal subspace is at most
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
243 % p(1), so that if NEIG is greater than p(1), p_eff will be equal to
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
244 % p(1). If the threshold criteria results in an empty signal subspace,
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
245 % once again we make p_eff = p(1).
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
246
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
247
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
248 % Use the signal space dimension or the threshold to separate the noise subspace eigenvectors
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
249 if length(p) == 2,
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
250 % The threshold will be the input threshold times the smallest eigenvalue
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
251 thresh = p(2)*eigenvals(end);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
252 indx = find(eigenvals > thresh);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
253 if ~isempty(indx)
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
254 p_eff = min( p(1), length(indx) );
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
255 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
256 p_eff = p(1);
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
257 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
258 else
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
259 p_eff = p;
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
260 end
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
261
f0afece42f48 Import.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
262 % [EOF] - music.m