Code covered by the BSD License  

Highlights from
Principal Component Spectral Analysis

image thumbnail

Principal Component Spectral Analysis

by

 

09 Apr 2009 (Updated )

A two-dimensional histogram of Power Spectral Densities.

pcsa(X, paxis, faxis, varargin)
function varargout = pcsa(X, paxis, faxis, varargin) 
%PCSA Principal Component Spectral Analysis.
%   N = PCSA(X, F, P) returns the PCSA matrix of the signal specified by
%   vector X in the matrix N. P contains psd value bins (in decibels) and F
%   contains the frequency bins used in forming the matrix N. PCSA is
%   calculated from the spectrogam as a two-variable histogram.
%
%   N = PCSA(X, [], []) uses the default F and P values.
%
%   [N,F,P] = PCSA(X, [], []) returns vectors P and F, where P contains 
%   the magnitude bins and F contains the frequency bins where the PCSA values 
%   are calculated. 
%
%   PCSA uses SPECTROGRAM to determine the frequencies in the signal X. 
%   Parameters for the SPECTROGRAM function can be added to the end of the 
%   PCSA call.
%
%   PCSA(...) with no output arguments plots the PCSA.
%
%   EXAMPLE 1: Display the PCSA of each segment of a sum of two 
%   linear chirps.
%     T = 0:0.0001:2;                  % 2 secs @ 10kHz sample rate
%     F = 0:0.5:100;
%     X = chirp(T,15,2,35)*0.3 + chirp(T,75,2,25);    
%     subplot(2,1,1)
%     pcsa(X,[],[],2048,1950,F,1E4);    % Display the PCSA
%     title('PCSA')
%     % Compare with spectrogram
%     subplot(2,1,2)
%     spectrogram(X,2048,1950,F,1E4, 'yaxis'); colorbar
%     title('Spectrogram')
%
%   See also SPECTROGRAM.

%   Author: M. Kesaniemi
%   Copyright 2009 The MathWorks, Inc.


error(nargchk(3,8,nargin,'struct'));
error(nargoutchk(0,3,nargout,'struct'));

[Y,F,T,P] = spectrogram(X, varargin{:}, 'yaxis');

P = 10*log10(abs(P));

if isempty(faxis)
    faxis = F(2:3:end);    
end;
if isempty(paxis)
    D = abs(max(P(:))-min(P(:)));
    paxis = linspace(min(P(:))+2*D/3, max(P(:))+D/20, ...
        max([numel(P)/(20*numel(faxis)), 8]));    
end;

N = hist3(cat(2, reshape(P, [numel(P) 1]), ...
    repmat(F, size(P, 2), 1)), ...
    'Edges', {paxis faxis});

switch nargout,
    case 0,
        surf(faxis, paxis, N, 'EdgeColor','none');
        axis xy; axis tight; colormap(jet); view(0,90);
        ylabel('Power');
        xlabel('Frequency (Hz)');
        colorbar;        
    case 1,
        varargout = {N};
    case 2,
        varargout = {N,faxis};
    case 3,
        varargout = {N,faxis,paxis};
end

Contact us