time_frec_distribution.m

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

end

if (nfreq < wlen)
  error('rectangular lag window must be less than or equal to data length!');
end
if (wlen > N)
  error('rectangular lag window must be less than or equal to the length of the signal!');
end
w = wlen/2;

% make the binomial kernel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ker = zeros(w);
ker(1,1) = 1;
for tau = 2:w,
  temp = ker(tau-1,1:tau-1);
  ker(tau,1:tau) = ([0 temp] + [temp 0])/2;
end

% Do the computations.
%%%%%%%%%%%%%%%%%%%%%%

% make the acf for positive tau
acf = lacf2(x, w);

% convolve with the kernel
acf2 = fft(acf.');
ker = [ker zeros(w,N-w)];
ker2 = fft(ker.');
gacf = ifft(acf2.*ker2);
gacf = gacf.';

% make the gacf for negative lags
gacf = [gacf ; zeros(nfreq-wlen+1,N) ; conj(flipud(gacf(2:w,:)))];

%compute the tfd
tfd = real(fft(gacf));
tfd = tfdshift(tfd)/nfreq;

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


function [tfd, t, f] = born_jordan2(x, fs, nfreq, wlen)
% born_jordan2 -- Compute samples of the type II Born_Jordan distribution.
%
%  Usage
%    [tfd, t, f] = born_jordan2(x, fs, nfreq, wlen)
%
%  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 twice the length of x)
%    wlen  length of the rectangular lag window on the auto-correlation
%          function, must be less than or equal to nfreq (optional, default
%          is twice the length of x)
%
%  Outputs
%    tfd  matrix containing the Born_Jordan 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 DiscreeTFDs/Copyright

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

error(nargchk(1, 4, nargin));
if (nargin < 4) 
  wlen = 2*N;
end
if (nargin < 3)
  nfreq = 2*N;
end
if (nargin < 2)
  fs = 1;
end

if (nfreq < wlen)
  error('rectangular lag window must be less than or equal to data length!');
end
if (wlen > N)
  error('rectangular lag window must be less than or equal to the length of the signal!');
end
w = wlen/2;

% make the born jordan kernel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ker = zeros(w);
for tau = 1:w,
  ker(tau,1:tau) = ones(1,tau)/tau;
end


% Do the computations.
%%%%%%%%%%%%%%%%%%%%%%

% make the acf for positive tau
acf = lacf2(x, w);

% convolve with the kernel
acf2 = fft(acf.');
ker = [ker zeros(w,N-w)];
ker2 = fft(ker.');
gacf = ifft(acf2.*ker2);
gacf = gacf.';

% make the gacf for negative lags
gacf = [gacf ; zeros(nfreq-wlen+1,N) ; conj(flipud(gacf(2:w,:)))];

%compute the tfd
tfd = real(fft(gacf));
tfd = tfdshift(tfd)/nfreq;

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


function lacf = lacf2(x, mlag)
% lacf2 -- Compute samples of the type II local acf.
%
%  Usage
%    lacf = lacf2(x, mlag)
%
%  Inputs
%    x     signal vector
%    mlag   maximum lag to compute.  must be <= length(x).
%           (optional, defaults to length(x))
%
%  Outputs
%    lacf  matrix containing the lacf of signal x.  If x has
%         length N, then lacf will be nfreq by N. (optional)
%
% This function has a tricky sampling scheme, so be careful if you use it.

% Copyright (C) -- see DiscreteTFDs/Copyright

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

error(nargchk(1, 2, nargin));
if (nargin < 2) 
  mlag = N;
end

if (mlag > N)
  error('mlag must be <= length(x)')
end

% make the acf for positive tau
lacf = zeros(mlag, N);
for t=1:N
  mtau = min(mlag, N-t+1);
  lacf(1:mtau, t) = conj(x(t))*x(t:t+mtau-1);
end
 
function [tfd, t, f] = levin2(x, fs, nfreq)
% levin2 -- Compute samples of the type II Levin distribution.
%
%  Usage
%    [tfd, t, f] = levin2(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 twice the length of x)
%
%  Outputs
%    tfd  matrix containing the binomial 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 = 2*N;
end
if (nargin < 2)
  fs = 1;
end

% make the acf for positive tau
gacf = zeros(N);
for t=1:N;
  gacf(1:t,t) = conj(x(t)) * x(t:-1:1);
end

% negative tau
gacf = [gacf ; zeros(nfreq-2*N+1,N) ; conj(flipud(gacf(2:N,:)))];

%compute the tfd
tfd = real(fft(gacf));
tfd = tfdshift(tfd)/nfreq;

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


function [tfd, t, f] = margenau_hill2(x, fs, nfreq);
% margenau_hill2 -- Compute samples of the type II Margenau-Hill 
%                    distribution.
%
%  Usage
%    [tfd, t, f] = margenau_hill2(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 Margenau-Hill 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

tfd = real(tfd);

 
function [tfd, t, f] = page2(x, fs, nfreq)
% page2 -- Compute samples of the type II Page distribution.
%
%  Usage
%    [tfd, t, f] = page2(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 twice the length of x)
%
%  Outputs
%    tfd  matrix containing the binomial 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 = 2*N;
end
if (nargin < 2)
  fs = 1;
end

% positive tau
acf = zeros(N);
for t=1:N;
  gacf(1:N-t+1,t) = x(t) * conj(x(t:N));
end

% negative tau
gacf = [gacf ; zeros(nfreq-2*N+1,N) ; conj(flipud(gacf(2:N,:)))];

%compute the tfd
tfd = real(fft(gacf));
tfd = tfdshift(tfd)/nfreq;

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


function [tfd, t, f] = qwigner2(x, fs, nfreq, wlen)
% qwigner2 -- Compute samples of the type II quasi-Wigner distribution.
%
%  Usage
%    [tfd, t, f] = qwigner2(x, fs, nfreq, wlen)
%
%  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 twice the length of x)
%    wlen  length of the rectangular lag window on the auto-correlation
%          function, must be less than or equal to nfreq (optional, default
%          is twice the length of x)
%
%  Outputs
%    tfd  matrix containing the quasi-Wigner 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)
%
% Note that the Wigner distribution does
% not exist for type II signals, but the quasi-Wigner distribution is 
% a reasonable approximation.

% Copyright (C) -- see DiscreteTFDs/Copyright


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

error(nargchk(1, 4, nargin));
if (nargin < 4) 
  wlen = 2*N;
end
if (nargin < 3)
  nfreq = 2*N;
end
if (nargin < 2)
  fs = 1;
end

if (nfreq < wlen)
  error('rectangular lag window must be less than or equal to data length!');
end
if (wlen > N)
  error('rectangular lag window must be less than or equal to the length of the signal!');
end
w = wlen/2;

% make the quasi-wigner kernel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ker = zeros(w);
for tau = 1:2:w,
  ker(tau,(tau+1)/2) = 1;
end
for tau = 2:2:w;
  ker(tau,tau/2) = .5;
  ker(tau,tau/2+1) = .5;
end

% Do the computations.
%%%%%%%%%%%%%%%%%%%%%%

% make the acf for positive tau
acf = lacf2(x, w);

% convolve with the kernel
acf2 = fft(acf.');
ker = [ker zeros(w,N-w)];
ker2 = fft(ker.');
gacf = ifft(acf2.*ker2);
gacf = gacf.';

% make the gacf for negative lags
gacf = [gacf ; zeros(nfreq-wlen+1,N) ; conj(flipud(gacf(2:w,:)))];

%compute the tfd
tfd = real(fft(gacf));
tfd = tfdshift(tfd)/nfreq;

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


function [tfd, t, f] = rihaczek2(x, fs, nfreq);
% rihaczek2 -- Compute samples of the type II Rihaczek distribution.
%