time_frec_distribution.m

主成分分析和偏最小二乘SquaresPrincipal成分分析（ PCA ）和偏最小二乘（ PLS ）

%  Usage
%    [tfd, t, f] = rihaczek2(x, fs, nfreq)
%
%  Inputs
%    x     signal vector
%    fs    sampling frequency of x (optional, default is 1 sample/second)
%    nfreq number of samples to compute in frequency (optional, default
%          is the length of x)
%
%  Outputs
%    tfd  matrix containing the Rihaczek distribution of signal x.  If x has
%         length N, then tfd will be nfreq by N. (optional)
%    t    vector of sampling times (optional)
%    f    vector of frequency values (optional)

% Copyright (C) -- see DiscreteTFDs/Copyright

error(nargchk(1, 3, nargin));
if (nargin == 1)
  [tfd, t, f] = rihaczek1(x);
elseif (nargin == 2)
  [tfd, t, f] = rihaczek1(x, fs);
elseif (nargin == 3)
  [tfd, t, f] = rihaczek1(x, fs, nfreq);
end
 

function [tfd, t, f] = spec2(x, fs, nfreq, decf, w, how)
% spec2 -- Compute samples of the type II spectrogram.
%
%  Usage
%    [tfd, t, f] = spec2(x, fs, nfreq, decf, w, how)
%
%  Inputs
%    x      signal vector
%    fs     sampling frequency of x (optional, default is 1 sample/second)
%    nfreq  number of samples to compute in frequency (optional, default
%           is 256)
%    decf   sub-sampling factor in time of the stft (optional, default 
%           is 1, i.e. no sub-sampling)
%    w      if length(w)==1 then the window is a guassian with duration 'w'
%           (see chirplets.m), otherwise 'w' is the window.  'w' must have an
%           odd length to have a center point. (optional, default is a 
%           gaussian with a duration of 5)
%    how    if how='short' then the stft is computed for the same times
%           as the signal and 
%               size(tfd,2) = length(x)
%           if how='long' then the stft will be computed for times before
%           and after the times of the signal and
%               size(tfd,2) = length(x) + length(w) - 1
%           (optional, default is 'short')
%
%  Outputs
%    tfd  matrix containing the spectrogram of signal x (optional)
%    t    vector of sampling times (optional)
%    f    vector of frequency values (optional)

% Copyright (C) -- see DiscreteTFDs/Copyright

error(nargchk(1, 6, nargin));
if (nargin == 1)
  [tfd, t, f] = stft2(x);
elseif (nargin == 2)
  [tfd, t, f] = stft2(x, fs);
elseif (nargin == 3)
  [tfd, t, f] = stft2(x, fs, nfreq);
elseif (nargin == 4)
  [tfd, t, f] = stft2(x, fs, nfreq, decf);
elseif (nargin == 5)
  [tfd, t, f] = stft2(x, fs, nfreq, decf, w);
elseif (nargin == 6)
  [tfd, t, f] = stft2(x, fs, nfreq, decf, w, how);
end

tfd = real(tfd.*conj(tfd));


function [tfd, t, f] = stft2(x, fs, nfreq, decf, w, how)
% stft2 -- Compute samples of the type II short-time Fourier transform.
%
%  Usage
%    [tfd, t, f] = stft2(x, fs, nfreq, decf, w, how)
%
%  Inputs
%    x      signal vector
%    fs     sampling frequency of x (optional, default is 1 sample/second)
%    nfreq  number of samples to compute in frequency (optional, default
%           is 256)
%    decf   sub-sampling factor in time of the stft (optional, default 
%           is 1, i.e. no sub-sampling)
%    w      if length(w)==1 then the window is a guassian with duration 'w'
%           (see chirplets.m), otherwise 'w' is the window.  'w' must have an
%           odd length to have a center point. (optional, default is a 
%           gaussian with a duration of 5)
%    how    if how='short' then the stft is computed for the same times
%           as the signal and 
%               size(tfd,2) = length(x)
%           if how='long' then the stft will be computed for times before
%           and after the times of the signal and
%               size(tfd,2) = length(x) + length(w) - 1
%           (optional, default is 'short')
%
%  Outputs
%    tfd   matrix containing the STFT of signal x
%    t     vector of sampling times (optional)
%    f     vector of frequency values (optional)

% Copyright (C) -- see DiscreteTFDs/Copyright

x = x(:);
N=length(x);

error(nargchk(1, 6, nargin));

if (nargin < 2) fs = 1; end
if (nargin < 3)
  M = 63;
  w = chirplets(M,[1 (M+1)/2 0 0 5]);
  nfreq = 256;
end
if (nargin < 4) decf = 1; end
if (nargin < 5)
  M = 63;
  w = chirplets(M,[1 (M+1)/2 0 0 5]);
end
if (nargin < 6) how = 'short'; end

if (length(w)==1)
  M = 2*round(32*w/5)-1;
  w = chirplets(M,[1 (M+1)/2 0 0 w]);
else
  w = w(:);
  M = length(w);
end

if (round(M/2) == M/2)
  error('window length must be odd')
end
if (N <= M)
  error('the signal must be longer than the window')
end

% Create a matrix filled with signal values.

x = [x; zeros(decf,1)];   % prevents going past the end of the array

if (how(1) == 'l')
  L = ceil((M+N-1)/decf);
  tfd=zeros(M,L);
  tfd(M,1)=x(1);
  for n=2:L,
    tfd(1:M-decf,n)=tfd(1+decf:M,n-1);
    t = (n-1)*decf+1;      % real time
    if (t-decf+1 <= N)    % we still have new data to add
      tfd(M-decf+1:M,n) = x(t-decf+1:t);
    end
  end
else % how == 'short'
  L = ceil(N/decf);
  tfd=zeros(M,L);
  tfd((M+1)/2:M,1)=x(1:(M+1)/2);
  for n=2:L,
    tfd(1:M-decf,n)=tfd(1+decf:M,n-1);
    t = (n-1)*decf+1;      % real time
    if ((M-1)/2+t-decf+1 <= N)
      tfd(M-decf+1:M,n)=x((M-1)/2+t-decf+1:(M-1)/2+t);
    end
  end
end

% Window the data.
w=w*ones(1,L);
tfd=tfd.*w;

% take care of the case if M > nfreq
if (M>nfreq)
  P = ceil(M/nfreq);
  tfd = [tfd ; zeros(P*nfreq-M,L)];
  tfd = reshape(tfd,nfreq,P,L);
  tfd = squeeze(sum(tfd,2));
end

% Perform the column ffts to get the stft.
tfd = fft(tfd, nfreq);
tfd = tfdshift(tfd)/sqrt(nfreq);

if (how(1) == 'l')
  t = 1/fs * ((0:N+M-2) - (M-1)/2);
  t = t(1:decf:end);
else
  t = 1/fs * (0:N-1);
  t = t(1:decf:end);
end
f = -fs/2:fs/nfreq:fs/2;
f = f(1:nfreq); 

function ptfd(tfd, t, f, fs)
% ptfd -- Display an image plot of a TFD with a linear amplitude scale.
%
%  Usage
%    ptfd(tfd, t, f, fs)
%
%  Inputs
%    tfd  time-frequency distribution
%    t    vector of sampling times (optional)
%    f    vector of frequency values (optional)
%    fs   font size of axis labels (optional)

% Copyright (C) -- see DiscreteTFDs/Copyright

error(nargchk(1, 4, nargin));

if (nargin < 4)
  fs = 10;
end
if (nargin < 3)
  f = [-0.5 0.5];
end
if (nargin < 2)
  t = [1 size(tfd,2)];
end

if isempty(t)
  t = [1 size(tfd,2)];
end
if isempty(f)
  f = [-0.5 0.5];
end

imagesc(t, f, abs(tfd)), axis('xy'), xlabel('time','FontSize',fs), ylabel('frequency','FontSize',fs), set(gca,'FontSize',fs); 

function ptfddb(tfd, dbs, t, f, fs)
% ptfddb -- Display an image plot of a TFD with a dB amplitude scale.
%
%  Usage
%    ptfddb(tfd, dbs, t, f, fs)
%
%  Inputs
%    tfd  time-frequency distribution
%    dbs  rabge in dBs (optional, default is 25)
%    t    vector of sampling times (optional)
%    f    vector of frequency values (optional)
%    fs   font size of axis labels (optional)

% Copyright (C) -- see DiscreteTFDs/Copyright

error(nargchk(1, 5, nargin));

if (nargin < 5)
  fs = 10;
end
if (nargin < 4)
  f = [-0.5 0.5];
end
if (nargin < 3)
  t = [1 size(tfd,2)];
end
if (nargin < 2)
  dbs = 25;
end

if isempty(t)
  t = [1 size(tfd,2)];
end
if isempty(f)
  f = [-0.5 0.5];
end

tfd=tfd./sum(sum(tfd));
tfd=20*log10(abs(tfd)+eps);
tfd = tfd - max(max(tfd));
imagesc(t, f, tfd), axis('xy'), xlabel('time','FontSize',fs), ylabel('frequency','FontSize',fs),set(gca,'FontSize',fs); 

function [tfd, t, f] = rihaczek1(x, fs, nfreq);
% rihaczek1 -- Compute samples of the type I Rihaczek distribution.
%
%  Usage
%    [tfd, t, f] = rihaczek1(x, fs, nfreq)
%
%  Inputs
%    x     signal vector
%    fs    sampling frequency of x (optional, default is 1 sample/second)
%    nfreq number of samples to compute in frequency (optional, default
%          is the length of x)
%
%  Outputs
%    tfd  matrix containing the Rihaczek distribution of signal x.  If x has
%         length N, then tfd will be nfreq by N. (optional)
%    t    vector of sampling times (optional)
%    f    vector of frequency values (optional)

% Copyright (C) -- see DiscreteTFDs/Copyright

% specify defaults
x = x(:);
N = length(x);

error(nargchk(1, 3, nargin));
if (nargin < 3)
  nfreq = N;
end
if (nargin < 2)
  fs = 1;
end
  
X = fft(x,nfreq);
tfd = (conj(X) * x.') .* exp(-j*[0:2*pi/nfreq:2*pi*(1-1/nfreq)]'*[0:N-1]);
tfd = tfdshift(tfd)/nfreq;

t = 1/fs * (0:N-1);
f = -fs/2:fs/nfreq:fs/2;
f = f(1:nfreq);

function out = tfdshift(in)
% tfdshift -- Shift the spectrum of a TFD by pi radians.
%
%  Usage
%    out = tfdshift(in)
%
%  Inputs
%    in   time-frequency distribution
%
%  Outputs
%    out  shifted time-frequency distribution

% Copyright (C) -- see DiscreteTFDs/Copyright

error(nargchk(1, 1, nargin));

N = size(in, 1);
M = ceil(N/2);
out = [in(M+1:N,:) ; in(1:M,:)]; 

function x = chirplets(N, P)
% chirplets -- make a signal that is a sum of chirplets
% 
%  Usage
%    x = chirplets(N, P)
%
%  Inputs
%    N    length of signal
%    P    matrix of parameters [amp time freq chirp_rate duration; ...]
%         (optional, default is [1 N/2+1 0 0 sqrt(N/4/pi)])
%
%  Outputs
%    x    signal
%
% d is the standard deviation of the guassian, d=sqrt(N/4/pi) gives an 
% atom with a circular Wigner distribution, and 2*sqrt(2)*d is the 
% Rayleigh limit.
%
%  Examples
%    N=128; x=chirplets(N, [1 N/2+1 0 0 sqrt(N/4/pi)]);
%    x=chirplets(128, [1 55 0 2*pi/128 12 ; 1 75 0 2*pi/128 12]);

% Copyright (C) -- see DiscreteTFDs/Copyright

error(nargchk(1,2,nargin));

if (nargin < 2)
  if (rem(N,2)==0)
    center = N/2+1;
  else
    center = (N+1)/2;
  end
  P = [1 center 0 0 sqrt(N/4/pi)];
end

if (size(P,2) ~= 5)
  error('Matrix P has the wrong number of columns.')
end

x = zeros(N,1);
n = (1:N)';
for p = 1:size(P,1),
  A = P(p,1);
  t = P(p,2);
  f = P(p,3);
  cr = P(p,4);
  d = P(p,5);
  am = exp(-((n-t)/2/d).^2) * sqrt(1/sqrt(2*pi)/d);
  x = x + A * am .* exp(j * (cr/2*(n-t).^2 + f*(n-t)));
end