Code covered by the BSD License

# Principal Component Spectral Analysis

by

### Martti K (view profile)

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')
%

%   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

```