www.gusucode.com > Matlab动力系统和时间序列分析工具箱 > Matlab动力系统和时间序列分析工具箱/lab432/toolbox/time_frec_distribution.m
function varargout = time_frec_distribution(varargin)
% time_frec_distribution creates time-frecuency distribution of time series
% [tfd,t,f]=time_frec_distribution(x,giu_mode,method,param);
% Parameters:
% x - data vector or structure with fields
% .data - data vector
% .time - time vector
% gui_mode - 'gui' or 'silent' (no gui)
% method -
% 1 - Born-Jordan distribution
% 2 - Binomial distribution
% 3 - quasi-Wigner distribution
% 4 - Levin distribution
% 5 - Rihaczek distribution
% 6 - Page distribution
% 7 - Margenau-Hill distribution
% 8 - Spectrogram
% 9 - Short-time Fourier Transform
% For methods 1-3 param is a length of rectangular lag window on the autocorr. function (even)
% For methods 4-7 no parameters required
% For methods 8-9 param is a duration of gaussian window (odd)
%
% Based on DiscreteTFD
% Jeffrey O'Neill <jeffo@bu.edu>
% 8 St. Mary's Street
% Boston, MA 02135
% Last Modified by GUIDE v2.5 22-Oct-2004 17:03:50
% last modified 8.01.04
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @time_frec_distribution_OpeningFcn, ...
'gui_OutputFcn', @time_frec_distribution_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before time_frec_distribution is made visible.
function time_frec_distribution_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to time_frec_distribution (see VARARGIN)
% Choose default command line output for time_frec_distribution
global tfdDATA
handles.output = hObject;
guidata(hObject, handles);
set(hObject,'name','Time-frequency distribution ');
set([handles.t_fsw handles.e_l_fsw],'visible','off');
if isempty(varargin)
error('Not enought input parameters');
end
if isa(varargin{1},'struct')
tfdDATA.data=varargin{1};
else
tfdDATA.data.data=varargin{1};
tfdDATA.data.time=1:length(varargin{1});
end
if rem(length(tfdDATA.data.data),2)==0
set(handles.e_l_tsw,'string',num2str(length(tfdDATA.data.data)));
else
set(handles.e_l_tsw,'string',num2str(length(tfdDATA.data.data)-1));
end
if length(varargin)>1
if strcmp(varargin{2},'gui')|strcmp(varargin{2},'silent')
tfdDATA.gui_mode=varargin{2};
else
warning('Wrong value of ''gui_mode'' parameter. ''gui_mode'' set to ''gui''');
tfdDATA.gui_mode='gui';
end
else
tfdDATA.gui_mode='gui';
end
tfdDATA.autocalc=0;
if length(varargin)>2
set(handles.method,'value',varargin{3});
end
if length(varargin)>=3
switch varargin{3} % method
case {1,2,3}
set([handles.t_fsw handles.e_l_fsw],'visible','off');
set([handles.t_tsw handles.t_l_tsw handles.e_l_tsw ],'visible','on');
if length(varargin)>3
tfdDATA.autocalc=1;
set(handles.e_l_tsw,'string',num2str(varargin{4}));
end
case {4,5,6,7}
set([handles.t_fsw handles.e_l_fsw],'visible','off');
set([handles.t_tsw handles.t_l_tsw handles.e_l_tsw ],'visible','off');
tfdDATA.autocalc=1;
case {8,9}
set([handles.t_fsw handles.e_l_fsw],'visible','on');
set([handles.t_tsw handles.t_l_tsw handles.e_l_tsw ],'visible','off');
if length(varargin)>3
tfdDATA.autocalc=1;
set(handles.e_l_fsw,'string',num2str(varargin{4}));
end
end
end
% UIWAIT makes time_frec_distribution wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = time_frec_distribution_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% handles structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
global tfdDATA tfdOUT
if tfdDATA.autocalc
compute_Callback(handles.compute, [], handles);
varargout{1} =tfdOUT.tfd;
varargout{2} =tfdOUT.t;
varargout{3} =tfdOUT.f;
if strcmp(tfdDATA.gui_mode,'silent')
try
close(handles.figure);
catch
end
end
else
uiwait;
try
varargout{1} =tfdOUT.tfd;
varargout{2} =tfdOUT.t;
varargout{3} =tfdOUT.f;
catch
end
end
% --- Executes on button press in cancel.
function cancel_Callback(hObject, eventdata, handles)
% hObject handle to cancel (see GCBO)
% handles structure with handles and user data (see GUIDATA)
close(gcf);
% --- Executes on button press in compute.
function compute_Callback(hObject, eventdata, handles)
% hObject handle to compute (see GCBO)
% handles structure with handles and user data (see GUIDATA)
global tfdDATA tfdOUT
I=get(handles.method,'value');
w=str2num(get(handles.e_l_fsw,'string'));
wlen=str2num(get(handles.e_l_tsw,'string'));
Dt=diff(tfdDATA.data.time);
if max(Dt-mean(Dt))>1e-10
warndlg('Time instants is not equally spaced. Frequency scale is wrong.','Warning');
end
fs=1/Dt(1);
nfreq=length(tfdDATA.data.data);
switch I
case 1
[tfd, t, f] = born_jordan2(tfdDATA.data.data, fs, nfreq, wlen);
case 2
[tfd, t, f] = binomial2(tfdDATA.data.data, fs, nfreq, wlen);
case 3
[tfd, t, f] = qwigner2(tfdDATA.data.data, fs, nfreq, wlen);
case 4
[tfd, t, f] = levin2(tfdDATA.data.data, fs, nfreq);
case 5
[tfd, t, f] = rihaczek2(tfdDATA.data.data, fs, nfreq);
case 6
[tfd, t, f] = page2(tfdDATA.data.data, fs, nfreq);
case 7
[tfd, t, f] = margenau_hill2(tfdDATA.data.data, fs, nfreq);
case 8
[tfd,t,f] = spec2(tfdDATA.data.data,fs,nfreq,1,w,'short');
case 9
[tfd,t,f] = stft2(tfdDATA.data.data,fs,nfreq,1,w,'short');
end
tfdOUT.tfd=tfd;
tfdOUT.t=t;
tfdOUT.f=f;
if strcmp(tfdDATA.gui_mode, 'gui')
axes(handles.ax);
colormap('jet');
if get(handles.las,'Value')
ptfd(tfd, t, f, 8);
else
ptfddb(tfd, 25,t, f, 8)
end
end
% --- Executes during object creation, after setting all properties.
function e_l_fsw_CreateFcn(hObject, eventdata, handles)
% hObject handle to e_l_fsw (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject,'string',num2str(5));
function e_l_fsw_Callback(hObject, eventdata, handles)
% hObject handle to e_l_fsw (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of e_l_fsw as text
% str2double(get(hObject,'String')) returns contents of e_l_fsw as a double
% --- Executes during object creation, after setting all properties.
function e_l_tsw_CreateFcn(hObject, eventdata, handles)
% hObject handle to e_l_tsw (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
% global tfdDATA
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
function e_l_tsw_Callback(hObject, eventdata, handles)
% hObject handle to e_l_tsw (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of e_l_tsw as text
% str2double(get(hObject,'String')) returns contents of e_l_tsw as a double
% --- Executes during object creation, after setting all properties.
function method_CreateFcn(hObject, eventdata, handles)
% hObject handle to method (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: popupmenu controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
str='Born-Jordan distribution'; % 1
str=[str '|Binomial distribution']; % 2
str=[str '|quasi-Wigner distribution']; % 3
str=[str '|Levin distribution']; % 4
str=[str '|Rihaczek distribution']; % 5
str=[str '|Page distribution']; % 6
str=[str '|Margenau-Hill distribution']; % 7
str=[str '|Spectrogram']; % 8
str=[str '|Short-time Fourier Transform']; % 9
set(hObject,'string',str);
% --- Executes on selection change in method.
function method_Callback(hObject, eventdata, handles)
% hObject handle to method (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: contents = get(hObject,'String') returns method contents as cell array
% contents{get(hObject,'Value')} returns selected item from method
I=get(hObject,'value');
switch I
case {1,2,3}
set([handles.t_fsw handles.e_l_fsw],'visible','off');
set([handles.t_tsw handles.t_l_tsw handles.e_l_tsw ],'visible','on');
case {4,5,6,7}
set([handles.t_fsw handles.e_l_fsw],'visible','off');
set([handles.t_tsw handles.t_l_tsw handles.e_l_tsw ],'visible','off');
case {8,9}
set([handles.t_fsw handles.e_l_fsw],'visible','on');
set([handles.t_tsw handles.t_l_tsw handles.e_l_tsw ],'visible','off');
end
% --- Executes on button press in das.
function das_Callback(hObject, eventdata, handles)
% hObject handle to das (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hint: get(hObject,'Value') returns toggle state of das
set(hObject,'Value',1);
set(handles.las,'value',~get(hObject,'Value'));
% --- Executes on button press in las.
function las_Callback(hObject, eventdata, handles)
% hObject handle to las (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hint: get(hObject,'Value') returns toggle state of las
set(hObject,'Value',1);
set(handles.das,'value',~get(hObject,'Value'));
%%% DiscreteTFD's functions (Type2,(3): discrete and aperiodic; continuous and periodic) %%%
% Copyright (C) 1997, 1998, 1999 Jeffrey C. O'Neill
% Copyright (C) 1998, 1999 Boston University
% Copyright (C) 1997, 1998 L'Ecole Normale Superieure de Lyon
% Copyright (C) 1997 The Regents of the University of Michigan
% All Rights Reserved
%
% Jeffrey O'Neill <jeffo@bu.edu>
% 8 St. Mary's Street
% Boston, MA 02135
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
function [tfd, t, f] = binomial2(x, fs, nfreq, wlen)
% binomial2 -- Compute samples of the type II binomial distribution.
%
% Usage
% [tfd, t, f] = binomial2(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 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, 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 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.
%
% 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