diff options
Diffstat (limited to 'macros/peig.sci')
| -rw-r--r-- | macros/peig.sci | 101 |
1 files changed, 58 insertions, 43 deletions
diff --git a/macros/peig.sci b/macros/peig.sci index 7529edc..b4157d7 100644 --- a/macros/peig.sci +++ b/macros/peig.sci @@ -14,47 +14,62 @@ function [S,f,v,e] = peig(varargin) // // Parameters: // x - int|double - vector|matrix - // Input signal. In case of a matrix, each row of x represents a - // seperate observation of the signal. If 'corr' flag is specified, + // Input signal. In case of a matrix, each row of x represents a + // seperate observation of the signal. If 'corr' flag is specified, // then x is the correlation matrix. - // If w is not specified in the input, it is determined by the - // algorithm. If x is real valued, then range of w is [0, pi]. + // If w is not specified in the input, it is determined by the + // algorithm. If x is real valued, then range of w is [0, pi]. // Otherwise, the range of w is [0, 2pi) // p - int|double - scalar|vector // p(1) is the dimension of the signal subspace - // p(2), if specified, represents a threshold that is multiplied by + // p(2), if specified, represents a threshold that is multiplied by // the smallest estimated eigenvalue of the signal's correlation matrix. // w - int|double - vector - // w is the vector of normalized frequencies over which the + // w is the vector of normalized frequencies over which the // pseuspectrogram is to be computed. // nfft - int - scalar (Default = 256) // Length of the fft used to compute pseudospectrum. The length of S // (and hence w/f) depends on the type of values in x and nfft. - // If x is real, length of s is (nfft/2 + 1) {Range of w = [0, pi]} if + // If x is real, length of s is (nfft/2 + 1) {Range of w = [0, pi]} if // nfft is even and (nfft+1)/2 {Range of w = [0, pi)} otherwise. // If x is complex, length of s is nfft. // fs - int|double - scalar (Default = 1) - // Sampling rate. Used to convert the normalized frequencies (w) to + // Sampling rate. Used to convert the normalized frequencies (w) to // actual values (f) and vice-versa. // nwin - int|double - scalar (int only)|vector (Default = 2*p(1)) // If nwin is scalar, it is the length of the rectangular window. // Otherwise, the vector input is considered as the window coefficients. // Not used if 'corr' flag present. - // If x is a vector, windowing not done in nwin in scalar. If x is a - // matrix, + // If x is a vector, windowing not done in nwin in scalar. If x is a + // matrix, // noverlap - int - scalar (Default = nwin-1) - // number of points by which successive windows overlap. noverlap not + // number of points by which successive windows overlap. noverlap not // used if x is a matrix // freqrange - string - // The range of frequencies over which the pseudospetrogram is + // The range of frequencies over which the pseudospetrogram is // computed. Three possible values - 'onesided', 'twosided', 'centered' // 'corr' flag - // Presence indicates that the primary input x is actually a + // Presence indicates that the primary input x is actually a // correlation matrix // // Examples: - // TODO: - // +//1. +// fs = 100; +// t = 0:1/fs:1-1/fs; +// s = 2*sin(2*%pi*25*t)+sin(2*%pi*35*t)+rand(1,100,"normal"); +// [S,w]=peig(s,2,512,fs,'half'); +// plot(w,S); + //OUTPUT: gives the plot of power vs normalised frequencies +////2. +// fs = 100; +//t = 0:1/fs:1-1/fs; +//s = 2*sin(2*%pi*25*t)+sin(2*%pi*35*t); +//[S,w]=peig(s,2,512,fs,'half'); +//plot(w,S); + + //EXECUTE FUNCTIONS subspaceMethodsInputPars.sci, musicBase.sci PRIOR THE EXECUTION OF THIS FUNCTION + + // See also // rooteig | pmusic | pmtm | pcov | pmcov | pburg | pyulear | pwelch | corrmtx // @@ -64,16 +79,16 @@ function [S,f,v,e] = peig(varargin) // References // [1] Petre Stoica and Randolph Moses, Introduction To Spectral // Analysis, Prentice-Hall, 1997, pg. 15 - // [2] S. J. Orfanidis, Optimum Signal Processing. An Introduction. + // [2] S. J. Orfanidis, Optimum Signal Processing. An Introduction. // 2nd Ed., Macmillan, 1988. - + funcprot(0); - - exec('subspaceMethodsInputParser.sci',-1); - exec('musicBase.sci',-1); - [data, msg, err_num] = subspaceMethodsInputParser(varargin); - + //exec('/home/shashi/Desktop/FOSSEE-Signal-Processing-Toolbox-master/macros/subspaceMethodsInputParser.sci',-1); + // exec('/home/shashi/Desktop/FOSSEE-Signal-Processing-Toolbox-master/macros/musicBase.sci',-1); + + [data, msg, err_num] = subspaceMethodsInputPars(varargin); + if length(msg)==0 then // no error occured else @@ -90,7 +105,7 @@ function [S,f,v,e] = peig(varargin) // computing the pseudospectrum [S,f] = pseudospectrum(musicData.noiseEigenvects, ... musicData.eigenvals,data.w,data.nfft, data.fs, data.freqrange,data.isFsSpecified); - + v = musicData.noiseEigenvects; e = musicData.eigenvals; @@ -100,15 +115,15 @@ function [pspec,w] = pseudospectrum(noiseEigenvects, eigenvals, freqvector, ... nfft, fs, freqrange,isFsSpecified) // disp("noise eigenvects in pseudospectrum - "); // disp(noiseEigenvects); - + // NOTE: the only difference from the pmusic code // Using the noise subspace eigenvalues as weights - weights = eigenvals(($-size(noiseEigenvects,2)+1):$); + weights = eigenvals(($-size(noiseEigenvects,2)+1):$); denominator = 0; - + isFreqGiven = %F; - + for i=1:size(noiseEigenvects,2); // disp("looping in pseudospectrum"); @@ -129,7 +144,7 @@ function [pspec,w] = pseudospectrum(noiseEigenvects, eigenvals, freqvector, ... pspec = pspec(:); - + if ~isFreqGiven then // correcting the range of pspec according to the user specification if strcmpi(freqrange, 'onesided')==0 then @@ -141,11 +156,11 @@ function [pspec,w] = pseudospectrum(noiseEigenvects, eigenvals, freqvector, ... end pspec = pspec(range); w = w(range); - + elseif strcmpi(freqrange,'centered')==0 then // convert two sided spectrum to centered rem = modulo(nfft,2); - + if rem then idx = [((nfft+1)/2+1):nfft 1:(nfft+1)/2]; else @@ -153,7 +168,7 @@ function [pspec,w] = pseudospectrum(noiseEigenvects, eigenvals, freqvector, ... end pspec = pspec(idx); w = w(idx); - + if rem then w(1:(nfft-1)/2) = - w(nfft:-1:((nfft+1)/2+1)); else @@ -165,8 +180,8 @@ function [pspec,w] = pseudospectrum(noiseEigenvects, eigenvals, freqvector, ... endfunction function [h,w] = computeFreqResponseByFFT(b,n,fs,isFsSpecified) - // returns the frequency response (h) and the corresponding frequency - // values (w) for a digital filter with numerator b. The evaluation of the + // returns the frequency response (h) and the corresponding frequency + // values (w) for a digital filter with numerator b. The evaluation of the // frequency response is done at n points in [0,fs) using fft algorithm // // Similar to MATLAB's freqz(b,a,n,'whole',fs) @@ -176,10 +191,10 @@ function [h,w] = computeFreqResponseByFFT(b,n,fs,isFsSpecified) w = linspace(0,2*%pi,n+1)'; w($) = []; w(1) = 0; // forcing the first frequency to be 0 - + // forcing b and a to be column vectors b = b(:); - + // zero padding for fft zeroPadLength = n - length(b); zeroPad = zeros(zeroPadLength,1); @@ -187,34 +202,34 @@ function [h,w] = computeFreqResponseByFFT(b,n,fs,isFsSpecified) h = fft(b); - + if isFsSpecified then w = w*fs/(2*%pi); end - + endfunction function [h,w] = computeFreqResponseByPolyEval(b,f,fs,isFsSpecified) // returns the frequency response (h) for a digital filter with numerator b. // The evaluation of the frequency response is done at frequency values f - + // disp(f); // disp(isFsSpecified); - + f = f(:); b = b(:); - + n = length(b); powerMatrix = zeros(length(f),n); powerMatrix(:,1) = 1; for i=2:n powerMatrix(:,i) = exp(f*(-i+1)*%i); end - + h = powerMatrix*b; - + if isFsSpecified then w = f * fs/(2*%pi); end - + endfunction |